Tutorial2019FitSlices.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 "OscLib/IOscCalc.h"
#include "Utilities/rootlogon.C"

#include "CAFAna/nue/Ana2019/joint_fit_2019_loader_tools.h"

#include "TCanvas.h"
#include "TBox.h"
#include "TColor.h"
#include "TGraph.h"
#include "TVectorD.h"
#include "TF1.h"
#include "TLegend.h"
#include "TText.h"
#include "TLatex.h"
#include "TPad.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 POTtag(bool both, bool FHCOnly, bool RHCOnly);

void MakeMaps ( std::vector <const IFitVar*> profvars,
                std::map<const IFitVar*, TGraph*> &profMap, 
                std::vector <const ISyst* > systs,
                std::map<const ISyst*, TGraph*> &systMap);

void SaveMaps ( TDirectory * dir, 
                TString prefix, 
                std::vector <const IFitVar*> profvars,
                std::vector <TString > profnames,
                std::map<const IFitVar*, TGraph*> &profMap, 
                std::vector <const ISyst* > systs,
                std::map<const ISyst*, TGraph*> &systMap);


void Tutorial2019FitSlices(bool createFile = false,         /// Do you want to make fit and write it into the file or just plot existing slice
                           bool corrSysts = true,           /// syst or stat only?
                           TString options="joint_realData_both",    /// what kind of analysis? numuOnly||nueOnly||joint ana for RHCOnly||RHCOnly||both modes
                           bool th23slice = false,          /// is it th23 slice?
                           bool octantSlice = true,         /// dou you want x2 curves split by octant?
                           bool dmsqSlice = false,          /// is it dmsq32 slice?
                           bool th13Slice = false)          /// is it th13 slice?
                                                            /// in case of th23slice = false && dmsqSlice = false && th13Slice = false it will be deltaCP slice

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

  auto suffix = options;
  if(corrSysts) suffix+="_systs";
  else suffix+="_stats";

/*  TString tag;
  if(FHCOnly) tag="FHCOnly/";
  if(RHCOnly) tag="RHCOnly/";
  if(both) tag = "RHC_and_FHC/";*/

  if(th23slice) suffix+="_th23";
  if(dmsqSlice) suffix+="_dmsq";
  if(th13Slice) suffix+="_th13";
  if(!octantSlice)suffix+="_noOct";

  TString outdir = "";
  TString outfilename (outdir +"hist_slices_2019_" + suffix);
  TString outFCfilename (outdir + "slices_FCInput_2019_" + suffix);   /// need for syst. debug profiles

  if(createFile){
    
    //////////////////////////////////////////////////
    // Load Nue and Numu experiments
    //////////////////////////////////////////////////
    //need numu only for prestage seeds

    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)});
        }
       }
    }
   
    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
    ////////////////////////////////////////////////////////////
    std::cout << "\nCreating multiexperiment\n" << std::endl;

    if (!th13Slice) {expts.push_back(WorldReactorConstraint2017()); std::cout << "Adding WorldReactorConstraint2017\n";}
    else {std::cout << "No reactor constraint, that's th13 slice\n";}
    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";

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

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

    osc::IOscCalcAdjustable* calc = DefaultOscCalc();

    SystShifts auxShifts = SystShifts::Nominal();

    // In case some systs need different seeds
    // If you run ana script for 2019 you will see seeds for some systematics for paricular plots
    // Comment for the workshop, but leave here for the education purposes
    
    std::vector <SystShifts> seedShifts = {};
    /*if(corrSysts && !th23slice && !dmsqSlice && !th13Slice && octantSlice && (both || FHCOnly)){   //Cherenkov seed for deltaCP and neutron as well :(
      cout<<"Add Cherenkov seeds \n";
      for (double systshift:{-0.5, +0.5}){                                     //for the FHC Only or FHC+RHC ana
        SystShifts tempShifts (&kAnaCherenkovSyst, systshift);
        seedShifts.push_back(tempShifts);
      }
      if(both){
        cout<<"Add neutron systematic seeds \n";                               //For the joint RHC+FHC ana only
        for (double systshift:{-0.5, +0.5}){
        SystShifts tempShifts (&kNeutronVisEScalePrimariesSyst2018, systshift);
        seedShifts.push_back(tempShifts);
        }
      }
    }
    */


//////////////////////////////////////////////////////////////////////
///////////////////////// Seed controller ////////////////////////////
//////////////////////////////////////////////////////////////////////

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

   std::vector <th23helper> th23seeds;

  if(octantSlice) {
    th23seeds.push_back( { {0.499, 0.45}, kFitSinSqTheta23LowerOctant, "LO"});
    th23seeds.push_back( { {0.501, 0.55}, kFitSinSqTheta23UpperOctant, "UO"});
   }
   else th23seeds.push_back({ {0.45, 0.5, 0.55}, &kFitSinSqTheta23, ""});

   std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
   std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
   double seed_dmsq = 2.5e-3;

