demoFitContours.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 "CAFAna/nue/Ana2018/joint_fit_2018_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);     
  }

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

void demoFitContours(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_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 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_2018_"    + 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)

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

    auto calc_fake = DefaultOscCalc();
    SetFakeCalc(calc_fake, 0.58, 2.51e-3, 0.17); /// use 2018 BF for making fake data

    /// And fill all vectors we created, the structure is: 
    // nue FHC, nue RHC, 4 numu FHC quartiles, 4 numu RHC quartiles - for the joint_both option selected
    // If you selected numuOnly|nueOnly|RHCOnly|FHCOnly then some of these predictions will be switched off

    if(!numuOnly) { /// fill nue part
      if(FHCOnly || both ) {
         preds.push_back(GetNuePrediction2018("combo", DefaultOscCalc(), corrSysts, "fhc", false));  // nue FHC
         cosmics.push_back(GetNueCosmics2018("fhc"));
      }
      if(RHCOnly || both ) {
         preds.push_back(GetNuePrediction2018("prop", DefaultOscCalc(), corrSysts, "rhc", false));  // nue FHC
         cosmics.push_back(GetNueCosmics2018("rhc"));
      } 
    }

    int nnumu = 4;
    if(!nueOnly) {  /// fill numu part and numu_only for preseeds
      if(FHCOnly || both ) {
      
        auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "fhc");     // numu FHC
        preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
        preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());

        auto numu_cosmics = GetNumuCosmics2018(nnumu, "fhc");
        cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
        cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
  
      }
      if(RHCOnly || both ) {
      
        auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "rhc");   // numu RHC
        preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
        preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());

        auto numu_cosmics = GetNumuCosmics2018(nnumu, "rhc");
        cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
        cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
       }
    }
    
    /// 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";

    for(int i = 0; i < int(preds.size()); ++i){
      double POT;
      if(FHCOnly) POT = kAna2018FHCPOT;
      if(RHCOnly) POT = kAna2018RHCPOT;
      if(both) {
        if(nueOnly || both) {
           if(i==0) POT = kAna2018FHCPOT;
           if(i==1) POT = kAna2018RHCPOT;
        }
        if(numuOnly) {
           if(i >= 0 && i < 4) POT = kAna2018FHCPOT;
           if(i >= 4 && i < 8) POT = kAna2018RHCPOT;
         }
         if(both) {
           if(i >= 2 && i < 6) POT = kAna2018FHCPOT;
           if(i >= 6 && i < 10) POT = kAna2018RHCPOT;
        }
      }
      /// done with POT, let's make fake data
      auto thisdata = GetFakeData (preds[i], calc_fake, POT, cosmics[i].first);
      cout<<i<<" POT "<<POT<<" tot MC "<<preds[i]->Predict(calc_fake).Integral(POT)<<" cos "<<cosmics[i].first->Integral()<<" cos er "<<cosmics[i].second<<" analyze data "<<thisdata->Integral(POT)<<endl;
      expts.push_back(new SingleSampleExperiment(preds[i],*thisdata, cosmics[i].first,cosmics[i].second)); // make an experiment with current predictions and data
      /// for debug purposes print spectra to the file
      auto debug=true;
      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 + (corrSysts?"syst/":"stat/")+ "debug_predictions_" + suffix + "_" + std::to_string(i) + ".pdf");
        gPad->Print(outdir +  "debug_predictions_" + suffix + "_" + std::to_string(i) + ".pdf");
      }//end debug
    }
     
    /// Out main preporation part is done - we have 10 experiments for fit
    /// one more step - do the same but with numu only predictions for pre-seeding 

    cout<<"\n\nMake numu only experiment to get the seeds"<<endl;
    //Follow the same scheme but with reduced number of predictions:
      for(int i = 0; i < (int) preds_numu_only.size(); ++i){
      double POT;
      if(FHCOnly) POT = kAna2018FHCPOT;
      if(RHCOnly) POT = kAna2018RHCPOT;
      if(both) {
           if(i >= 0 && i < 4) POT = kAna2018FHCPOT;
           if(i >= 4 && i < 8) POT = kAna2018RHCPOT;
      }
      auto thisdata = GetFakeData (preds_numu_only[i],calc_fake, POT, cosmics_numu_only[i].first);
       cout<<i<<" POT "<<POT<<" tot MC "<<preds_numu_only[i]->Predict(calc_fake).Integral(POT)<<" cos "<<cosmics_numu_only[i].first->Integral()<<" cos er "<<cosmics_numu_only[i].second<<" analyze data "<<thisdata->Integral(POT)<<endl;
      expts_numu_only.push_back(new SingleSampleExperiment(preds_numu_only[i],*thisdata, cosmics_numu_only[i].first,cosmics_numu_only[i].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);
    /// and multiexperiment for numu pre-seeds
    std::cout << "Creating Multiexperiment of numu only SimpleExp with a total of "
              << expts_numu_only.size() << " experiments\n\n" << std::endl;
    auto exptThis_numu_only = new MultiExperiment(expts_numu_only);

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

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

   /// Now, let's prepare systematics
   std::vector< const ISyst * > systs =  {};
   std::vector< const ISyst * > systs_numu_only =  {};
   if (corrSysts){      // syst case
       systs =  {&kAnaCalibrationSyst};                //here will be your own syst, which you made previously
       systs_numu_only =  {&kAnaCalibrationSyst};      //here will be your own syst, which you made previously
   }

    /// And that's how it looks like in the analysis - you need to set proper correlations for all simple experiments

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

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

////////////////////////////////////////////////////////////
/// 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
   /// There is numu pre-seed part here as well. Idea is to choose th23 and dmsq23 seeds
   /// from the numu only fast fit. It is more important for FC, but for similarity we put this
   /// part in gaussian fit as well

    cout<<"------------------- Start prestage seeds --------------------------"<<endl;

    double minchi_numu = 1E20;
    double pre_seed_th23;
    double pre_seed_dmsq;
    ResetOscCalcToDefault(calc);
    auxShifts.ResetToNominal();

    double maxmix = 0.514;  // from the numu results
    double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
    /// it simply searches BF in NH and use BF th23 and dmsq values as seeds
    for(int hie:{1}){
        std::cout << "\n\nFinding test best fit " << (hie>0? "NH " : "IH ") << "\n\n";
        Fitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32Scaled}, {});

        auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
                                        {{&kFitDmSq32Scaled,{hie*2.5}},
                                        {&kFitSinSqTheta23, {0.4} }}, {},// seedShifts
                                        Fitter::kVerbose);

        if(thisminchi < minchi_numu) minchi_numu = thisminchi;
    }
    /// our BF values
    pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
    pre_seed_dmsq = kFitDmSq32Scaled.GetValue(calc);
    /// make two seeds in LO and UO with the help of pre_seed_dmsq
    numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
    numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));

    ResetOscCalcToDefault(calc);
    auxShifts.ResetToNominal();

    cout<<"------------------- End prestage seeds --------------------------"<<endl;

    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, numu_pre_seedLONH}, kFitSinSqTheta23LowerOctant, "LO"});
    th23seeds.push_back( { {0.501, numu_pre_seedUONH}, kFitSinSqTheta23UpperOctant, "UO"});
    /// Put three seeds - LO, UO, maxmixing
    std::vector<double> th23_all_seeds_NH = {numu_pre_seedLONH, 0.5, numu_pre_seedUONH};
    /// delta CP seeds
    std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};

