joint_fit_2017_contours.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Analysis/Plots.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Fit/FrequentistSurface.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 "3FlavorAna/Prediction/PredictionSystNue2017.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 "3FlavorAna/Ana2017/nue/joint_fit_2017_tools.h"

#include "TCanvas.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>

using namespace ana;

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

TLegend * ContourLegend(int hie, bool fillContour, bool fccorr,
                        Int_t kFillColor1, Int_t kFillColor2,
                        Int_t kFillColor3, Int_t kDarkColor,
                        bool bestFit);
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  DrawContourCanvas (TString surfName, double minx = 0,
                         double maxx = 2, double miny = 0, double maxy = 1){
  auto c1 = new TCanvas(ana::UniqueName().c_str());
  c1->SetFillStyle(0);
  c1->SetLeftMargin(0.14);
  c1->SetBottomMargin(0.15);
  TH2* axes = new TH2F();
  TString title;
  if(surfName.Contains("delta_th23")) title=";#delta_{CP};sin^{2}#theta_{23}";
  if(surfName.Contains("th13_delta")) title=";sin^{2}2#theta_{13};#delta_{CP}/#pi";
  if(surfName.Contains("th23_dm32"))  title=";sin^{2}#theta_{23};#Deltam^{2}_{32} (10^{-3}eV^{2})";
  axes = new TH2F("",title,100,minx,maxx,100,miny,maxy);

  CenterTitles(axes);
  axes->Draw();
  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
  axes->GetXaxis()->SetTitleOffset(0.8);
  axes->GetYaxis()->SetTitleOffset(0.8);
  TGaxis::SetMaxDigits(3);
  if(surfName.Contains("th23_dm32")) {
    axes->GetYaxis()->SetDecimals() ;
    axes->GetYaxis()->SetTitleOffset(0.85);
    axes->GetYaxis()->SetLabelOffset(0.002);
    axes->GetYaxis()->SetLabelSize(0.058);
    axes->GetYaxis()->SetTitleSize(0.078);
  }
  if(surfName.Contains("delta_th23"))  XAxisDeltaCPLabels(axes);
  c1->RedrawAxis();
}

