Tutorial2019FitContours.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Fit.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/NuePlotStyle.h"
#include "CAFAna/Analysis/Surface.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/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/FitVars.h"
#include "Utilities/rootlogon.C"
#include "OscLib/IOscCalc.h"

#include "CAFAna/nue/Ana2019/joint_fit_2019_loader_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;

std::vector <Spectrum * > GetMockData(TString ana_options);

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);     
  }

/// Start the Demo script (it's a bit simplified nue / Ana2018 / FitandFC / joint_fit_2018_contours.C):

void Tutorial2019FitContours(bool createFile=false,               /// Do you want to fit and write it into the file (true) or make plots (false)?
                             bool corrSysts = true,               /// Stat only (false) or with systs (true)?
                             TString options="joint_realData_both_onlyIH", /// Write exactly what fit is it? joint/nueOnly/numuOnly + both/FHCOnly/RHCOnly + onlyIH/onlyNH
                             bool dmsqSurf = false,               /// What kind of plot do you want, is it dmsq vs th23 surface?
                             bool th13Surf = false,               /// What kind of plot do you want, is it th13 vs delta CP surface? 
                                                                  /// If dmsqSurf = false && th13Surf = false then it will be th23 vs delta CP surface
                             bool zoomIn   = true,                /// Since it's pretty universal code axes are very extanded for cases of nueOnly || numuOnly || RHCOnly || FHCOnly. 
                                                                  /// Thus sometimes you'll want to zoom axis (especially in case of FHC&&RHC)
                             bool fillContour = false)            /// Do you want to fill it with color or stay with border lines only?
{
  /// Some general info about further fit 

  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";

  ///How do you want to organize and sort your outputs? 

  /*TString outdir="";

  if(th13Surf) outdir += "th13_delta/";
  else if (dmsqSurf) outdir += "th23_dmsq/";
  else outdir += "delta_th23/";

  TString outfilename   (outdir + (corrSysts?"syst/":"stat/") + "hist_contours_2018_"    + suffix);*/

  /// For tute purposes save all plots in current dir
  TString outdir="";
  TString outfilename   (outdir + "hist_contours_2019_"    + suffix);

  /// Enter the fit and file creation part

  if(createFile){

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

    /// Create vectors for all predictions we need - MC predictions, cosmic predictions, data or fake data
    /// Then we'll need vectors of experiments for making a multiexperiment 
    /// And we'll need separate set of all items for making numu-only preseeds (since we follow 2018 ana procedure)

    struct predictions {
           const string name;
           const IPrediction * pred;
           std::pair <TH1D*, double> cos;
           double pot;
    };

    std::vector <predictions> preds;
    std::vector <Spectrum*> data;
    std::vector <const IExperiment * > expts;

    auto calc_fake = DefaultOscCalc();
    if(fake2019) SetFakeCalc(calc_fake, 0.565, 2.48e-3, 1.99); /// use 2019 BF for making fake data

    if(!numuOnly) {
      if(FHCOnly || both ) {
         preds.push_back({"nue_fhc", GetNuePrediction2019("combo", DefaultOscCalc(), corrSysts, "fhc", false), GetNueCosmics2019("fhc"), GetPOT()}); /// prediction downloading functions are defined in the ..._tools.h header, GetPOT () returns the POT exposure value for rhc or fhc mode
      }
      if(RHCOnly || both ) {
         preds.push_back({"nue_fhc", GetNuePrediction2019("prop", DefaultOscCalc(), corrSysts, "rhc", false), GetNueCosmics2019("rhc"), GetPOT(false)});
      }
    }

   int nnumu = 4;
   if(!nueOnly) {
      if(FHCOnly || both ) {
        auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "fhc");
        auto numu_cosmics = GetNumuCosmics2019(nnumu, "fhc");
        for (int i = 0; i< nnumu; i++ ){
            preds.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
        }
      }
      if(RHCOnly || both ) {
        auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "rhc");
        auto numu_cosmics = GetNumuCosmics2019(nnumu, "rhc");
        for (int i = 0; i< nnumu; i++ ){
            preds.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
        }
       }
    }    

    /// Now assign proper POT for each of these predictions following the same structure
    /// nue FHC, nue RHC, 4 numu FHC and 4 numu RHC quartiles
    /// end make fake data. Feed everything to the experiments
    
    cout<<"\n\n Predictions used for the fit:\n";

    if(realData){
      if(!numuOnly){
        if(FHCOnly || both ) {
           auto temp = GetMockData("nue_fhc");
           data.insert(data.end(),temp.begin(), temp.end());
        }
        if(RHCOnly || both ) {
           auto temp = GetMockData("nue_rhc");
           data.insert(data.end(),temp.begin(), temp.end());
        }
      }
      if(!nueOnly){
        if(FHCOnly || both ){
          auto temp = GetMockData("numu_fhc");
          data.insert(data.end(), temp.begin(), temp.end());}
        if(RHCOnly || both ){
          auto temp = GetMockData("numu_rhc");
          data.insert(data.end(), temp.begin(), temp.end());}
      }
    }

    for(int i = 0; i < int(preds.size()); ++i){
        double POT = preds[i].pot;
        auto thisdata = GetFakeData (preds[i].pred, calc_fake, POT, preds[i].cos.first);
        if(realData) thisdata = data[i];
        cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds[i].cos.first->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
      expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, preds[i].cos.first, preds[i].cos.second));
     
    }
     

