plot_joint_fit_2020_contours.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Plots.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Experiment/Dmsq32Constraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/FC/FCSurface.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/FitVars.h"
#include "Utilities/rootlogon.C"
#include "OscLib/IOscCalc.h"

#include "../joint_fit_2020_loader_tools.h"
#include "../joint_fit_2020_datarelease_tools.h"

#include "TCanvas.h"
#include "TMarker.h"
#include "TBox.h"
#include "TLatex.h"
#include "TColor.h"
#include "TGraph.h"
#include "TVectorD.h"
#include "TF1.h"
#include "TLegend.h"
#include "TLine.h"
#include "TMarker.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TGaxis.h"

#include <algorithm>
#include <vector>
#include <string>

using namespace ana;

void DrawHieTag(int hie, bool th13=true, bool vert=false);

void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
                const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
                const Int_t  kDarkColor, const TString surftype);

void DrawReactorConstraint(double mean, double err);

// kCMYK analog
void SetPaletteDeltaCyclicNew()
  {
    Double_t white[9]  = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
    Double_t Red[9]   = {  61./255.,  99./255., 136./255., 181./255., 213./255., 225./255., 198./255., 99./255., 61./255.};
    Double_t Green[9] = { 149./255., 140./255.,  96./255.,  83./255., 132./255., 178./255., 190./255., 140./255., 149./255};
    Double_t Blue[9]  = { 214./255., 203./255., 168./255., 135./255., 110./255., 100./255., 111./255., 203./255., 214./255.};
    Double_t Length[9] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};

    TColor::CreateGradientColorTable(9, Length, Red, Green, Blue, 255, 1.0);
    gStyle->SetNumberContours(255);
  }