void joint_fit_2017_contours(std::string decomp = "prop",
                             bool createFile=false,
                             bool corrSysts = false,
                             TString options="joint_fakeNHLO",
                             bool dmsqSurf = false,
                             bool th13Surf = false,
                             bool fccorr=false,
                             bool zoomIn   = true,
                             bool fillContour = false,
                             bool smooth = true)
{
  bool debug = 1;
  bool nueOnly = options.Contains("nueOnly");
  bool numuOnly = options.Contains("numuOnly");
  bool joint = options.Contains("joint");
  assert (nueOnly || numuOnly || joint);

  bool fakeNHLO = options.Contains("fakeNHLO");
  bool fakeNHUO = options.Contains("fakeNHUO");
  bool fake2017 = options.Contains("fake2017");
  bool realData = options.Contains("realData");

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

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

  TString outdir="/nova/ana/nu_e_ana/Ana2017/Results/contours/";
  if(th13Surf) outdir += "th13_delta/";
  else if (dmsqSurf) outdir += "th23_dmsq/";
  else outdir += "delta_th23/";

  //TString outdir = "contours/";

  TString outfilename   (outdir + "hist_contours_2017_"    + suffix);
  TString outFCfilename (outdir + "contours_FCInput_2017_" + decomp + ".root");

  if(createFile){

    //////////////////////////////////////////////////
    // Load Nue and Numu experiments
    //////////////////////////////////////////////////

    std::vector <const IPrediction * > preds;
    std::vector <std::pair <TH1D*, double > > cosmics;
    std::vector <Spectrum * > data;
    std::vector <const IExperiment * > expts;

    auto calc_fake = DefaultOscCalc();
    if(fakeNHLO) SetFakeCalc(calc_fake, 0.404, 2.7e-3, 1.48);
    else if(fakeNHUO) SetFakeCalc(calc_fake, 0.623, 2.7e-3, 0.74);
    else if(fake2017) SetFakeCalc(calc_fake, 0.56, 2.43e-3, 1.21);
    else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}

    if(!numuOnly) {
      preds.push_back(GetNuePrediction2017(decomp, DefaultOscCalc(), corrSysts));
      cosmics.push_back(GetNueCosmics2017());
    }
    int nnumu = 4;
    if(!nueOnly) {
      auto numu_preds = GetNumuPredictions2017(nnumu, corrSysts);
      preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
      auto numu_cosmics = GetNumuCosmics2017(nnumu);
      cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
    }
    if(realData){
      if(!numuOnly){
        data.push_back(GetNueData2017());
      }
      if(!nueOnly){
        auto numu_data = GetNumuData2017(nnumu);
        data.insert(data.end(),numu_data.begin(), numu_data.end());
      }
    }
    for(int i = 0; i < int(preds.size()); ++i){
      auto thisdata = GetFakeData (preds[i],calc_fake, kAna2017POT,
                                   cosmics[i].first);
      if(realData) thisdata = data[i];
      expts.push_back(new SingleSampleExperiment(preds[i],*thisdata,
                                                 cosmics[i].first,cosmics[i].second));

      if(debug){
        new TCanvas();
        DataMCComparison (*thisdata, preds[i], DefaultOscCalc());
        cosmics[i].first->SetMarkerStyle(34);
        cosmics[i].first->SetMarkerColor(kCosmicBackgroundColor);
        cosmics[i].first->Draw("hist p same");
        gPad->Print(outdir + "debug_predictions_" + suffix + "_" + std::to_string(i) + ".pdf");
      }//end debug
    }

    //////////////////////////////////////////////////
    // Add constraints, make experiments
    //////////////////////////////////////////////////
    std::cout << "\nCreating multiexperiment\n" << std::endl;
    if(!th13Surf){
      std::cout << "Adding WorldReactorConstraint2017\n";
      expts.push_back(WorldReactorConstraint2017());
    }
    if(nueOnly) {
      std::cout << "Adding Dmsq32ConstraintPDG2017\n";
      expts.push_back(&kDmsq32ConstraintPDG2017);
    }
    std::cout << "Creating Multiexperiment with a total of "
              << expts.size() << " experiments\n\n" << std::endl;
    auto exptThis = new MultiExperiment(expts);

    ////////////////////////////////////////////////////////////
    // Systematics
    ////////////////////////////////////////////////////////////
    std::cout << "Systematics for the fit:\n";
    auto systs = GetJointFitSystematicList(corrSysts, !numuOnly, !nueOnly, true);

//    std::cout << "HACK MECQ0\n";
//    systs = {kMECq0ShapeSyst2017};

    std::cout << "\n\nSystematic correlations...\n";
    if(!nueOnly && ! numuOnly && corrSysts){
      exptThis->SetSystCorrelations(0, GetCorrelations(true));
      auto notfornumu = GetCorrelations(false);
      for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
    }

    ////////////////////////////////////////////////////////////
    // Fit
    ////////////////////////////////////////////////////////////
    std::cout << "Starting the fit" << std::endl;

    osc::IOscCalcAdjustable* calc = DefaultOscCalc();

    std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits,&kFitSinSqTheta23,
          &kFitDmSq32Scaled,&kFitSinSq2Theta13};
    if (th13Surf) fitvars = {&kFitDeltaInPiUnits,&kFitSinSqTheta23,&kFitDmSq32Scaled};

    //Find the best fit
    MinuitFitter fit23(exptThis,fitvars,systs);

    SystShifts auxShifts = SystShifts::Nominal();
    // In case some systs need different seeds
    std::vector <SystShifts> seedShifts = {};
    if(corrSysts) {
      for (double systshift:{+0.5,-0.5}){
        SystShifts tempShifts (&kMECq0ShapeSyst2017,systshift);
        seedShifts.push_back(tempShifts);
      }
    }

    double minchi23 = 1E20;
    for(double thseed:{0.45,0.55}){
      for(int hie:{-1,1}){
        std::cout << "\n\nFinding best fit "
                  << (thseed<0.5?"LO ":"UO ")
                  << (hie>0? "NH " : "IH ")
                  << "\n\n";
        auto thisminchi = fit23.Fit(calc,auxShifts ,{
            {&kFitDmSq32,{hie*fabs(calc->GetDmsq32())}},
            {&kFitSinSqTheta23,{thseed}},
            {&kFitDeltaInPiUnits,{0,0.5,1,1.5}}
          },  seedShifts
          );
        minchi23= min(minchi23, thisminchi);
        ResetOscCalcToDefault(calc);
        auxShifts.ResetToNominal();
      }//end hie
    }//end th23
    std::cout << "\nFound overall minchi " << minchi23 << std::endl;

    //We'll save to this file, make pretty plots later
    TFile * outfile = new TFile(outfilename+".root","recreate");
    //The FC also needs to know about the fitted syst values, save those in
    //a separate file
    TFile * outFCfile = new TFile(outFCfilename+".root","recreate");

    outfile->cd();
    TVectorD v(1);
    v[0] = minchi23;
    v.Write("overall_minchi");

    //HACK HACK
    int steps = 30;//10;// 5;
    if(dmsqSurf){
      std::cerr << "\n WARNING Using 20 bins for dmsq surface\n\n";
      steps = 20;
    }
    //Now create all the surfaces
    for(int hie: {-1, +1}){

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

      std::cout << "Starting surface " << (hie>0? "NH ": "IH") << "\n\n";
      std::string hieStr = hie>0 ? "NH" : "IH";
        
      if(!th13Surf && ! dmsqSurf){
        ResetOscCalcToDefault(calc);
        calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));

        FrequentistSurface surf23(exptThis, calc,
                       &kFitDeltaInPiUnits,steps, 0, 2,
                       &kFitSinSqTheta23,  steps, (nueOnly?0:0.3), (nueOnly?1:0.7),
                       {&kFitDmSq32Scaled,&kFitSinSq2Theta13},systs,
                       {},seedShifts);

        auto surf1=surf23.ToTH2(minchi23);
        outfile->cd();
        surf1->Write((TString)"delta_th23_"+hieStr);
        surf23.SaveTo(outfile, (TString)"surf_delta_th23_"+hieStr);

        outFCfile->cd();
        std::vector<TH2*> profhists = surf23.GetProfiledHists();
        //Osc parameters are special case
        profhists[0]->Write((hieStr+"_DmSq32").c_str());
        profhists[1]->Write((hieStr+"_SinSq2Theta13").c_str());
        //Also need to save a hist for each syst (new for 2017)
        for (int i = 0; i < int(systs.size()); i++){
          profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
        }
      }
      if(!th13Surf && dmsqSurf){
        calc = DefaultOscCalc();
        calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));

        FrequentistSurface surf23m(exptThis, calc,
                        &kFitSinSqTheta23, steps,0.3,0.7,
                        &kFitDmSq32Scaled, steps,
                        (hie>0?2.1:-2.8),(hie>0?2.8:-2.1),
                        {&kFitSinSq2Theta13,&kFitDeltaInPiUnits},systs,
                        {{&kFitDeltaInPiUnits,{0,0.5,1,1.5}}},seedShifts);
        auto surf6=surf23m.ToTH2(minchi23);
        outfile->cd();
        surf6->Write((TString)"th23_dm32_"+(hie>0?"NH":"IH"));
        surf23m.SaveTo(outfile, (TString)"surf_th23_dm32_"+(hie>0?"NH":"IH"));

        outFCfile->cd();
        std::vector<TH2*> profhists = surf23m.GetProfiledHists();
        //Osc parameters are special case
        profhists[0]->Write((hieStr+"_SinSq2Theta13").c_str());
        profhists[1]->Write((hieStr+"_DeltaCPpi").c_str());
        //Also need to save a hist for each syst (new for 2017)
        for (int i = 0; i < int(systs.size()); i++){
          profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
        }

      }
      if(th13Surf){
        calc = DefaultOscCalc();
        calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));

        FrequentistSurface surf13(exptThis, calc,
                       &kFitSinSq2Theta13,  steps, 0, 0.35,
                       &kFitDeltaInPiUnits,steps, 0, 2,
                       {&kFitSinSqTheta23, &kFitDmSq32Scaled}, systs,
                       {{&kFitSinSqTheta23,{0.45,0.55}}}, seedShifts);
        auto surf4 = surf13.ToTH2(minchi23);
        outfile->cd();
        surf4->Write((TString)"th13_delta_"+(hie>0?"NH":"IH"));
        surf13.SaveTo(outfile, (TString)"surf_th13_delta_"+(hie>0?"NH":"IH"));

        outFCfile->cd();
        std::vector<TH2*> profhists = surf13.GetProfiledHists();
        //Osc parameters are special case
        profhists[0]->Write((hieStr+"_SinSqTheta23").c_str());
        profhists[1]->Write((hieStr+"_DmSq32").c_str());
        //Also need to save a hist for each syst (new for 2017)
        for (int i = 0; i < int(systs.size()); i++){
          profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
        }

      }
    } //end hie

    delete outfile;
    std::cout << "\nSurfaces saved to " << outfilename << std::endl;
    return;
  }//end createfile