//////////////////////////////////////////////////
/// Add constraints, make experiments
//////////////////////////////////////////////////

    /// Lets make a multiexperiment for simultaneous fit 
    /// For different contours/fits we will use different additional information 

    std::cout << "\nAll experiments are done \nCreating multiexperiment\n" << std::endl;
    if(!th13Surf){ /// We use reactor constraint as 11th "experiment" - th13 will be varied in PDG's error range with PDG central value
                   /// but we wouldn't use such constraint if we want to make our own measurement, right?
      std::cout << "Adding WorldReactorConstraint2017\n";
      expts.push_back(WorldReactorConstraint2017());
    }
    if(nueOnly) { /// If we don't use NOvA's result in the fit, then we will use PDG's constraint
      std::cout << "Adding Dmsq32ConstraintPDG2017\n";
      expts.push_back(&kDmsq32ConstraintPDG2017);
    }
    /// Let's make a multiexperiment
    std::cout << "Creating Multiexperiment with a total of "
              << expts.size() << " experiments\n\n" << std::endl;
    auto exptThis = new MultiExperiment(expts);

    std::cout<<"Clean verctors"<<std::endl;
    preds.clear();

////////////////////////////////////////////////////////////
/// Systematics
////////////////////////////////////////////////////////////

   /// Now, let's prepare systematics
   std::vector< const ISyst * > systs =  {};
   if (corrSysts) systs =  {&kAnaCalibrationSyst};


////////////////////////////////////////////////////////////
/// Fit
////////////////////////////////////////////////////////////
    /// And finally let's start fitting

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

    osc::IOscCalcAdjustable* calc = DefaultOscCalc();
    
    SystShifts auxShifts = SystShifts::Nominal();
    // In case some systs need different seeds
    std::vector <SystShifts> seedShifts = {};

//////////////////////////////////////////////////////////////////////
/////////////////////////// Seed controller ////////////////////////////
////////////////////////////////////////////////////////////////////////
   /// This section contains all seeds that will be used further

    struct th23helper{
      std::vector<double> seeds;
      const IFitVar * var;
      TString label;
     };

    std::vector <th23helper> th23seeds;
    //for NH:
    /// make two seeds in both UO and LO (near maxmixing and th23 from pre-seeding stage)
    th23seeds.push_back( { {0.499, 0.45}, kFitSinSqTheta23LowerOctant, "LO"});
    th23seeds.push_back( { {0.501, 0.55}, kFitSinSqTheta23UpperOctant, "UO"});

    /// Put three seeds - LO, UO, maxmixing
    std::vector<double> th23_all_seeds_NH = {0.45, 0.5, 0.55};
    /// delta CP seeds
    std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
    double seed_dmsq = 2.5e-3;