/////////////////////////////////////////////////////////////////////////
////////////////////////// Different best fit searches //////////////////
/////////////////////////////////////////////////////////////////////////

    // Find the best fit points
    double minchi23 = 1E20;
    for(int hie:{-1, 1}){
      for (auto & thseed:th23seeds){

        std::cout << "\n\nFinding best fit " << (hie > 0 ? "NH " : "IH ")
                  << thseed.label << ", "
                  << "ssth23 seeds ";
        for (auto const & s:thseed.seeds) std::cout << s << ", ";
        std::cout << std::endl;

        std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
                                                thseed.var, 
                                                &kFitDmSq32, 
                                                &kFitSinSq2Theta13};
        
        Fitter fit23(exptThis, fitvars, systs);
        auto thisminchi = fit23.Fit(calc, auxShifts,
                                    {
                                      { &kFitDmSq32,   {hie*seed_dmsq} },
                                      { thseed.var,          thseed.seeds },
                                      { &kFitDeltaInPiUnits, delta_seeds }
                                    }, seedShifts, 
                                    Fitter::kVerbose);
        minchi23= min(minchi23, thisminchi);

        ResetOscCalcToDefault(calc);
        auxShifts.ResetToNominal();
      }
    }
    std::cout << "\nFound overall minchi " << minchi23 << "\n\n";