//////////////////////////////////////////////////////////////////////
// Make plots
//////////////////////////////////////////////////////////////////////
  else{ //plot only
    TFile * infile = new TFile (outfilename+".root","read");
    TFile * official =0;
    if(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();
    }

    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., 0.225, 0.725);
          else DrawContourCanvas(surfName, 0., 2., 0., 1.);
        }
        if(surfName.Contains("th23_dm32")) {
          if (hie>0) DrawContourCanvas(surfName, 0.3, 0.7, 2.05, 2.75);
          else DrawContourCanvas(surfName, 0.3, 0.7, -2.95, -2.25);
        }
        if(surfName.Contains("th13_delta")) {
          if (zoomIn) DrawContourCanvas(surfName, 0., 0.38, 0, 2);
          else DrawContourCanvas(surfName, 0., 0.5, 0, 2);
        }
        //Load the surface; fill and FC corrections if available
        auto surf = *ISurface::LoadFrom(infile, "surf_" + surfName + (hie > 0 ? "NH" : "IH"));

        TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
        TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
        TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);

        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/Ana2017/FC/";
          TString fcFileName ((TString)"FCCol_"+(hie>0?"NH":"IH")+"_surf_ssth23delta.root");
          FCSurface *fcXH = 0;
          auto fccol = FCCollection::LoadFromFile((fcFolder + fcFileName).Data()).release();
          fcXH = new FCSurface(30,0,2,30,0.3,0.7);
          fcXH->Add(*fccol);
          surf_1Sigma = fcXH->SurfaceForSignificance(0.6827);
          surf_2Sigma = fcXH->SurfaceForSignificance(0.9545);
          surf_3Sigma = fcXH->SurfaceForSignificance(0.9973);
          if(smooth){
            surf_1Sigma->Smooth();
            surf_2Sigma->Smooth();
            surf_3Sigma->Smooth();
          }
        }

        Color_t k2SigmaConfidenceColor = ((hie>0)? k2SigmaConfidenceColorNH:k2SigmaConfidenceColorIH);
        Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.70);
        Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.35);
        Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.20);
        Int_t kDarkColor =  kGray+3;

        if(fillContour){
          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") + (fccorr?"_fccorr":""));
        }
        else{
          surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
          surf.DrawContour(surf_2Sigma, k2SigmaConfidenceStyle, kFillColor2,minchi23);
          surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
        }

        gPad->RedrawAxis();

        bool bestfit = 0; //todo
        auto leg = ContourLegend(hie, fillContour, fccorr,
                                 kFillColor1,kFillColor2,kFillColor3,kDarkColor,
                                 bestfit);
        if(th13Surf) {
          leg->SetFillStyle(1001);
          leg->SetNColumns(1);
          leg->SetX1(0.76); leg->SetX2(0.89);
          leg->SetY1(0.52); leg->SetY2(0.85);
          leg->SetHeader("No FC");
        }
        leg->Draw();
        DrawHieTag(hie,false);

        for(TString ext: {".pdf",".root"}){
          gPad->Print(outdir + "contour_" +  suffix + "_" +
                      surfName +
                      (hie>0?"NH":"IH") + (fccorr?"_fccorr":"") +
                      ((fccorr && smooth)?"_smooth":"") +
                      (zoomIn?"_zoom":"") +
                      ext);
        }//end ext

        if(debug && !fccorr){
          auto hists = surf.GetProfiledHists();
          int nSysts = corrSysts ? 30 : 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")) {
              SetPaletteBlueRedCyclic();
              hists[i]->GetZaxis()->SetRangeUser(0,2);
            }
            else SetPaletteBlueRedWhite();
            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)) {
              SetPaletteWhiteBlueDark();
              hists[i]->GetZaxis()->SetRangeUser(0.1,20);
              gPad->SetLogz();
            }

            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(outdir + "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 plot only
}//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(kPrimColorNH);
    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(kPrimColorIH);
    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}");
  }
}