/////////////////////////////////////////////////////////////////////////
//////////////////////////// Different best fit searches //////////////////
///////////////////////////////////////////////////////////////////////////
    /// Now, let's find BF value, it searches BF in IH/NH and UO/LO and chooses the best chisq
    double minchi23 = 1E20;
    for(int hie:{-1,1}){
       /// will use th23seeds which we defined previously for LO and UO
       for (auto & thseed: th23seeds){
        std::cout << "\n\nFinding best fit "
                  << thseed.label
                  << (hie>0? " NH " : " IH ")
                  << "\n\n";
        for (auto const & s:thseed.seeds) std::cout << s << ", ";
        std::cout << std::endl;
        std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
                                                thseed.var,
                                                &kFitDmSq32Scaled,
                                                &kFitSinSq2Theta13};

        Fitter fit23(exptThis, fitvars, systs);
        auto thisminchi = fit23.Fit(calc,auxShifts ,{
                                   {&kFitDmSq32Scaled,{hie*seed_dmsq}},
                                   {thseed.var, thseed.seeds},
                                   {&kFitDeltaInPiUnits, delta_seeds},
                                   {&kFitSinSq2Theta13, {calc->GetTh13()}} },
                                   {}, Fitter::kVerbose
                                   );
        /// compare with global BF
        minchi23= min(minchi23, thisminchi);
        ResetOscCalcToDefault(calc);
        auxShifts.ResetToNominal();
      }//end th23
    }//end hie
    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");

    outfile->cd();
    TVectorD v(1);
    v[0] = minchi23;
    cout<<"minchi wrote into the file "<<minchi23<<endl;
    v.Write("overall_minchi");

    /// How much point on the grid do you want? It will be steps*steps
    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
    /// Generally need both IH and NH, but it's too long to run them all together
    for(int hie: {-1, +1}){

      if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;  /// Will run only one hie which you chose at the beginning

      std::cout << "Starting surface " << (hie>0? "NH ": "IH") << "\n\n";
      cout<<"the delta seeds are ";
      for (auto const & s:delta_seeds) std::cout << s << ", ";
      std::string hieStr = hie>0 ? "NH" : "IH";

      /// It will be deltaCP vs Th23 surface
      if(!th13Surf && ! dmsqSurf){
        ResetOscCalcToDefault(calc);
        calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
        /// Set Y axis's low and up values
        double low = (nueOnly||RHCOnly)?0:0.3;
        double high = (nueOnly||RHCOnly)?1:0.7;
        cout<<"low "<<low<<" high "<<high<<endl;
        /// It will minimize dmsq23 and th13 only, no seeds here
        Surface surf23(exptThis, calc,
                       &kFitDeltaInPiUnits,steps, 0, 2,       
                       &kFitSinSqTheta23,  steps, low, high, 
                       {&kFitDmSq32Scaled,&kFitSinSq2Theta13},systs//,
                       //{},seedShifts
                       );

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

      /// This part will make th23 vs dmsq surface
      if(!th13Surf && dmsqSurf){
        calc = DefaultOscCalc();
        calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
        /// Put Y axis's low and up values
        double low = (nueOnly||RHCOnly)?0.2:0.35;
        double high = (nueOnly||RHCOnly)?1:0.7;
        cout<<"low "<<low<<" high "<<high<<endl;
        /// It will minimize deltaCP and th13, use seeds for deltaCP
        /// also, there is some logic for x axis's values as well
        /// for different cases of fit (RHC/FHC/both) contours can be larger or smaller. To much pain in these values ...   
        Surface surf23m(exptThis, calc,
                        &kFitSinSqTheta23, steps,low, high,
                        &kFitDmSq32Scaled, steps,
                        (hie>0?(RHCOnly?2.0:2.2):(RHCOnly?-3.5:-2.8)),(hie>0?(RHCOnly?3.2:2.85):(RHCOnly?-2.0:-2.3)),
                        {&kFitSinSq2Theta13, &kFitDeltaInPiUnits},systs,
                        {{&kFitDeltaInPiUnits, delta_seeds}});
        auto surf6=surf23m.ToTH2(minchi23);
        outfile->cd();
        surf6->Write((TString)"th23_dm32_"+(hie>0?"NH":"IH"));
        surf23m.SaveTo(outfile->mkdir((TString)"surf_th23_dm32_"+(hie>0?"NH":"IH")));
      }

      /// Also, we can do th13 vs deltaCP contour, why not 
      /// (we don't show them to public since there are some issues with FC for th13 measurements)
      /// but gaussian fit is ok
      if(th13Surf){
        calc = DefaultOscCalc();
        calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
        /// Minimize dmsq23, th23, put seeds for th23 as (smth in LO, 0.5, smth in UO)
        Surface surf13(exptThis, calc,
                       &kFitSinSq2Theta13,  steps, 0, (RHCOnly?1:0.35),
                       &kFitDeltaInPiUnits,steps, 0, 2,
                       {&kFitSinSqTheta23, &kFitDmSq32Scaled}, systs,
                       {{&kFitSinSqTheta23, th23_all_seeds_NH}});//, seedShifts);
        auto surf4 = surf13.ToTH2(minchi23);
        outfile->cd();
        surf4->Write((TString)"th13_delta_"+(hie>0?"NH":"IH"));
        surf13.SaveTo(outfile->mkdir((TString)"surf_th13_delta_"+(hie>0?"NH":"IH")));

      }
    } //end hie

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

/// Done with fit! Now play with style and make pretty plots.

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

  else{ //plot only

    /// read file from the previous section
    TFile * infile = new TFile (outfilename+".root","read");
    /// Get BF
    auto mins =* (TVectorD*)infile->Get("overall_minchi");
    double minchi23 = mins[0];
    /// Some logic to select contour which you will plot
    std::vector <TString> surfNames;
    if(!th13Surf && ! dmsqSurf)surfNames = {"delta_th23_"};
    if (dmsqSurf) surfNames.push_back("th23_dm32_");
    if (th13Surf) surfNames.push_back("th13_delta_");
    for(TString surfName:surfNames){
    /// again, need two hierarchies, but all files were done with just one type
      for (int hie:{-1,+1}){

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

        //Setup canvas for each contour type, zoom if necessary 
        if(surfName.Contains("delta_th23")) {
          if (zoomIn) DrawContourCanvas(surfName, 0., 2., (RHCOnly?0.15:0.275), (RHCOnly?0.8:0.725));
          else DrawContourCanvas(surfName, 0., 2., 0., 1.);
        }
        if(surfName.Contains("th23_dm32")) {
          if (hie>0) DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?1.8:2.05), (RHCOnly?3.3:2.85));
          else DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?-3.5:-2.95), (RHCOnly?-2.0:-2.25));
        }
        if(surfName.Contains("th13_delta")) {
          if (zoomIn) DrawContourCanvas(surfName, 0., 0.38, 0, 2);
          else DrawContourCanvas(surfName, 0., 1.0, 0, 2);
         }
        /// Read your surface and set 1,2,3 sigma curves 
        auto surf = *Surface::LoadFrom(infile->GetDirectory("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);
        /// colors
        Color_t k2SigmaConfidenceColor = ((hie>0)? k2SigmaConfidenceColorNH:k2SigmaConfidenceColorIH);
        Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);
        Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.38);
        Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.18);
        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"));
        }
        else{
          cout<<"minchi from file is "<<minchi23<<endl;
            
            surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
            surf.DrawContour(surf_2Sigma, k2SigmaConfidenceStyle, kFillColor2,minchi23);
            surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
        }

        gPad->RedrawAxis();
   
        /// Set legend
        auto leg = ContourLegend(hie, fillContour, false,
                                 kFillColor1,kFillColor2,kFillColor3,kDarkColor,
                                 false);
        /// Put PDG constraint line on the th13 vs deltaCP surface
        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);
          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);
       
        /// and finally (!) save the plot to the file
        for(TString ext: {".pdf",".root"}){
          gPad->Print(outdir + "contour_" +  suffix + "_" +
                      surfName +
                      (hie>0?"NH":"IH")  +
                      (zoomIn?"_zoom":"") +
                      ext);
        }//end ext

        /// in case if you want (of course you want) to see debug plots for each syst or osc parameter put true here
        auto debug = true; 
        if(debug){
          auto hists = surf.GetProfiledHists();
          int nSysts = corrSysts ? 51 : 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(outdir + "debug_contour_" +  suffix + "_" +
                        surfName +
                        (hie > 0 ? "NH":"IH")  +
                        (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}");
  }
}