void plot_joint_fit_2020_contours(bool corrSysts = false,
                                  bool runOnGrid = false,
                                  TString options="joint_fake2019_both_onlyIH",
                                  bool dmsqSurf = false,
                                  bool th13Surf = false,
                                  bool fccorr=false,
                                  bool zoomIn   = true,
                                  bool fillContour = false,
                                  bool smooth = true,
                                  bool fc_overlay = false)
{
  bool debug = 1;
  bool nueOnly = options.Contains("nueOnly");
  bool numuOnly = options.Contains("numuOnly");
  bool joint = options.Contains("joint");
  assert (nueOnly || numuOnly || joint);

  bool FHCOnly = options.Contains("FHCOnly");
  bool RHCOnly = options.Contains("RHCOnly");
  bool both = options.Contains("both");
  assert (FHCOnly || RHCOnly || both);

  bool fake2019 = options.Contains("fake2019");
  bool realData = options.Contains("realData");

  bool onlyNH = options.Contains("onlyNH");
  bool onlyIH = options.Contains("onlyIH");

  auto suffix = options;
  if(dmsqSurf) suffix+= "_dmsq";
  if(th13Surf) suffix+= "_th13";
  if(corrSysts) suffix+="_systs";
  else suffix+="_stats";

  TString outdir = "/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/";
  if (runOnGrid) outdir = "./";
  TString outFCfilename (outdir +"contours_FCInput_2020_" + suffix);
  TString outfilename   (outdir +  "hist_contours_2020_"    + suffix);


  //////////////////////////////////////////////////////////////////////
  // Make plots
  //////////////////////////////////////////////////////////////////////

  /*TString infilename = "/nova/ana/nu_e_ana/Ana2019/Results/RHC_and_FHC/contours/";
  if(th13Surf) infilename += "th13_delta/";
  else if (dmsqSurf) infilename += "th23_dmsq/";
  else infilename += "delta_th23/";
  infilename += corrSysts?"syst/":"stat/";*/
  TString infilename = outdir;
  infilename += "hist_contours_2020_";
  infilename += options;
  if (dmsqSurf) infilename += "_dmsq";
  infilename += corrSysts?"_systs.root":"_stats.root";

  TFile * infile = new TFile (infilename,"read");
  TFile * official =0;
  TString officialpreffix = "./";
  bool makedatarelease = true;

  if(!makedatarelease && corrSysts && options.Contains("joint") && fccorr && 0) {
     official = new TFile("NOvA_nue_official_results.root","recreate");
     TH2* axes = new TH2F("axes",";#delta_{CP}/#pi;sin^{2}#theta_{23}",100,0,2,100,0.225,0.725);
     //CenterTitles(axes);
     axes->Write();
  }
  else{
    if (makedatarelease){
       official = MakeDataReleaseFileContour(officialpreffix, dmsqSurf);
    }
  }

  auto mins =* (TVectorD*)infile->Get("overall_minchi");
  double minchi23 = mins[0];

  std::vector <TString> surfNames;// = {"delta_th23_"};<
  if(!th13Surf && ! dmsqSurf)surfNames = {"delta_th23_"};
  if (dmsqSurf) surfNames.push_back("th23_dm32_");
  if (th13Surf) surfNames.push_back("th13_delta_");

  for(TString surfName:surfNames){
      for (int hie:{-1,+1}){

        if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;

        //Setup canvas
        if(surfName.Contains("delta_th23")) {
          if (zoomIn) DrawContourCanvas(surfName, 0., 2., (RHCOnly?0:0.275), (RHCOnly?1.0:0.725));
          else DrawContourCanvas(surfName, 0., 2., 0., 1.);
        }
        if(surfName.Contains("th23_dm32")) {
          if (hie>0) DrawContourCanvas(surfName, ((RHCOnly)?0.2:0.3), ((RHCOnly)?0.8:0.7), (RHCOnly?2.0:2.075), (RHCOnly?3.3:2.725));
          else DrawContourCanvas(surfName, ((RHCOnly)?0.2:0.3), ((RHCOnly)?0.8:0.7), (RHCOnly?-3.0:-2.8), (RHCOnly?-2.0:-2.2));
        }
        if(surfName.Contains("th13_delta")) {
          if (zoomIn) DrawContourCanvas(surfName, 0., 0.38, 0, 2);
          else DrawContourCanvas(surfName, 0., 1.0, 0, 2);
        }
        //Load the surface; fill and FC corrections if available
        auto surf = *FrequentistSurface::LoadFrom(infile, "surf_"+surfName+(hie>0?"NH":"IH"));
        TH2 * hsurf = (TH2*)infile->Get(surfName+(hie>0?"NH":"IH"));
        TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
        TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
        TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
        TH2 * fc_surf_1Sigma = new TH2F();
        TH2 * fc_surf_2Sigma = new TH2F();
        TH2 * fc_surf_3Sigma = new TH2F();
        TH2 * fc_surf_90CL = new TH2F();
        if(fccorr){
          // if(!joint || !corrSysts || decomp != "combo" || dmsqSurf || th13Surf)
          //      { std::cerr << "\nWe don't have FC\n\n"; exit(1);}

          TString fcFolder="/nova/ana/nu_e_ana/Ana2020/FC/";
          //TString fcFolder = outdir;
          TString fcFileName;
          if(dmsqSurf) fcFileName = "contours/ssth23dmsq/ssth23dmsq_";
          else if(!dmsqSurf && !th13Surf)
            fcFileName = "contours/deltassth23/deltassth23_";
          fcFileName += hie > 0? "nh":"ih";
          fcFileName += ".root";

          FCSurface *fcXH = 0;
          auto fccol = FCCollection::LoadFromFile((fcFolder + fcFileName).Data()).release();
          int NX = hsurf->GetXaxis()->GetNbins();
          int NY = hsurf->GetYaxis()->GetNbins();
          double minX = hsurf->GetXaxis()->GetBinLowEdge(1);
          double maxX = hsurf->GetXaxis()->GetBinUpEdge(NX);
          double minY = hsurf->GetYaxis()->GetBinLowEdge(1);
          double maxY = hsurf->GetYaxis()->GetBinUpEdge(NY);

          std::cout << "Making FCSurface: \n";
          std::cout << "-------------------------------\n";
          std::cout << "NX: " << NX;
          std::cout << "\tNY: " << NY;
          std::cout << "\tminX: " << minX;
          std::cout << "\tmaxX: " << maxX;
          std::cout << "\tminY: " << minY;
          std::cout << "\tmaxY: " << maxY;
          std::cout << "\n------------------------------\n";

          fcXH = new FCSurface(NX, minX, maxX, NY, minY, maxY);
          std::string plot_name = "deltassth23";
          if(dmsqSurf) plot_name = "ssth23dmsq32";
          if(th13Surf) plot_name = "deltath13";

          fcXH->Add(*fccol, plot_name);
          fc_surf_1Sigma = fcXH->SurfaceForSignificance(0.6827);
          fc_surf_2Sigma = fcXH->SurfaceForSignificance(0.9545);
          fc_surf_3Sigma = fcXH->SurfaceForSignificance(0.9973);
          fc_surf_90CL   = fcXH->SurfaceForSignificance(0.90);
          if(smooth){
            fc_surf_1Sigma->Smooth();
            fc_surf_2Sigma->Smooth();
            fc_surf_90CL->Smooth();
            if(dmsqSurf || th13Surf)
              fc_surf_3Sigma->Smooth();
            else
              CylindricalSmoothing(fc_surf_3Sigma);
          }
        }
        Color_t k2SigmaConfidenceColor = ((hie>0)? cnh:cih);
        Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.99);
        Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.5);
        Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.2);
        Int_t kFCFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.65: 0.75);
        Int_t kFCFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.4: 0.38);
        Int_t kFCFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.18: 0.18);
        Int_t kDarkColor =  kGray+3;

        if(fillContour){
          if(!fccorr){
            auto g1 = surf.GetGraphs(surf_1Sigma,minchi23);
            auto g2 = surf.GetGraphs(surf_2Sigma,minchi23);
            auto g3 = surf.GetGraphs(surf_3Sigma,minchi23);
            FillGraphs(g1, g2, g3,
                       kFillColor1, kFillColor2, kFillColor3,
                       kDarkColor, suffix + surfName + (hie>0?"NH":"IH"));
          }
          else {
            auto fc_g1 = surf.GetGraphs(fc_surf_1Sigma,minchi23);
            auto fc_g2 = surf.GetGraphs(fc_surf_2Sigma,minchi23);
            auto fc_g3 = surf.GetGraphs(fc_surf_3Sigma,minchi23);
            FillGraphs(fc_g1, fc_g2, fc_g3,
                       kFillColor1, kFillColor2, kFillColor3,
                       kDarkColor, suffix + surfName + (hie>0?"NH":"IH")+ "fccorr");
          }
        }
        else{
          cout<<"minchi from file is "<<minchi23<<endl;
          if(!fccorr || fc_overlay){

            if(fc_overlay) {
              auto g1 = surf.GetGraphs(surf_1Sigma,minchi23);
              auto g2 = surf.GetGraphs(surf_2Sigma,minchi23);
              auto g3 = surf.GetGraphs(surf_3Sigma,minchi23);
              FillGraphs(g1, g2, g3,
                         kFillColor1, kFillColor2, kFillColor3,
                         kDarkColor, suffix + surfName + (hie>0?"NH":"IH"));
            }
            surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
            surf.DrawContour(surf_2Sigma, k2SigmaConfidenceStyle, kFillColor2,minchi23);
            surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
          }
          if(fccorr) {
            surf.DrawContour(fc_surf_1Sigma, kSolid, kFCFillColor1, minchi23);
            surf.DrawContour(fc_surf_2Sigma, k2SigmaConfidenceStyle, kFCFillColor2, minchi23);
            surf.DrawContour(fc_surf_3Sigma, kSolid, kFCFillColor3, minchi23);
          }
        }

        gPad->RedrawAxis();

        bool bestfit = 1;
        double  delta, th23, dmsq32, th13;
        TMarker bf;
        if(bestfit && both && corrSysts && hie>0){
           TFile * bffile = new TFile("/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/bestfits_joint_realData_both_systs.root","READ"); 
           auto calc_tmp = LoadFrom<osc::IOscCalc>(bffile, "NH_UO_calc").release();
           auto calc = (osc::IOscCalcAdjustable*)calc_tmp->Copy(); 
           if(!th13Surf){
             delta = kFitDeltaInPiUnits.GetValue(calc);
             th23 = kFitSinSqTheta23.GetValue(calc);
             dmsq32 = kFitDmSq32Scaled.GetValue(calc);
             th13 = kFitSinSq2Theta13.GetValue(calc);
           }
           else{delta = 0.678516; th13 = 0.0957894; cout<<"NB! Old 2019 values"<<endl;} /// NB! Old 2019 values
           bf.SetMarkerColor(kDarkColor);
          // bf.SetMarkerStyle(kFullStar);
           bf.SetMarkerStyle(43);
           bf.SetMarkerSize(2);
           if(!dmsqSurf && ! th13Surf){
              bf.SetX(delta);
              bf.SetY(th23);
              bf.DrawClone();
           }
           if(dmsqSurf){
              bf.SetX(th23);
              bf.SetY(dmsq32);
              bf.DrawClone();
           }
           if(th13Surf){
             bf.SetX(th13);
             bf.SetY(delta);
             bf.DrawClone();
           }
        }
        auto leg = ContourLegend(hie, fillContour, fccorr,
                                 kFillColor1,kFillColor2,kFillColor3,kDarkColor,
                                 bestfit && both);

        if(th13Surf) {
          DrawReactorConstraint(0.082,0.004);
          leg->SetFillStyle(1001);
          leg->SetNColumns(1);
          leg->SetX1(0.7); leg->SetX2(0.8);
          leg->SetY1(0.4); leg->SetY2(0.85);
          if(!fccorr) leg->SetHeader("No FC");

          TGraph * dummy = new TGraph;
          dummy->SetLineWidth(2);
          dummy->SetLineColor(kGray+2);
          dummy->SetFillColor(kGray);
          leg->AddEntry(dummy,"#splitline{Reactor}{68%C.L.}","f");
        }
        leg->Draw();
        Preliminary();
        DrawHieTag(hie,false);


        for(TString ext: {".png", ".pdf",".root"}){
          gPad->Print("./contour_" +  suffix + "_" +
                      surfName +
                      (hie>0?"NH":"IH") + (fccorr?"_fccorr":"") +
                      ((fccorr && smooth)?"_smooth":"") +
                      (zoomIn?"_zoom":"") +
                      (fc_overlay?"_overlay":"") +
                      ext);
        }//end ext
     
        if(official) {
          auto fc_g1 = surf.GetGraphs(fc_surf_1Sigma,minchi23);
          auto fc_g2 = surf.GetGraphs(fc_surf_2Sigma,minchi23);
          auto fc_g3 = surf.GetGraphs(fc_surf_3Sigma,minchi23);
          auto fc_g9 = surf.GetGraphs(fc_surf_90CL,minchi23);
          GraphsToDataRelease(official, fc_g1,"1sigma", hie>0?"NH":"IH");
          GraphsToDataRelease(official, fc_g2,"2sigma", hie>0?"NH":"IH");
          GraphsToDataRelease(official, fc_g3,"3sigma", hie>0?"NH":"IH");
          GraphsToDataRelease(official, fc_g9,"90CL", hie>0?"NH":"IH");
          if(hie>0)  //gah
            BestFitToDataRelease(official, bf, "best_fit",hie>0?"NH":"IH");
        }

        debug = false; 
        if(debug && !fccorr){
          auto hists = surf.GetProfiledHists();
          int nSysts = corrSysts ? 70 : 0;
          for (int i = 0; i < nSysts; ++i){
            hists.push_back((TH2*) infile->Get((TString) "surf_" + surfName +
                                               (hie > 0 ? "NH" : "IH") +
                                               "/profHists/hist" + std::to_string(i)));
          }
          hists.push_back(surf.ToTH2(minchi23));

          for(uint i = 0; i < hists.size(); ++i){
            if(!hists[i]) continue;
            if(TString(hists[i]->GetTitle()).Contains("delta")) {
              SetPaletteDeltaCyclicNew();
              hists[i]->GetZaxis()->SetRangeUser(0,2);
            }
            else gStyle->SetPalette(87);
            if( i > 3)  {
              double minz = hists[i]->GetBinContent(hists[i]->GetMinimumBin());
              double maxz = hists[i]->GetBinContent(hists[i]->GetMaximumBin());
              double lim = std::max(fabs(minz), fabs(maxz));
              hists[i]->GetZaxis()->SetRangeUser(-lim,+lim);
            }
            if(i==(hists.size()-1)) {
              gStyle->SetPalette(56);
              hists[i]->GetZaxis()->SetRangeUser(0.,11.83);
            }

            hists[i]->Draw("colz");
            hists[i]->GetYaxis()->SetTitleOffset(0.8);
            gPad->SetRightMargin(0.14);
            surf.DrawContour(surf_1Sigma, kSolid, kDarkColor,minchi23);
            surf.DrawContour(surf_2Sigma, kSolid, kDarkColor,minchi23);
            surf.DrawContour(surf_3Sigma, kSolid, kDarkColor,minchi23);
            gPad->Print("./debug_contour_" +  suffix + "_" +
                        surfName +
                        (hie > 0 ? "NH":"IH") + (fccorr?"_fccorr":"") +
                        (zoomIn ? "_zoom" : "") +
                        (std::to_string(i)) + ".pdf");

          }//end profiles
        }//end debug
      }//end hie
    }//end surf
}//end all


