mre_blessed.C

Go to the documentation of this file.

#include <fstream>
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Cuts/SpillCuts.h"

using namespace ana;

#include "OscLib/IOscCalc.h"
#include "OscLib/OscCalcDumb.h"
#include "OscLib/OscCalcPMNSOpt.h"

#include "TCanvas.h"
#include "TObject.h"
#include "TArrow.h"
#include "TColor.h"
#include "TFile.h"
#include "TLine.h"
#include "TH1.h"
#include "TText.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TPad.h"
#include "TGraphAsymmErrors.h"

void Prelim()
{
   TLatex* prelim = new TLatex(.93, .9, "NOvA Preliminary");
   prelim->SetTextColor(kBlue);
   prelim->SetNDC();
   prelim->SetTextSize(2/30.);
   prelim->SetTextAngle(270);
   prelim->SetTextAlign(12);
   prelim->Draw();
}

TCanvas *MakeMePlot(std::vector<TH1*> Histos, std::vector<TString> Names, bool isFHC, string yname = "none"){

    TCanvas *c = new TCanvas("c", "canvas", 800, 700);
    TPad* p1 = new TPad("", "", 0, 0, 1, 1);
    TPad* p2 = new TPad("", "", 0, 0, 1, 1);
    p1->SetBottomMargin(.3);
    p2->SetTopMargin(.7);

    for(auto p:{p1,p2}){
        p->SetLeftMargin(0.1);
        p->SetFillStyle(0);
        p->Draw();
     }
     p1->cd();

     auto denom = (TH1*)Histos[0]->Clone();
     denom->SetMaximum(denom->GetMaximum()*1.4);
     denom->SetMinimum(0);
     denom->GetYaxis()->CenterTitle();
     denom->GetXaxis()->CenterTitle();
     denom->GetXaxis()->SetLabelOffset(999);
     if(yname !="none") denom->GetYaxis()->SetTitle(yname.c_str());

     denom->Draw("hist E1");
     for(int i=1; i<Histos.size(); i++){
         Histos[i]->Draw("hist E1 same");
     }
     denom->Draw("hist E1 same");
     p1->RedrawAxis();
     gPad->Update();
     TText *tt;
     if(isFHC){
        auto *tt = new TText(gPad->GetUxmax()*0.7, 1.02*gPad->GetUymax(), "FHC beam");
        tt->SetTextColor(kBlue);
        tt->SetTextSize(0.05);
        tt->Draw();}
      else{
        auto *tt = new TText(gPad->GetUxmax()*0.7, 1.02*gPad->GetUymax(), "RHC beam");
        tt->SetTextColor(kBlue);
        tt->SetTextSize(0.05);
        tt->Draw();}

      TLegend *leg = new TLegend(0.52,0.76,0.78,0.87);
      leg->SetTextSize(0.025);
      for(int i=0; i<Histos.size(); i++) leg->AddEntry(Histos[i], Names[i], "l");
      leg->Draw();

      gPad->RedrawAxis();
      p2->cd();
      auto temp = (TH1*)Histos[1]->Clone();
      temp->Divide(denom);
      temp->Draw();
      for(int i=2; i<Histos.size(); i++){
          auto temp = (TH1*)Histos[i]->Clone();
          temp->Divide(denom);
          temp->Draw("same");
       }
       temp->GetYaxis()->SetRangeUser(0.8, 1.25);
       temp->GetYaxis()->SetTitle("Ratio");
       p2->RedrawAxis();

       return c;

}