/// This one needs some work by hands with putting JoinGraphs properly for each contour
/// I didn't make it for the tutorial cases
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("c f");
        JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("c f");
        JoinGraphs(g2[0], g2[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");}
      }
    }

    if(isJoint && !isNH){
      if(isRHC)  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
      else{
          JoinGraphs(g3[0], g3[2], kFillColor3)->Draw("f c");
          g3[1]->SetFillColor(kFillColor3); g3[1]->Draw("f c");
          g3[3]->SetFillColor(kFillColor3); g3[3]->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();
}


std::vector <Spectrum * > GetMockData(TString ana_options){

  cout<<"\n\n ============= You're loading Mock data from the file ============="<<endl;
  std::vector <Spectrum *> spec;
  std::string Dir = "./";
  bool RHCmode = ana_options.Contains("rhc");
  bool FHCmode = ana_options.Contains("fhc");
  bool NueData = ana_options.Contains("nue");
  bool NumuData = ana_options.Contains("numu");
  if (NueData){
     if (FHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_0").release());
     if (RHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_1").release());
  }
  if (NumuData){
     if (FHCmode){
        for (int i = 2; i < 6; ++i) {
             auto temp = LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_"+std::to_string(i)).release();
             spec.push_back(temp);
        }
     }
     if (RHCmode){
        for (int i = 6; i < 10; ++i) {
             auto temp = LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_"+std::to_string(i)).release();
             spec.push_back(temp);
        }
     }

  }
  return spec;

}