mre_comp_split.C

Go to the documentation of this file.

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/LoadFromFile.h"

#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"

#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "3FlavorAna/Vars/NumuVars.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TROOT.h"
#include "TRint.h"
#include "TColor.h"
#include "TFile.h"
#include "TLine.h"
#include "TText.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TPad.h"

#include <fstream>
#include <iomanip>
#include <iostream>

#include "Utilities/func/MathUtil.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();
}

using namespace ana;

void mre_comp_split(bool plot = true, bool isFHC = true){

  const std::string mode_tag = isFHC?"fhc":"rhc";
  std::string file_name = "mre_pred_no_extrap_nd_cuts_"+mode_tag+"_mc_vs_data.root";

  const int kNumVars = 3;

  struct HistDef
  {
        std::string name;
        HistAxis axis;
  };
  
  const Var kSlcE([](const caf::SRProxy *sr){
        return sr->slc.calE;
  });

  const HistDef defs[kNumVars] = {
        {"energy", {"#nu_{e} Energy(GeV)", Binning::Simple(120,0,12), kNueEnergy2018}},
        {"caloE", {"Calorimetric energy", Binning::Simple(120,0,12), kSlcE}},
        {"cvne", {"CVNe", Binning::Simple(120,-0.1,1.1), kCVNSSe}},
  };

  const int kNumSels=2;

  const Cut sels[kNumSels] = {
        kNue2018NDPresel && kNue2017MRParentSliceCut,
        kNue2018NDCVNSsb && kNue2017MRParentSliceCut
  };

  const std::string selNames[kNumSels] = {
        "nd-presel-mre",
        "nd-full-mre"
  };
    
  if (!plot){

    std::string fname = "";
    std::string fname_data = "";

    // FHC MC
   if(isFHC){
        fname = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_fhc_full_v1 with stride 1";
        fname_data = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_numi_fhc_full_v1_goodruns";
     }
    // RHC MC
    if(!isFHC) {
       fname = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_genie_rhc_full_v1 with stride 1";
       fname_data = "defname: prod_mrecaf_R17-11-14-prod4mre.b_nd_numi_rhc_full_v1_goodruns";
     }

    
    SpectrumLoader loader(fname);
    SpectrumLoader loader_data(fname_data);

    loader.SetSpillCut(kStandardSpillCuts);
    loader_data.SetSpillCut(kStandardSpillCuts);
 
    
    Spectrum* spec_data[kNumSels][kNumVars];
    Spectrum* spec[kNumSels][kNumVars];
    Spectrum* spec_dis[kNumSels][kNumVars];
    Spectrum* spec_res[kNumSels][kNumVars];
    Spectrum* spec_qe[kNumSels][kNumVars];
    Spectrum* spec_mec[kNumSels][kNumVars];
    Spectrum* spec_coh[kNumSels][kNumVars];
    
    for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
        for (int varIdx=0; varIdx<kNumVars;varIdx++){
            const HistAxis& axis = defs[varIdx].axis;
            spec[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx], kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // All int
            spec_dis[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsDIS, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // DIS
            spec_res[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsRes, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // RES
            spec_qe[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsQE, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // QE
            spec_mec[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsDytmanMEC, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // MEC
            spec_coh[selIdx][varIdx] = new Spectrum(loader, axis, sels[selIdx] && kIsCoh, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018); // COH
            spec_data[selIdx][varIdx] = new Spectrum(loader_data, axis, sels[selIdx]);
        }
    }
    
    loader.Go();
    loader_data.Go();
    
    TFile* file = new TFile(file_name.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[selIdx][varIdx]->SaveTo(d, TString::Format("%s", name));
            spec_dis[selIdx][varIdx]->SaveTo(d, TString::Format("dis_%s", name));
            spec_res[selIdx][varIdx]->SaveTo(d, TString::Format("res_%s", name));
            spec_qe[selIdx][varIdx]->SaveTo(d, TString::Format("qe_%s", name));
            spec_mec[selIdx][varIdx]->SaveTo(d, TString::Format("mec_%s", name));
            spec_coh[selIdx][varIdx]->SaveTo(d, TString::Format("coh_%s", name));
            spec_data[selIdx][varIdx] ->SaveTo(d, TString::Format("data_%s", name));
        }
    }
    file->Close();
  } // end of !plot    

  if (plot){
     Spectrum* spec_data[kNumSels][kNumVars];
     Spectrum* spec[kNumSels][kNumVars];
     Spectrum* spec_dis[kNumSels][kNumVars];
     Spectrum* spec_res[kNumSels][kNumVars];
     Spectrum* spec_qe[kNumSels][kNumVars];
     Spectrum* spec_mec[kNumSels][kNumVars];
     Spectrum* spec_coh[kNumSels][kNumVars];

     TH1* hMC[kNumSels][kNumVars];
     TH1* hMC_dis[kNumSels][kNumVars];
     TH1* hMC_res[kNumSels][kNumVars];
     TH1* hMC_qe[kNumSels][kNumVars];
     TH1* hMC_mec[kNumSels][kNumVars];
     TH1* hMC_coh[kNumSels][kNumVars];
     TH1* hData[kNumSels][kNumVars];

     double POT = -1;
  
     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();

            spec[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/%s", selName, varName).Data()).release();
            spec_dis[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/dis_%s", selName, varName).Data()).release();
            spec_res[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/res_%s", selName, varName).Data()).release();
            spec_qe[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/qe_%s", selName, varName).Data()).release();
            spec_mec[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/mec_%s", selName, varName).Data()).release();
            spec_coh[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/coh_%s", selName, varName).Data()).release();
            spec_data[selIdx][varIdx] = LoadFromFile<Spectrum>(file_name, TString::Format("%s/data_%s", selName, varName).Data()).release();

            POT = spec_data[selIdx][varIdx]->POT();
            cout<<"data pot is "<<POT<<endl;

            hMC[selIdx][varIdx]=spec[selIdx][varIdx]->ToTH1(POT, kPink, kSolid, ana::EExposureType::kPOT);
            hMC_dis[selIdx][varIdx]=spec_dis[selIdx][varIdx]->ToTH1(POT, kGray+1, kSolid, ana::EExposureType::kPOT);
            hMC_res[selIdx][varIdx]=spec_res[selIdx][varIdx]->ToTH1(POT, kGreen+2, kSolid, ana::EExposureType::kPOT);
            hMC_qe[selIdx][varIdx]=spec_qe[selIdx][varIdx]->ToTH1(POT, kAzure+7, kSolid, ana::EExposureType::kPOT);
            hMC_mec[selIdx][varIdx]=spec_mec[selIdx][varIdx]->ToTH1(POT, kOrange-3, kSolid, ana::EExposureType::kPOT);
            hMC_coh[selIdx][varIdx]=spec_coh[selIdx][varIdx]->ToTH1(POT, kViolet-1, kSolid, ana::EExposureType::kPOT);
            hData[selIdx][varIdx]=spec_data[selIdx][varIdx]->ToTH1(POT, kBlack, kSolid, ana::EExposureType::kPOT);

        } // end for varIdx
    } // end for selIdx  


    // now let's make one plot for blessing - CVNe for presel with different int. types
    int selIdx = 0;
    int varIdx = 2;

    cout<<"integrals "<<hMC[selIdx][varIdx]->Integral()<<" "<<hMC_dis[selIdx][varIdx]->Integral()<<" "<<hMC_res[selIdx][varIdx]->Integral()<<" "<<hMC_qe[selIdx][varIdx]->Integral()<<" "<<hMC_mec[selIdx][varIdx]->Integral()<<" "<<hMC_coh[selIdx][varIdx]->Integral()<<" "<<hData[selIdx][varIdx]->Integral()<<endl;

     TCanvas *c = new TCanvas("c", "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();
    
     for (auto &h:{hMC[selIdx][varIdx], hMC_dis[selIdx][varIdx], hMC_res[selIdx][varIdx], hMC_qe[selIdx][varIdx], hMC_mec[selIdx][varIdx], hMC_coh[selIdx][varIdx], hData[selIdx][varIdx]}){
             //h->SetLineWidth(2);
             h->Rebin(2);
             h->Scale(1/10000.0);
     } 
     hMC[selIdx][varIdx]->SetMaximum(hMC[selIdx][varIdx]->GetMaximum()*1.2);
     hMC[selIdx][varIdx]->SetMinimum(0);
     hMC[selIdx][varIdx]->GetYaxis()->CenterTitle();
     hMC[selIdx][varIdx]->GetXaxis()->CenterTitle();
     hMC[selIdx][varIdx]->GetYaxis()->SetTitleSize(0.045);
     hMC[selIdx][varIdx]->GetYaxis()->SetTitleOffset(1.17);
     hMC[selIdx][varIdx]->GetYaxis()->SetTitle(TString::Format("10^{4} events / %1.2g #times 10^{20} POT", POT/1E20).Data());

     hMC[selIdx][varIdx]->Draw("hist");
     hMC_dis[selIdx][varIdx]->Draw("hist same");
     hMC_res[selIdx][varIdx]->Draw("hist same");
     hMC_qe[selIdx][varIdx]->Draw("hist same");
     hMC_mec[selIdx][varIdx]->Draw("hist same");
     hMC_coh[selIdx][varIdx]->Draw("hist same");
     hData[selIdx][varIdx]->SetMarkerStyle(20);
     hData[selIdx][varIdx]->Draw("E1 same");

     c->RedrawAxis();
     gPad->Update();
   
     TText *tt;
     TText *tt2;
     TLine *l2;
     if(isFHC){
        auto *tt = new TText(gPad->GetUxmax()*0.53, 1.05*gPad->GetUymax(), "Neutrino beam, MRE");
        tt->SetTextSize(0.04);
        tt->Draw();
        auto *l2 = new TLine(0.84,0,0.84,gPad->GetUymax());
        l2->SetLineColor(kBlack);
        l2->Draw();
        auto *tt2 = new TText(0.82, 0.6*gPad->GetUymax(), "Cut position");
        tt2->SetTextAngle(90.0);
        tt2->SetTextSize(0.035);
        tt2->Draw();
     }
     else{
        auto *tt = new TText(gPad->GetUxmax()*0.46, 1.05*gPad->GetUymax(), "Antineutrino beam, MRE");
        tt->SetTextSize(0.04);
        tt->Draw();
        auto *l2 = new TLine(0.89,0,0.89,gPad->GetUymax());
        l2->SetLineColor(kBlack);
        l2->Draw();
        auto *tt2 = new TText(0.87, 0.6*gPad->GetUymax(), "Cut position");
        tt2->SetTextAngle(90.0);
        tt2->SetTextSize(0.035);
        tt2->Draw();
     }

     TLegend *leg = new TLegend(0.2,0.6,0.5,0.85);
     leg->SetTextSize(0.03);
     leg->AddEntry(hData[selIdx][varIdx], "Data", "p");
     leg->AddEntry(hMC[selIdx][varIdx], "Monte Carlo, Total", "l");
     leg->AddEntry(hMC_qe[selIdx][varIdx], "Monte Carlo, QE", "l");
     leg->AddEntry(hMC_res[selIdx][varIdx], "Monte Carlo, RES", "l");
     leg->AddEntry(hMC_mec[selIdx][varIdx], "Monte Carlo, MEC", "l");
     leg->AddEntry(hMC_dis[selIdx][varIdx], "Monte Carlo, DIS", "l");
     leg->AddEntry(hMC_coh[selIdx][varIdx], "Monte Carlo, COH", "l");
     leg->Draw();

     Prelim();
   
     p2->cd();
     auto temp = (TH1*)hData[selIdx][varIdx]->Clone();
     temp->Divide(hMC[selIdx][varIdx]);
     temp->SetLineColor(kPink);
     temp->SetMarkerStyle(1);
     temp->GetXaxis()->SetTitle("#nu_{e} CVN classifier");
     temp->Draw();
 
     TLine *l = new TLine(-0.1,1,1.1,1);
     l->SetLineColor(kBlack);
     l->Draw();
     temp->GetYaxis()->SetRangeUser(0.75, 1.25);
     temp->GetYaxis()->CenterTitle();
     temp->GetXaxis()->CenterTitle();
     temp->GetYaxis()->SetTitle("Data / MC");
     temp->GetYaxis()->SetTitleSize(0.04);
     temp->GetYaxis()->SetTitleOffset(1.2);
     p2->RedrawAxis();

     c->SaveAs(("blpl/"+mode_tag+"_plot-data-mc_cvne_presel.pdf").c_str());
  }    
}