/////////////////////////////////////////////////////////////////////////
//////////////////////////// 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;
        cout<<"pre seed from numu is "<<pre_seed_dmsq<<endl;
        /// replace th23 variable by kFitSinSqTheta23LowerOctant/kFitSinSqTheta23UpperOctant
        // We minimize deltaCP, th23, dmsq23, th13 (which is constrained)
        std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
                                                thseed.var,
                                                &kFitDmSq32Scaled,
                                                &kFitSinSq2Theta13};

        Fitter fit23(exptThis, fitvars, systs);
        /// Minimize and use seeds for dmsq (the one from pre-seed stage * 1 or -1 for NH and IH)
        /// th23 will use seeds from prestage for kFitSinSqTheta23LowerOctant/kFitSinSqTheta23UpperOctant
        /// th13 seed is PDG's central value
        /// deltaCP seeds at 0, pi/2, pi, 3pi/2
        auto thisminchi = fit23.Fit(calc,auxShifts ,{
                                   {&kFitDmSq32Scaled,{hie*pre_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
        FrequentistSurface 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 ...   
        FrequentistSurface 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)
        FrequentistSurface 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 will you plot
    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){
    /// again, need two hierarchies, but all files were done with just one
      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 = *FrequentistSurface::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 kFCFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);
        Int_t kFCFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.38);
        Int_t kFCFillColor3 = 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();
   
        /// Make pretty point at BF value
        bool bestfit = 1;
        double  delta, th23, dmsq32, th13;
        if(bestfit && both && corrSysts && hie>0){
           if(!th13Surf){
             delta = 0.16576;
             th23 = 0.584509;
             dmsq32 = 2.50515;
             th13 = 0.0824663; 
           }
           else{delta = 0.678516; th13 = 0.0957894;}
           TMarker bf;
           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();
           }
        }
        /// Set legend
        auto leg = ContourLegend(hie, fillContour, false,
                                 kFillColor1,kFillColor2,kFillColor3,kDarkColor,
                                 bestfit && both);
        /// Put PDG's 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}");
  }
}


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