TLegend * ContourLegend(int hie, bool fillContour, bool fccorr,
                        Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor,
                        bool bestFit){
  TLegend * leg = new TLegend(0.16,fccorr?0.19:0.17,0.65,fccorr?0.26:0.28);
  leg->SetTextSize(kBlessedLabelFontSize);
  leg->SetFillStyle(0);
  leg->SetMargin(1.3*leg->GetMargin());
  leg->SetNColumns(3);

  TGraph * dummy = new TGraph;
  dummy->SetLineWidth(3);

  if(fillContour){
    dummy->SetLineColor(kDarkColor);
    dummy->SetFillColor(kFillColor1);
    leg->AddEntry(dummy->Clone(),"1 #sigma","f");
    dummy->SetFillColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2 #sigma","f");
    dummy->SetFillColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3 #sigma","f");
  }
  else{
    leg->SetMargin(1.4*leg->GetMargin());
    dummy->SetLineColor(kFillColor1);
    dummy->SetLineStyle(kSolid);
    leg->AddEntry(dummy->Clone(),"1 #sigma","l");
    dummy->SetLineStyle(k2SigmaConfidenceStyle);
    dummy->SetLineColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2 #sigma","l");
    dummy->SetLineStyle(kSolid);
    dummy->SetLineColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3 #sigma","l");
  }
  if(bestFit){
    leg->SetNColumns(4);
//    leg->SetColumnSeparation(0.08);
    leg->SetX2(0.75);
    dummy->SetMarkerStyle(kFullCircle);
    dummy->SetMarkerColor(kDarkColor);
    dummy->SetLineColor(kDarkColor);
    dummy->SetLineStyle(7);
    if(hie>0) leg->AddEntry(dummy->Clone(),"Best Fit","p");
    if(hie<0) leg->AddEntry((TObject*)0, "#color[0]{Best Fit}", "");
    dummy->SetLineStyle(1);
  }
  if(!fccorr) leg->SetHeader("No Feldman-Cousins");
  return leg;
}

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 isNH = surftype.Contains("NH");
  bool isSysts =  surftype.Contains("systs");
  bool isFccorr =  surftype.Contains("fccorr");

  if (surftype.Contains("delta_th23")){
    if(isJoint && isNH ){
      JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
      JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
      JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
      if(isSysts && !isFccorr) JoinGraphs(g1[0], g1[1], kWhite)->Draw("f");
      if(isSysts && !isFccorr) JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("f");
      if(!isSysts || isFccorr) for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("f");}
      if(!isSysts || isFccorr) for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("f");}
    }

    if(isJoint && !isNH){
      for (auto &g:g3) {
        g->SetFillColor(kFillColor3); g->Draw("f");}
      for (auto &g:g2) {
        g->SetFillColor(kWhite); g->DrawClone("f");
        g->SetFillColor(kFillColor2); g->Draw("f");
      }
    }
  }

  if (surftype.Contains("th23_dm32")){
    for (auto &g:g3) {
      g->SetFillColor(kFillColor3); g->Draw("f");}
    for (auto &g:g2) {
      g->SetFillColor(kWhite); g->DrawClone("f");
      g->SetFillColor(kFillColor2); g->Draw("f");
    }
    for (auto &g:g1) {
      g->SetFillColor(kWhite); g->DrawClone("f");
      g->SetFillColor(kFillColor1); g->Draw("f");
    }
  }

  if (surftype.Contains("th13_delta")) {
    JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
    JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
    JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
    JoinGraphs(g1[0], g1[1], kWhite)->Draw("f");
    JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("f");
  }


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

}