void DrawHieTag(int hie,bool th13, bool vert){
  if(hie>0){
    TLatex * ltxNH = new TLatex(0.16,0.89,"#bf{NH}");
    ltxNH->SetTextColor(cnh);
    ltxNH->SetNDC();
    ltxNH->SetTextAlign(22);
    ltxNH->SetTextSize(kBlessedLabelFontSize);
    if(th13)    ltxNH->Draw();
    else if(!vert) ltxNH->DrawLatex(.85,.20,"#bf{NH}");
    else ltxNH->DrawLatex(.85,.85,"#bf{NH}");
  }
  else{
    TLatex * ltxIH = new TLatex (0.16,0.46,"#bf{IH}");
    ltxIH->SetTextColor(cih_dark);
    ltxIH->SetNDC();
    ltxIH->SetTextAlign(22);
    ltxIH->SetTextSize(kBlessedLabelFontSize);
    if(th13)    ltxIH->Draw();
    else if(!vert)  ltxIH->DrawLatex(.85,.20,"#bf{IH}");
    else ltxIH->DrawLatex(.85,.85,"#bf{IH}");
  }
}


void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
                const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
                const Int_t  kDarkColor, const TString surftype)
{
  bool isJoint = surftype.Contains("joint");
  bool isBoth = surftype.Contains("both");
  bool isRHC = surftype.Contains("RHCOnly");
  bool isNH = surftype.Contains("NH");
  bool isSysts =  surftype.Contains("systs");
  bool isFccorr =  surftype.Contains("fccorr");

  if (surftype.Contains("delta_th23")){
    cout<<"delta"<<endl;
    if(isJoint && isNH ){
      cout<<"joint nh"<<endl;
      if(isRHC){
         JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
         JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
         JoinGraphs(g2[1], g2[2], kWhite)->Draw("f");
         JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f same");
         JoinGraphs(g2[1], g2[2], kFillColor2)->Draw("f same");
      }
      else{
        JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
        JoinGraphs(g3[0], g3[1], kWhite)->Draw("l");
        JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
        JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
        JoinGraphs(g2[0], g2[1], kWhite)->Draw("l");
        //JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("c f");
        //JoinGraphs(g1[0], g1[1], kWhite)->Draw("c");
        if(isBoth){
           for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("c f");}
           for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("c f");}
           for (auto & g:g1) { g->Draw("c");}
        }
        else{
          JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("c f");
          JoinGraphs(g1[0], g1[1], kWhite)->Draw("c");
        }
      }
    }

    if(isJoint && !isNH){
      if(isRHC)  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
      else{
          //JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f c");
          g3[1]->SetFillColor(kFillColor3); g3[1]->Draw("f c");
          g3[0]->SetFillColor(kFillColor3); g3[0]->Draw("f c");
          for (auto &g:g3) {
              g->SetLineWidth(3);
              g->Draw("c");
          }
       }
        for (auto &g:g2) {
          g->SetFillColor(kWhite); g->DrawClone("f c");
          g->SetFillColor(kFillColor2); g->Draw("f c");
          g->SetLineColor(kDarkColor);
          g->SetLineWidth(3);
          g->Draw("c");
      }
       for (auto &g:g1) {
          g->SetFillColor(kWhite); g->DrawClone("f c");
          g->SetFillColor(kFillColor1); g->Draw("f c");
          g->SetLineColor(kDarkColor);
          g->SetLineWidth(3);
          g->Draw("c");
      }
       for(auto & gs:{g3, g2, g1}){
         for(auto & g:gs) {
           g->SetLineColor(kDarkColor);
           g->SetLineWidth(3);
           g->Draw("c");
         }
       }
    }
 }

  if (surftype.Contains("th23_dm32")){
    for (auto &g:g3) {
      g->SetFillColor(kFillColor3); g->Draw("cf");
    }
    for (auto &g:g2) {
      g->SetFillColor(kWhite); g->DrawClone("cf");
      g->SetFillColor(kFillColor2); g->Draw("cf");
    }
    for (auto &g:g1) {
      g->SetFillColor(kWhite); g->DrawClone("cf");
      g->SetFillColor(kFillColor1); g->Draw("cf");
    }
    for (auto  & gs:{g3,g2,g1}){
       for (auto & g:gs) {
            g->SetLineColor(kDarkColor);
            g->SetLineWidth(3);
            g->Draw("c");
       }
     }

  }

  if (surftype.Contains("th13_delta")) {
    if(isBoth) {
       JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
       if(isNH) {
         JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
         JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
       }
       else{
         for (auto & g:g2) { g->SetFillColor(kWhite); g->DrawClone("f");}
         for (auto & g:g2) { g->SetFillColor(kFillColor2); g->Draw("f");}
       }
       for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("f");}
       for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("f");}
    }
    else{
       //for (auto &g:g3) {
       //    g->SetFillColor(kFillColor3); g->Draw("f");}
        //JoinGraphs(g3[0], g3[1], kWhite)->Draw("f");
       JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f c");
       JoinGraphs(g2[0], g2[1], kWhite)->Draw("f c");
       JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f c");
       /*for (auto &g:g2) {
           g->SetFillColor(kWhite); g->DrawClone("f");
           g->SetFillColor(kFillColor2); g->Draw("f");
       }*/
       JoinGraphs(g1[0], g1[1], kWhite)->Draw("f c");
       JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("f c");
   }

     for (auto  & gs:{g3,g2,g1}){
         for (auto & g:gs) {
              g->SetLineColor(kDarkColor);
              g->SetLineWidth(3);
              g->Draw("l");
         }
     }
  }

}

void DrawReactorConstraint(double mean, double err)
{
  TBox* box = new TBox(mean-err, 0, mean+err, 2);
  box->SetFillColorAlpha(kGray,0.7);

  TLine* linLo = new TLine(mean-err, 0, mean-err, 2);
  linLo->SetLineColor(kGray+2);
  linLo->SetLineWidth(2);

  TLine* linHi = new TLine(mean+err, 0, mean+err, 2);
  linHi->SetLineColor(kGray+2);
  linHi->SetLineWidth(2);

  box->Draw();
  linLo->Draw();
  linHi->Draw();
  gPad->RedrawAxis();
}

// ssth23dmsq32 NH
// cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyNH true false true true false true false
// ssth23dmsq32 IH
// cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyIH true false true true false true false

// deltassth23 NH
// cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyNH false false true true false true false
// deltassth23 IH
// cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyIH false false true true false true false