void mre_blessed(bool plot = true, bool isFHC = true, TString drawing_options = "red"){
    
    const std::string mode_tag = isFHC?"fhc":"rhc";
    std::string file_mc = "mre_pred_nd_cuts_"+mode_tag+"_mc.root";
    
    const Var kSlcE([](const caf::SRProxy *sr){
        return sr->slc.calE;
    });

    const Var kCalE1erProng([](const caf::SRProxy *sr){
        if ( !sr->vtx.elastic.IsValid ) return -1000.f;
        if (sr->vtx.elastic.fuzzyk.npng < 1) return -1000.f;
        return float(sr->vtx.elastic.fuzzyk.png[0].calE);});

    const Var kHadCalENo1erProng = kSlcE - kCalE1erProng; 

    const int kNumVars = 4;
   
    struct HistDef
    {
        std::string name;
        HistAxis axis;        
        int lowedge;
        int highedge;
        int rebin;
    };
 
     const HistDef defs[kNumVars] = {
        {"energy", {"#nu_{e} Energy (GeV)", Binning::Simple(120,0,12), kNueEnergy2018}, 0, 45, 1},
        {"caloE", {"Calorimetric energy (GeV)", Binning::Simple(120,0,12), kSlcE}, 0, 45, 1},
        {"caloE_1pr", {"Calorimetric energy, 1st prong", Binning::Simple(120,0,12), kCalE1erProng}, 0, 45, 1},
        {"calE_no1pr", {"Calorimetric energy, not 1st prong", Binning::Simple(120,0,12), kHadCalENo1erProng}, 0, 45, 1},
    };
    

    const Cut kNue2018FHCCVNCutLowBin = kCVNSSe >= kNue2018CVNFHCLowEdge && kCVNSSe < kNue2018CVNFHCHighEdge;
    const Cut kNue2018RHCCVNCutLowBin = kCVNSSe >= kNue2018CVNRHCLowEdge && kCVNSSe < kNue2018CVNRHCHighEdge;
    const Cut kNue2018FHCCVNCutHighBin = kCVNSSe >= kNue2018CVNFHCHighEdge;
    const Cut kNue2018RHCCVNCutHighBin = kCVNSSe >= kNue2018CVNRHCHighEdge;

    const Cut kNue2018CVNCutLowBin([kNue2018FHCCVNCutLowBin, kNue2018RHCCVNCutLowBin](const caf::SRProxy* sr) {
          if(kIsRHC(sr)) return kNue2018RHCCVNCutLowBin(sr);
          else return kNue2018FHCCVNCutLowBin(sr);
      });

    const Cut kNue2018CVNCutHighBin([kNue2018FHCCVNCutHighBin, kNue2018RHCCVNCutHighBin](const caf::SRProxy* sr) {
          if(kIsRHC(sr)) return kNue2018RHCCVNCutHighBin(sr);
          else return kNue2018FHCCVNCutHighBin(sr);
      });

    const  int kNumSels = 6;

    const Cut sels[kNumSels] = {
        kNue2018NDPresel && kNue2017MRParentSliceCut,
        kNue2018NDCVNSsb && kNue2017MRParentSliceCut,
        kNue2018NDPresel && kNue2017MRParentSliceCut,
        kNue2018NDPresel && kNue2017MRParentSliceCut && kNue2018CVNCutLowBin,
        kNue2018NDPresel && kNue2017MRParentSliceCut,
        kNue2018NDPresel && kNue2017MRParentSliceCut && kNue2018CVNCutHighBin,
    };
    
    const std::string selNames[kNumSels] = {
        "nd-presel-mre",
        "nd-full-mre",
        "nd-presel-mre2",
        "nd-full-mre-low",
        "nd-presel-mre3",
        "nd-full-mre-high"
    };
    
    if(!plot){

       std::string fname_data = "";
       std::string fname = "";
       std::string fname_calibxy_pos = "";
       std::string fname_calibxy_neg = "";
       std::string fname_calib = "";
       std::string fname_light_up = "";
       std::string fname_light_down = "";
       std::string fname_ckv = "";
       
       if(isFHC){
         fname_data = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_numi_fhc_full_v1_goodruns";
         fname = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_fhc_full_v1 with stride 1";
         fname_calibxy_pos = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1 with stride 1";
         fname_calibxy_neg = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1 with stride 1";
         fname_calib = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1 with stride 1";
         fname_light_up = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1 with stride 1";
         fname_light_down = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1 with stride 1";
         fname_ckv = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1 with stride 1";
       }
       if(!isFHC) {
         fname_data = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_numi_rhc_full_v1_goodruns";
         fname = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_rhc_full_v1 with stride 1";
         fname_calibxy_pos = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_rhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1 with stride 1";
         fname_calibxy_neg = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_rhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1 with stride 1";
         fname_calib = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_rhc_nova_v08_full_calib-shift-nd-func_v1 with stride 1";
         fname_light_up = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_rhc_nova_v08_full_lightmodel-lightup-calibdown_v1 with stride 1";
         fname_light_down = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_rhc_nova_v08_full_lightmodel-lightdown-calibup_v1 with stride 1";
         fname_ckv = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_nonswap_rhc_nova_v08_full_ckv-proton-shift-down_v1 with stride 1";
       }

       SpectrumLoader loader_data(fname_data);
       SpectrumLoader loader(fname);
       SpectrumLoader loader_calibxy_pos(fname_calibxy_pos);
       SpectrumLoader loader_calibxy_neg(fname_calibxy_neg);
       SpectrumLoader loader_calib(fname_calib);
       SpectrumLoader loader_light_up(fname_light_up);
       SpectrumLoader loader_light_down(fname_light_down);
       SpectrumLoader loader_ckv(fname_ckv);

       loader_data.SetSpillCut(kStandardSpillCuts);
       loader.SetSpillCut(kStandardSpillCuts);
       loader_calibxy_pos.SetSpillCut(kStandardSpillCuts);
       loader_calibxy_neg.SetSpillCut(kStandardSpillCuts);
       loader_calib.SetSpillCut(kStandardSpillCuts);
       loader_light_up.SetSpillCut(kStandardSpillCuts);
       loader_light_down.SetSpillCut(kStandardSpillCuts);
       loader_ckv.SetSpillCut(kStandardSpillCuts);

       Spectrum* spec_data[kNumSels][kNumVars];
       Spectrum* spec[kNumSels][kNumVars];
       Spectrum* spec_calibxy_pos[kNumSels][kNumVars];
       Spectrum* spec_calibxy_neg[kNumSels][kNumVars];
       Spectrum* spec_calib[kNumSels][kNumVars];
       Spectrum* spec_light_up[kNumSels][kNumVars];
       Spectrum* spec_light_down[kNumSels][kNumVars];
       Spectrum* spec_ckv[kNumSels][kNumVars];

       for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
        for (int varIdx=0; varIdx<kNumVars;varIdx++){
            const HistAxis& axis = defs[varIdx].axis;
            spec_data[selIdx][varIdx] = new Spectrum(loader_data, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
            spec[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
            spec_calibxy_pos[selIdx][varIdx] = new Spectrum(loader_calibxy_pos, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
            spec_calibxy_neg[selIdx][varIdx] = new Spectrum(loader_calibxy_neg, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
            spec_calib[selIdx][varIdx] = new Spectrum(loader_calib, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
            spec_light_up[selIdx][varIdx] = new Spectrum(loader_light_up, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
            spec_light_down[selIdx][varIdx] = new Spectrum(loader_light_down, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
            spec_ckv[selIdx][varIdx] = new Spectrum(loader_ckv, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
          }
      }

       loader_data.Go();
       loader.Go();
       loader_calibxy_pos.Go();
       loader_calibxy_neg.Go();
       loader_calib.Go();
       loader_light_up.Go();
       loader_light_down.Go();
       loader_ckv.Go();

       TFile* file = new TFile(file_mc.c_str(),"recreate");

       for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
          TDirectory* d = file->mkdir(selNames[selIdx].c_str());
          for (int varIdx=0; varIdx<kNumVars;varIdx++){
            const char* name = defs[varIdx].name.c_str();
            spec_data[selIdx][varIdx]->SaveTo(d, TString::Format("data_%s", name));
            spec[selIdx][varIdx]->SaveTo(d, TString::Format("%s", name));
            spec_calibxy_pos[selIdx][varIdx]->SaveTo(d, TString::Format("calibxy_pos_%s", name));
            spec_calibxy_neg[selIdx][varIdx]->SaveTo(d, TString::Format("calibxy_neg_%s", name));
            spec_calib[selIdx][varIdx]->SaveTo(d, TString::Format("calib_%s", name));
            spec_light_up[selIdx][varIdx]->SaveTo(d, TString::Format("light_up_%s", name));
            spec_light_down[selIdx][varIdx]->SaveTo(d, TString::Format("light_down_%s", name));
            spec_ckv[selIdx][varIdx]->SaveTo(d, TString::Format("ckv_%s", name));
          }
       }
       file->Close();
    }

    if(plot){
   
    Spectrum *data[kNumSels][kNumVars];
    Spectrum *mc[kNumSels][kNumVars];
    Spectrum *mc_calibxy_pos[kNumSels][kNumVars];
    Spectrum *mc_calibxy_neg[kNumSels][kNumVars];
    Spectrum *mc_calib[kNumSels][kNumVars];
    Spectrum *mc_light_up[kNumSels][kNumVars];
    Spectrum *mc_light_down[kNumSels][kNumVars];
    Spectrum *mc_ckv[kNumSels][kNumVars];
 
    for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
        const char* selName = selNames[selIdx].c_str();
        for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
            const char* varName = defs[varIdx].name.c_str();
            data[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/data_%s", selName, varName).Data()).release();
            mc[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/%s", selName, varName).Data()).release();  
            mc_calibxy_pos[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/calibxy_pos_%s", selName, varName).Data()).release();
            mc_calibxy_neg[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/calibxy_neg_%s", selName, varName).Data()).release();
            mc_calib[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/calib_%s", selName, varName).Data()).release();
            mc_light_up[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/light_up_%s", selName, varName).Data()).release();
            mc_light_down[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/light_down_%s", selName, varName).Data()).release();
            mc_ckv[selIdx][varIdx] = LoadFromFile<Spectrum>(file_mc, TString::Format("%s/ckv_%s", selName, varName).Data()).release();
            
        } // end for varIdx
    } // end for selIdx

    TH1* hData[kNumSels][kNumVars];
    TH1* hMC[kNumSels][kNumVars];
    TH1* hMC_calibxy_pos[kNumSels][kNumVars];
    TH1* hMC_calibxy_neg[kNumSels][kNumVars];
    TH1* hMC_calib[kNumSels][kNumVars];
    TH1* hMC_light_up[kNumSels][kNumVars];
    TH1* hMC_light_down[kNumSels][kNumVars];
    TH1* hMC_ckv[kNumSels][kNumVars];

    double POTmc = data[0][0]->POT();
    cout<<POTmc<<endl;

    for(int selIdx = 0; selIdx < 2; ++selIdx){
        const char* selName = selNames[selIdx].c_str();
        for(int varIdx = 0; varIdx < 2; ++varIdx){
            const char* varName = defs[varIdx].name.c_str();
            hData[selIdx][varIdx]=data[selIdx][varIdx]->ToTH1(POTmc, kBlack, kSolid, ana::EExposureType::kPOT);
            hMC[selIdx][varIdx]=mc[selIdx][varIdx]->ToTH1(POTmc, kBlack, kSolid, ana::EExposureType::kPOT);
            hMC_calibxy_pos[selIdx][varIdx]=mc_calibxy_pos[selIdx][varIdx]->ToTH1(POTmc, kBlue, kSolid, ana::EExposureType::kPOT);
            hMC_calibxy_neg[selIdx][varIdx]=mc_calibxy_neg[selIdx][varIdx]->ToTH1(POTmc, kRed, kSolid, ana::EExposureType::kPOT);
            hMC_calib[selIdx][varIdx]=mc_calib[selIdx][varIdx]->ToTH1(POTmc, kGreen+2, kSolid, ana::EExposureType::kPOT);
            hMC_light_up[selIdx][varIdx]=mc_light_up[selIdx][varIdx]->ToTH1(POTmc, kBlue, kSolid, ana::EExposureType::kPOT);
            hMC_light_down[selIdx][varIdx]=mc_light_down[selIdx][varIdx]->ToTH1(POTmc, kRed, kSolid, ana::EExposureType::kPOT);
            hMC_ckv[selIdx][varIdx]=mc_ckv[selIdx][varIdx]->ToTH1(POTmc, kGreen+2, kSolid, ana::EExposureType::kPOT);
             

            for(auto &h:{hData[selIdx][varIdx], hMC[selIdx][varIdx], hMC_calibxy_pos[selIdx][varIdx], hMC_calibxy_neg[selIdx][varIdx], hMC_calib[selIdx][varIdx], hMC_light_up[selIdx][varIdx], hMC_light_down[selIdx][varIdx], hMC_ckv[selIdx][varIdx]}){
                if(defs[varIdx].rebin>1) h->Rebin(defs[varIdx].rebin);
                h->GetXaxis()->SetRange(defs[varIdx].lowedge, defs[varIdx].highedge/defs[varIdx].rebin);
            }

            double nom_mean = hMC[selIdx][varIdx]->GetMean();
            double pos_mean = hMC_calibxy_pos[selIdx][varIdx]->GetMean();
            double neg_mean = hMC_calibxy_neg[selIdx][varIdx]->GetMean();
            double shape_mean = hMC_calib[selIdx][varIdx]->GetMean();
            string nom_name = "Nominal (mean "+std::to_string(nom_mean)+")";
            string pos_name = "Calib XY pos (mean "+std::to_string(pos_mean)+")";
            string neg_name = "Calib XY neg (mean "+std::to_string(neg_mean)+")";
            string shape_name = "Calib shape (mean "+std::to_string(shape_mean)+")";

            auto c = MakeMePlot({hMC[selIdx][varIdx], hMC_calibxy_pos[selIdx][varIdx], hMC_calibxy_neg[selIdx][varIdx], hMC_calib[selIdx][varIdx]}, {nom_name, pos_name, neg_name, shape_name}, isFHC);
            //c->SaveAs(("blpl/"+mode_tag+"_plot-calib-"+selName+"-"+varName+".pdf").c_str());
            delete c;

            auto c1 = MakeMePlot({hMC[selIdx][varIdx], hMC_light_up[selIdx][varIdx], hMC_light_down[selIdx][varIdx], hMC_ckv[selIdx][varIdx]}, {"Nominal", "Light up", "Light down", "Cherenkov"}, isFHC);
            //c1->SaveAs(("blpl/"+mode_tag+"_plot-light-ckv-"+selName+"-"+varName+".pdf").c_str());
            delete c1;

            if(selIdx>0 &&  (selIdx%2!=0)){

                TH1* hEff_mc = (TH1*)hMC[selIdx][varIdx]->Clone();

                TH1* hEff_mc_denom = (TH1*)hMC[selIdx-1][varIdx]->Clone();
                TH1* hEff_data = (TH1*)hData[selIdx][varIdx]->Clone();
                TH1* hEff_data_denom = (TH1*)hData[selIdx-1][varIdx]->Clone();
                TH1* hEff_mc_calibxy_pos = (TH1*)hMC_calibxy_pos[selIdx][varIdx]->Clone();
                TH1* hEff_mc_calibxy_pos_denom = (TH1*)hMC_calibxy_pos[selIdx-1][varIdx]->Clone();
                TH1* hEff_mc_calibxy_neg = (TH1*)hMC_calibxy_neg[selIdx][varIdx]->Clone();
                TH1* hEff_mc_calibxy_neg_denom = (TH1*)hMC_calibxy_neg[selIdx-1][varIdx]->Clone();
                TH1* hEff_mc_calib = (TH1*)hMC_calib[selIdx][varIdx]->Clone();
                TH1* hEff_mc_calib_denom = (TH1*)hMC_calib[selIdx-1][varIdx]->Clone();
                TH1* hEff_mc_light_up = (TH1*)hMC_light_up[selIdx][varIdx]->Clone();
                TH1* hEff_mc_light_up_denom = (TH1*)hMC_light_up[selIdx-1][varIdx]->Clone();
                TH1* hEff_mc_light_down = (TH1*)hMC_light_down[selIdx][varIdx]->Clone();
                TH1* hEff_mc_light_down_denom = (TH1*)hMC_light_down[selIdx-1][varIdx]->Clone();
                TH1* hEff_mc_ckv = (TH1*)hMC_ckv[selIdx][varIdx]->Clone();
                TH1* hEff_mc_ckv_denom = (TH1*)hMC_ckv[selIdx-1][varIdx]->Clone();

                hEff_mc->Divide(hEff_mc_denom);
                hEff_data->Divide(hEff_data_denom);
                hEff_mc_calibxy_pos->Divide(hEff_mc_calibxy_pos_denom);
                hEff_mc_calibxy_neg->Divide(hEff_mc_calibxy_neg_denom);
                hEff_mc_calib->Divide(hEff_mc_calib_denom);
                hEff_mc_light_up->Divide(hEff_mc_light_up_denom);
                hEff_mc_light_down->Divide(hEff_mc_light_down_denom);
                hEff_mc_ckv->Divide(hEff_mc_ckv_denom);

                auto c2 = MakeMePlot({hEff_mc, hEff_mc_calibxy_pos, hEff_mc_calibxy_neg, hEff_mc_calib}, {"Nominal", "Calib XY pos", "Calib XY neg", "Calib shape"}, isFHC, "Selection efficiency");
                //c2->SaveAs(("blpl/"+mode_tag+"_plot-efficiency-calib-"+selName+"-"+varName+".pdf").c_str());
                delete c2;
                 
                auto c3 = MakeMePlot({hEff_mc, hEff_mc_light_up, hEff_mc_light_down, hEff_mc_ckv}, {"Nominal", "Light up", "Light down", "Cherenkov"}, isFHC, "Selection efficiency");
                //c3->SaveAs(("blpl/"+mode_tag+"_plot-efficiency-light-"+selName+"-"+varName+".pdf").c_str());
                delete c3;
          
                hEff_mc_light_up->SetLineColor(kViolet);
                hEff_mc_light_down->SetLineColor(kPink);
                hEff_mc_calib->SetLineColor(kOrange);
                auto c4 = MakeMePlot({hEff_mc, hEff_mc_light_up, hEff_mc_light_down, hEff_mc_ckv, hEff_mc_calibxy_pos, hEff_mc_calibxy_neg, hEff_mc_calib}, {"Nominal", "Light up", "Light down", "Cherenkov", "Calib XY pos", "Calib XY neg", "Calib shape"}, isFHC, "Selection efficiency");
                //c4->SaveAs(("blpl/"+mode_tag+"_plot-efficiency-"+selName+"-"+varName+".pdf").c_str());
                delete c4;

               // blessing plot - Data vs MC eff band + syst. band               
    
                TCanvas *c5 = new TCanvas("c5", "canvas", 1000, 1000);
                TPad* p1 = new TPad("", "", 0, 0, 1, 1);
                TPad* p2 = new TPad("", "", 0, 0, 1, 1);
                p1->SetBottomMargin(.3);
                p2->SetTopMargin(.7);

                for(auto p:{p1,p2}){
                    p->SetLeftMargin(0.105);
                    p->SetFillStyle(0);
                    p->Draw();
                }
                p1->cd();

                /*Int_t ci = TColor::GetFreeColorIndex();
                TColor *color = new TColor(ci, 155/255.0, 233/255.0, 233/255.0);
                Int_t line = TColor::GetFreeColorIndex(); 
                TColor *color_line = new TColor(line, 0/255.0, 39/255.0, 34/255.0);
                Int_t line_data = TColor::GetFreeColorIndex();
                TColor *color_data = new TColor(line_data, 255/255.0, 13/255.0, 109/255.0);*/ // too cheerful
                Int_t ci = TColor::GetFreeColorIndex();
                TColor *color = new TColor(ci, 216/255.0, 216/255.0, 216/255.0); // lighter gray
                //ci = kGray;
                Int_t line = kGray+2;
                Int_t line_data = kBlack;
                if(drawing_options.Contains("red")) {
                   ci= kRed-10;
                   line = kPink;//kRed;
                   line_data = kBlack;}

                hEff_data->SetLineColor(line_data);
                hEff_mc->SetLineColor(line);
                hEff_data->SetLineWidth(2);
                hEff_mc->SetLineWidth(2);
                hEff_mc->SetMaximum(1.0);
                hEff_mc->SetMinimum(0);
                hEff_mc->GetXaxis()->SetRange(0, 45);
                hEff_mc->GetYaxis()->CenterTitle();
                hEff_mc->GetXaxis()->CenterTitle();
                hEff_mc->GetYaxis()->SetTitle("Selection efficiency");
                hEff_mc->Draw("hist E1");
                std::vector<TH1*> syst = {hEff_mc_light_up, hEff_mc_light_down, hEff_mc_ckv, hEff_mc_calibxy_pos, hEff_mc_calibxy_neg, hEff_mc_calib};
                TGraphAsymmErrors* g = new TGraphAsymmErrors;
                TH1F* hUp =  (TH1F*)hEff_mc->Clone();
                TH1F* hDown =  (TH1F*)hEff_mc->Clone();
                for(int binIdx = 0; binIdx < hEff_mc->GetNbinsX()+2; ++binIdx){
                    const double y = hEff_mc->GetBinContent(binIdx);
                    g->SetPoint(binIdx, hEff_mc->GetXaxis()->GetBinCenter(binIdx), y);
                    const double w = hEff_mc->GetXaxis()->GetBinWidth(binIdx);
                    double errUp = 0;
                    double errDn = 0;
                    for(unsigned int systIdx = 0; systIdx < syst.size(); ++systIdx){
                        double shift = syst[systIdx]->GetBinContent(binIdx)-y;
                        if(shift < 0) errDn += shift*shift;
                        if(shift >= 0) errUp+= shift*shift;
                    }
                    g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));       
                    hUp->SetBinContent(binIdx, hEff_mc->GetBinContent(binIdx)+sqrt(errUp));
                    hDown->SetBinContent(binIdx, hEff_mc->GetBinContent(binIdx) - sqrt(errDn));
                }
                g->SetFillColor(ci);
                g->Draw("e2 same");
                hEff_mc->SetLineWidth(3);
                hEff_mc->Draw("hist E1 same");
                hEff_data->SetMarkerStyle(8);
                hEff_data->SetMarkerSize(2);
                hEff_data->Draw("E1 same");
                c5->RedrawAxis();
                gPad->Update();
                TText *tt; 
                double systname_x = 0.02;
                double leg_x = 0.63;
                double leg_x_max = 0.85; 
                if (varIdx == 0){
                   systname_x = 0.55;
                   leg_x = 0.15;
                   leg_x_max = 0.35;
                   if (!isFHC){
                      leg_x = 0.12;
                      leg_x_max = 0.3;
                   }
                }
                if(isFHC){
                  auto *tt = new TText(gPad->GetUxmax()*0.55, 1.05*gPad->GetUymax(), "Neutrino beam, MRE");
                  tt->SetTextSize(0.04);
                  tt->Draw();
                  tt = new TText(gPad->GetUxmax()*systname_x, 0.9*gPad->GetUymax(), "Calibration + lightlevel syst.");
                  tt->SetTextSize(0.03);
                  tt->Draw();}
                else{
                  auto *tt = new TText(gPad->GetUxmax()*0.5, 1.05*gPad->GetUymax(), "Antineutrino beam, MRE");
                  tt->SetTextSize(0.04);
                  tt->Draw();
                  tt = new TText(gPad->GetUxmax()*systname_x, 0.9*gPad->GetUymax(), "Calibration + lightlevel syst.");
                  tt->SetTextSize(0.03);
                  tt->Draw();}
 
                 TLegend *leg = new TLegend(leg_x,0.73,leg_x_max,0.85);
                 leg->SetTextSize(0.033);
                 leg->AddEntry(hEff_data, "Data", "p");
                 leg->AddEntry(hEff_mc, "Monte Carlo", "l");
                 leg->AddEntry(g, "1#sigma syst. error", "f");
                 leg->Draw();

                 Prelim();

                 gPad->RedrawAxis();

                 p2->cd();
                 auto temp = (TH1*)hEff_mc->Clone();
                 temp->Divide(hEff_data);
                 temp->Draw();
                 auto temp_up = (TH1*)hUp->Clone();
                 auto temp_down = (TH1*)hDown->Clone();
                 temp_up->Divide(hEff_data);
                 temp_down->Divide(hEff_data);

                 TGraphAsymmErrors* g_ratio = new TGraphAsymmErrors;
                 for(int binIdx = 0; binIdx < temp->GetNbinsX()+2; ++binIdx){
                     const double w = temp->GetXaxis()->GetBinWidth(binIdx);
                     g_ratio->SetPoint(binIdx, temp->GetXaxis()->GetBinCenter(binIdx), temp->GetBinContent(binIdx));
                     g_ratio->SetPointError(binIdx, w/2, w/2, temp->GetBinContent(binIdx) - temp_down->GetBinContent(binIdx), temp_up->GetBinContent(binIdx) - temp->GetBinContent(binIdx));      
                 }

                 g_ratio->SetFillColor(ci); 
                 g_ratio->Draw("e2 same");
                 temp->Draw("same");
                 TLine *l = new TLine(-0.1,1,4.5,1);
                 l->SetLineColor(kBlack);
                 l->Draw();
                 temp->GetYaxis()->SetRangeUser(0.5, 1.5);
                 temp->GetYaxis()->SetTitle("MC / Data");
                 p2->RedrawAxis();

                 c5->SaveAs(("blpl/"+mode_tag+"_plot-efficiency-error-band-"+selName+"-"+varName+"_"+drawing_options+".pdf"));
                } // effic
         } // var
    } // sel

     cout<<"Table for "<<mode_tag<<endl;
     cout<<"Data "<<hData[0][0]->Integral()<<" "<<hData[1][0]->Integral()<<" "<<hData[1][0]->Integral()/hData[0][0]->Integral()<<endl;
     cout<<"Data "<<hMC[0][0]->Integral()<<" "<<hMC[1][0]->Integral()<<" "<<hMC[1][0]->Integral()/hMC[0][0]->Integral()<<endl;
     cout<<"difference "<<(hMC[1][0]->Integral()/hMC[0][0]->Integral() - hData[1][0]->Integral()/hData[0][0]->Integral())/(hMC[1][0]->Integral()/hMC[0][0]->Integral())<<endl;
    
    // bl plot with Data vs MC at presel and full sel.

    for(int varIdx = 0; varIdx < 2; ++varIdx){
        const char* varName = defs[varIdx].name.c_str();   
        TCanvas *c6 = new TCanvas("c6", "canvas", 1000, 800);
        c6->SetLeftMargin(0.105);

        Int_t ci = TColor::GetFreeColorIndex();
        TColor *color = new TColor(ci, 216/255.0, 216/255.0, 216/255.0); // lighter gray
        //ci = kGray;
        Int_t line = kGray+2;
        Int_t line_data = kBlack;
        if(drawing_options.Contains("red")) {
           ci= kRed-10;
           line = kPink;//kRed;
           line_data = kBlack;
        }
        
        auto hMC_presel = hMC[0][varIdx];
        auto hMC_full = hMC[1][varIdx];
        auto hData_presel = hData[0][varIdx];
        auto hData_full = hData[1][varIdx];

        hMC_presel->SetLineStyle(7);
        hData_presel->SetLineStyle(7); 
        hData_presel->SetLineColor(line);
        hData_full->SetLineColor(line_data);
        hMC_presel->SetLineColor(line);
        hMC_full->SetLineColor(line_data);
        for (auto &h:{hMC_presel, hMC_full, hData_presel, hData_full}){
             h->SetLineWidth(2);
             h->Scale(1/10000.0);
        }
        hData_presel->SetMarkerStyle(8);
        hData_presel->SetMarkerColor(line);
        hData_presel->SetMarkerSize(2);
        hData_full->SetMarkerStyle(8);
        hData_full->SetMarkerSize(2);
        hData_presel->SetMaximum(hData_presel->GetMaximum()*1.2);
        hData_presel->SetMinimum(0);
        hData_presel->GetXaxis()->SetRange(0, 45);
        hData_presel->GetYaxis()->CenterTitle();
        hData_presel->GetXaxis()->CenterTitle();
        hData_presel->GetYaxis()->SetTitle(TString::Format("10^{4} events / %1.2g #times 10^{20} POT", POTmc/1E20).Data());
        hData_presel->Draw("E1");
        hMC_presel->Draw("hist E1 same");
        hData_full->Draw("E1 same");
        hMC_full->Draw("hist E1 same");
        c6->RedrawAxis();
        gPad->Update();
        TText *tt;
        if(isFHC){
                auto *tt = new TText(gPad->GetUxmax()*0.6, 1.05*gPad->GetUymax(), "Neutrino beam, MRE");
                tt->SetTextSize(0.04);
                tt->Draw();
        }
        else{
                 auto *tt = new TText(gPad->GetUxmax()*0.55, 1.05*gPad->GetUymax(), "Antineutrino beam, MRE");
                 tt->SetTextSize(0.04);
                 tt->Draw();
        }

        TLegend *leg = new TLegend(0.5, 0.65, 0.8, 0.85);
        leg->SetTextSize(0.033);
        leg->AddEntry(hData_presel, "Data, preselected", "p");
        leg->AddEntry(hMC_presel, "Monte Carlo, preselected", "l");
        leg->AddEntry(hData_full, "Data, fully selected", "p");
        leg->AddEntry(hMC_full, "Monte Carlo, fully selected", "l");
        leg->Draw();

        Prelim();

        gPad->RedrawAxis();

        c6->SaveAs(("blpl/"+mode_tag+"_plot-data-mc-presel-and-fullsel-"+varName+"_"+drawing_options+".pdf"));

     }

    }// is plot    
}