/////////////////////////////////////////////////////////////////////////////
/////////////////////// Now, slices! ////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////

    ///We'll save to this file, make pretty plots later
    TFile * outfile = new TFile(outfilename+".root","recreate");
    TFile * outFCfile = new TFile(outFCfilename+".root","recreate"); 
    /// again, save BF into the file
    outfile->cd();
    TVectorD v(1);
    v[0] = minchi23;
    v.Write("overall_minchi");
    // how many points do you want?
    int steps = 60;// 80;

    ///Now create all the slices
    for(int hie: {-1, +1}){

      ResetOscCalcToDefault(calc);
      calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));

      std::cout << "Starting profile " << (hie>0? "NH ": "IH") << "\n\n";
      ///Structure defining the slice 
      struct ProfDef{
        const IFitVar * fitvar;   // which var to fix in each point? Or what kind of slice it is?
        double minX;
        double maxX;
        std::vector <const IFitVar * > profvars; /// which vars to profile?
        std::vector <TString> profvarnames;      ///names of these vars
        std::map  <const IFitVar *, std::vector <double> >profseeds; /// seeds for osc. params
        TString shortname;
      };

       /// will loop once if you selected w/o octant split and twice if you want LO and UO curves
      for(auto const & thseed:th23seeds){  
        ProfDef profdef;
        const IFitVar* kCurTheta23 =  thseed.var; 
         
        /// deltaCP slice
        if(!th23slice && !dmsqSlice && !th13Slice){ 
          profdef.fitvar =  &kFitDeltaInPiUnits;
          profdef.minX = 0;
          profdef.maxX = 2;
          profdef.profvars =  {&kFitDmSq32, &kFitSinSq2Theta13, kCurTheta23};  ///takes th23 vars from th23seeds (defined in seed controller) - it can be th23 var constrained in LO or UO, or ordinary th23
          profdef.profvarnames = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
          profdef.profseeds = {{kCurTheta23, thseed.seeds}, { &kFitDmSq32, {hie*seed_dmsq}} };  /// th23 seed from the seed controller
          profdef.shortname="delta";
        }
 
        ///th23 slice - no split into octants   
        if(th23slice){
          profdef.fitvar = &kFitSinSqTheta23;
          profdef.minX = (nueOnly ? 0:0.3);
          profdef.maxX = (nueOnly ? 1:0.7);
          profdef.profvars =  {&kFitDmSq32, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
          profdef.profvarnames =  {"DmSq32","SinSq2Theta13","DeltaCPpi"};
          profdef.profseeds = { {&kFitDeltaInPiUnits,delta_seeds}};
          profdef.shortname="th23";
        }
        ///th13 slice
        if(th13Slice){
          profdef.fitvar = &kFitSinSq2Theta13;
          profdef.minX = 0.01;
          profdef.maxX = 0.2;
          profdef.profvars =  {&kFitDmSq32, kCurTheta23, &kFitDeltaInPiUnits};
          profdef.profvarnames =  {"DmSq32","SinSqTheta23","DeltaCPpi"};
          profdef.profseeds = { {kCurTheta23, thseed.seeds}, {&kFitDeltaInPiUnits,delta_seeds}, { &kFitDmSq32, {hie*seed_dmsq}} };
          profdef.shortname="th13";
        }
        ///dmsq32 slice
        if(dmsqSlice){
          profdef.fitvar = &kFitDmSq32;
          profdef.minX = hie > 0 ? 2e-3 : -3e-3;
          profdef.maxX = hie > 0 ? 3e-3 : -2e-3;
          profdef.profvars =   {kCurTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
          profdef.profvarnames = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
          profdef.profseeds = { {kCurTheta23, thseed.seeds},
                                {&kFitDeltaInPiUnits,delta_seeds}};
          profdef.shortname="dmsq";
        }
        // these three lines are arguments for SqrtProfile, they are used for FC
        std::map<const IFitVar*, TGraph*> profVarsMap;   
        std::map<const ISyst*, TGraph*> systMap;
        MakeMaps (profdef.profvars, profVarsMap, systs, systMap);

        
        std::cout << " Profile "<< thseed.label << "\n";
        ///make a profile       
        auto slice = SqrtProfile(exptThis, calc, 
                                 profdef.fitvar, steps, profdef.minX, profdef.maxX,
                                 minchi23,
                                 profdef.profvars,
                                 systs, 
                                 profdef.profseeds,
                                 seedShifts,
                                 profVarsMap, systMap);
        //Save!
        TString hieroctstr = (hie > 0 ? "NH" : "IH") + thseed.label;
        //save to the file
        outfile->cd();
        slice->Write((TString ) "slice_" + profdef.shortname + "_" + hieroctstr);
        SaveMaps (outFCfile, hieroctstr, 
                  profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
                
      }//end seeds
    }//end hie
  
    delete outfile;
    std::cout << "\nProfiles saved to " << outfilename << std::endl;
    return;
  }//end createfile


///Done with fitting!
/// Now style work

  //////////////////////////////////////////////////////////////////////
  // Make plots
  //////////////////////////////////////////////////////////////////////
  else{ 
    std::cout << "Making plots from " << outfilename+".root" << std::endl;
    TFile * infile = new TFile (outfilename+".root","read");

     //read minichi
    auto mins =* (TVectorD*)infile->Get("overall_minchi");
    double minchi23 = mins[0];
    
    std::vector <TString> sliceNames = {};
    if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {"slice_delta_"};
    if (th23slice) sliceNames = {"slice_th23_"};
    if (dmsqSlice) sliceNames = {"slice_dmsq_"};
    if (th13Slice) sliceNames = {"slice_th13_"};

    // it shouldn't play any role, but copy old style seeding here for sake of naming, everywhere is just .lable
   
    struct th23helper{
      std::vector<double> seeds;
      const IFitVar * var;
      TString label;
    };
    std::vector <th23helper> th23seeds;
    if(octantSlice) {
       th23seeds.push_back( { {0.45}, kFitSinSqTheta23LowerOctant, "LO"});
       th23seeds.push_back( { {0.55}, kFitSinSqTheta23UpperOctant, "UO"});
     }
    else th23seeds.push_back({ {0.45, 0.5, 0.55}, &kFitSinSqTheta23, ""});

   //end of old-style-seeding

    /// loop over all slices in this file (IH, NH, LO, UO)
    for(TString sliceName:sliceNames){
      std::vector < std::pair <TGraph*, TString> > graphs;
      
      if(sliceName.Contains("delta")) {
        for (TString hie:{"NH","IH"}){
          for(auto const & thseed:th23seeds){
            auto h = (TGraph*) infile ->Get( sliceName + hie + thseed.label);
            if(!h) {
              std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; 
              exit(1);
            }
            graphs.push_back({h, hie + " " + thseed.label}); 
          }
        }
        DrawSliceCanvas(sliceName, 7, 0, 2);
      }
      else {
        for (TString hie:{"NH","IH"}) {
          for(auto const & thseed:th23seeds){
            auto h = (TGraph*) infile ->Get( sliceName + hie + thseed.label);
            if(!h) {std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; }
            graphs.push_back({h,hie + " " + thseed.label});
          }
        }
        
        if(sliceName.Contains("th23")) {
          DrawSliceCanvas(sliceName, 5, 
                          (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
        }
        
        if(sliceName.Contains("dmsq")) {
          DrawSliceCanvas(sliceName, 5, 2.15, 2.75);
        }
        if(sliceName.Contains("th13")) {
          DrawSliceCanvas(sliceName, 5, 0.05, 0.2);
        }
      }

      //loop over all slices and set colors
      for (auto &gr:graphs) {
        if(gr.second.Contains("IH")) gr.first->SetLineColor(kPrimColorIH);
        if(gr.second.Contains("NH")) gr.first->SetLineColor(k3SigmaConfidenceColorNH);
        if(gr.second.Contains("LO")) gr.first->SetLineStyle(7);
        gr.first->SetLineWidth(3);
        if(dmsqSlice){ 
          for (int i =0; i < gr.first->GetN(); ++i){
            gr.first->SetPoint(i,1000*fabs(gr.first->GetX()[i]),gr.first->GetY()[i]);
          }
        }
        gr.first->Draw("l same");
      }
      
      //POTtag(both, FHCOnly, RHCOnly); 
      PreliminarySide();
      gPad->RedrawAxis();
      
      auto leg = SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
      leg->Draw();
     
      //      outdir = "test_delta_slice/";
      //gPad->SetGrid();
      for(TString ext: {".pdf",".root"}) {
        gPad->Print(outdir + "slice_" +  suffix +
                    ext);
      }//end ext

    }//end slicename

    
    bool debug = true;

    if(debug){

      TFile * inFCfile = new TFile(outFCfilename+".root","read");
      std::cout << "\n\nOpen profiled " << outFCfilename+".root\n" << std::endl;

      std::vector <TString> strprof = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
      if(th23slice) strprof = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
      if(th13Slice) strprof = {"DmSq32","DeltaCPpi", "SinSqTheta23"};
      if(dmsqSlice) strprof = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};

      std::vector< const ISyst * > systs =  {};
      if (corrSysts){      // syst case
          systs =  {&kAnaCalibrationSyst}; // replace the only one syst
      }
      
      for (const ISyst* s : systs){
        strprof.push_back(s->ShortName());
      }
      
      for(auto const & str:strprof){
        new TCanvas;
        double miny=-1.5, maxy=1.5;
        if(str=="DeltaCPpi")     {miny=0; maxy=2;}
        if(str=="DmSq32")        {miny=2.2; maxy=2.9;}
        if(str=="SinSq2Theta13") {miny=0.080; maxy=0.086;}
        if(str=="SinSqTheta23")  {miny=0.3; maxy=0.8;}
        
        auto axes = new TH2F("",";#delta_{CP}/#pi;" + str, 100, 0, 2, 100, miny, maxy);
        if(dmsqSlice) axes = new TH2F("",";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" + str, 100, 2.2, 2.8, 100, miny, maxy);
        if(th23slice) axes = new TH2F("",";sin^{2}#theta_{23};" + str, 100, 0.3, 0.7, 100, miny, maxy);
        if(th13Slice) axes = new TH2F("",";sin^{2}#theta_{13};" + str, 100, 0.01, 0.2, 100, miny, maxy);
        axes->Draw("hist");

        for (TString hie:{"NH","IH"}){
          for(auto const & thseed:th23seeds){
            auto gr = (TGraph*)  inFCfile->Get(hie + thseed.label + "_" + str);
            if(!gr) {
              std::cerr << "Problem " << hie + thseed.label + "_" + str << "\n"; 
              exit(1);
            }
            if(str=="DmSq32"){
              for (int i = 0; i < (int) gr->GetN(); ++i){
                gr->SetPoint(i, gr->GetX()[i], 1000* fabs(gr->GetY()[i]));
              }
            }
            if(dmsqSlice){
              for (int i = 0; i < (int) gr->GetN(); ++i){
                gr->SetPoint(i, 1000*fabs(gr->GetX()[i]), gr->GetY()[i]);
              }
            }
            gr->SetLineColor(hie=="NH"?k3SigmaConfidenceColorNH:kPrimColorIH);
            if(thseed.label.Contains("LO")) gr->SetLineStyle(7);
            gr->SetLineWidth(3);
            gr->SetMarkerColor(gr->GetLineColor());
            gr->Draw("l same");
          }//end th23
        }//end hie
        gPad->SetGrid();
        gPad->Print(outdir + "debug_slice_prof_" + suffix + "_"+ str + ".pdf");
      }//end profile vars
    }//end debug
    
  }//end plots
}//end all


void MakeMaps ( std::vector <const IFitVar*> profvars,
                std::map<const IFitVar*, TGraph*> &profMap,
                std::vector <const ISyst* > systs,
                std::map<const ISyst*, TGraph*> &systMap)
{
  for (const IFitVar * var:profvars)
              profMap.insert(std::pair<const IFitVar*, TGraph*> (var, new TGraph()));
  for (const ISyst* syst : systs)
    systMap.insert(std::pair<const ISyst*,   TGraph*> (syst, new TGraph()));
}
void SaveMaps ( TDirectory * dir,
                TString prefix,
                std::vector <const IFitVar*> profvars,
                std::vector <TString > profnames,
                std::map<const IFitVar*, TGraph*> &profMap,
                std::vector <const ISyst* > systs,
                std::map<const ISyst*, TGraph*> &systMap)
{
  TDirectory *tmp = gDirectory;
  dir->cd();
  for (int i = 0; i < (int) profvars.size(); ++i){
    profMap[profvars[i]]->Write((prefix + "_" + profnames[i]));
  }
  for (const ISyst* s : systs)
    systMap[s]->Write((prefix + "_" + s->ShortName()));
  tmp->cd();
}


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;

}