plot_ND_DataMC_energybinning.C

Go to the documentation of this file.

#include "CAFAna/Analysis/DataMCPair.h"
#include "3FlavorAna/Plotting/SharedPlotLists.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "Utilities/rootlogon.C"

#include "TCanvas.h"
#include "TGraph.h"
#include "TGaxis.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TFile.h"
#include "TSystem.h"


#include <iostream>
#include <string>
#include <sstream>
#include <fstream>
#include <iomanip>

using namespace ana;

TString sNumu_Nue = "nue";
std::vector<TString> hcNames = {"fhc", "rhc"}; const unsigned int nHC = hcNames.size();
TString inDir = "/nova/ana/users/essmith/ND_DataMC2020/";
//TString inDir = "/nova/app/users/essmith/plotlists_dev/3FlavorAna/Ana2020/ND_DataMC/";
std::string outDir = "plots_energybinning/";

void SetAxis(TH1* h, bool visible, int color, int style = 1, bool scale = true);
TString POTstr(double pot);
double DerivePowScale(TH1* h);
TH1* MakeRatio(TH1* num, TH1* denom, int color=1);
void OneLine();
void Zero();
void Legend(TH1* h_to_clone);
void RatioLegend();
void NeutrinoMode(TString mode);

void plot_ND_DataMC_energybinning() {

  std::vector<HistDef> defs;
  AddHistDefRecoND(defs);
  AddHistDefSlice(defs);
  AddHistDefNueND(defs);
  AddHistDefNueShower(defs);
  AddHistDefNuePID(defs);
  AddHistDefNueSelectionExtras(defs);
  AddHistDefProngCVN(defs);
  AddHistDefNueEnergy(defs);
  unsigned int nVars = defs.size(); 
  std::vector<std::string> varNames;
  for(unsigned int var = 0; var<nVars; ++var){
    varNames.emplace_back(defs[var].name);
  }

  std::vector<std::string> selNames = {"Presel", "CVNeLow", "CVNeHigh"};
  unsigned int nSels = selNames.size();
  std::vector<std::string> compNames = {"nue", "numu", "nc", "antinue", "antinumu"};
  unsigned int nComps = compNames.size();

  std::unique_ptr<Spectrum> specData;
  double datapot[nHC][nVars][nSels];
  TH1* hData[nHC][nVars][nSels];
  TH1* hMC[nHC][nVars][nSels][nComps];
  TH1* hDecomp[nHC][nVars][nSels][nComps];

  unique_ptr<Spectrum> tempspec;

  for(unsigned int HC = 0; HC<nHC; ++HC){
    std::cout<<"Loading files for "<<hcNames[HC]<<std::endl;
    TString nueDataName = inDir+hcNames[HC]+"_nue_data.root";
    TString nueMCName = inDir+hcNames[HC]+"_nue_mc.root";
    TString nueMCDecompName = inDir+hcNames[HC]+"_nue_mc_decomp.root";

    TFile* nueDataFile = TFile::Open(nueDataName);
    TFile* nueMCFile = TFile::Open(nueMCName);
    TFile* nueMCDecompFile = TFile::Open(nueMCDecompName);

    for(unsigned int var = 0; var<nVars; ++var){
      for(unsigned int sel = 0; sel<nSels; ++sel){
        TString name = varNames[var]+"_"+selNames[sel];
        if(!name.Contains("CVN2d_CVNeLow")) continue;

        specData = Spectrum::LoadFrom(nueDataFile, name);
        datapot[nHC][nVars][nSels] = specData->POT();
        hData[HC][var][sel] = specData->ToTH1(datapot[nHC][nVars][nSels]);
        for(unsigned int comp = 0; comp<nComps; ++comp){
          tempspec = Spectrum::LoadFrom(nueMCFile, name+"_"+compNames[comp]);
          hMC[HC][var][sel][comp] = tempspec->ToTH1(datapot[nHC][nVars][nSels]);

          tempspec = Spectrum::LoadFrom(nueMCDecompFile, name+"_"+compNames[comp]);
          hDecomp[HC][var][sel][comp] = tempspec->ToTH1(datapot[nHC][nVars][nSels]);
        }
      }
    }

  // all up the components to get total MC
    for(unsigned int var = 0; var<nVars; ++var){
      for(unsigned int sel = 0; sel<nSels; ++sel){
        TString name = varNames[var]+"_"+selNames[sel];
        if(!name.Contains("CVN2d_CVNeLow")) continue;
        TString forfile = "_"+sNumu_Nue+"_"+hcNames[HC];

        TH1* totalMC = (TH1*)hMC[HC][var][sel][0]->Clone("totalMC");
        TH1* totalDecomp = (TH1*)hDecomp[HC][var][sel][0]->Clone("totalDecomp");

        for(unsigned int comp = 1; comp<nComps; ++comp){
          totalMC->Add(hMC[HC][var][sel][comp]);
          totalDecomp->Add(hDecomp[HC][var][sel][comp]);
        }

        TH1* rMC;
        TH1* rDecomp;
        // histograms for ratios
        rMC = MakeRatio(hData[HC][var][sel], totalMC, 17);
        rDecomp = MakeRatio(hData[HC][var][sel], totalDecomp);

        // Lets start drawing
        // Split plot
        TCanvas* c = new TCanvas(name+forfile, name+forfile);
        double ysplit=0.35;
        TPad* p1 =new TPad("p1", "", 0,0,1,1);
        TPad* p2 =new TPad("p2","",0,0,1,1);
        p1->SetBottomMargin(ysplit);
        p2->SetTopMargin(1-ysplit);
        p1->SetFillStyle(0);
        p2->SetFillStyle(0);

        p1->cd();


        std::vector <int> colors = {kMagenta,
                                    kGreen+2,
                                    kAzure,
                                    kMagenta-2,
                                    kGreen-3,
                                    kRed};
        std::vector<std::string> labels = {"#nu_{e} CC",
                                           "#nu_{#mu} CC",
                                           "NC",
                                           "#bar{#nu}_{e} CC",
                                           "#bar{#nu}_{#mu} CC"};

        totalDecomp->GetYaxis()->SetTitle(POTstr(datapot[nHC][nVars][nSels]));
        totalDecomp->GetYaxis()->CenterTitle();

        SetAxis(totalDecomp, true, colors[nComps]);
        double maxY = totalDecomp->GetMaximum();
        totalDecomp->GetYaxis()->SetRangeUser(0, 1.30*maxY);
        totalDecomp->Draw("hist");
        Legend(totalDecomp);

        Nue2018ThreeBinDivisions(true);
        Nue2018ThreeBinLabels(1.42*maxY, 0.06, kGray+3, true);        

                                 
        SetAxis(totalMC, false, colors[nComps], 7);
        totalMC->Draw("hist same");
        for(unsigned int comp = 0; comp<nComps; comp++){
          SetAxis(hDecomp[HC][var][sel][comp], false, colors[comp]);
          hDecomp[HC][var][sel][comp]->Draw("hist same");
          SetAxis(hMC[HC][var][sel][comp], false, colors[comp], 7);
          hMC[HC][var][sel][comp]->Draw("hist same");
        }
        hData[HC][var][sel]->SetMarkerStyle(kFullCircle); 
        SetAxis(hData[HC][var][sel], false, 1);
        hData[HC][var][sel]->Draw("ep same");

        Nue2018ThreeBinAxis(totalDecomp, false, true, true);
          
        totalDecomp->GetXaxis()->SetLabelSize(0.);
        totalDecomp->GetXaxis()->SetTitleSize(0.);
        totalDecomp->GetYaxis()->SetTitleSize(0.75*totalDecomp->GetYaxis()->GetTitleSize());
        totalDecomp->GetYaxis()->SetRangeUser(0.001,totalDecomp->GetMaximum());

        // The empty bins of the histograms still get drawn at y=0. Try to
        // draw over these colored lines with black
        Zero();

        // draw the ratio
        p2->cd();

        SetAxis(rDecomp, true, 1, 1, false);
        rDecomp->Draw("ep");

        SetAxis(rMC, false, 17, 1, false);
        rMC->Draw("ep same");
        OneLine();
        Nue2018ThreeBinDivisions(true);
        Nue2018ThreeBinAxis(rDecomp, true, true, true);

        rDecomp->GetYaxis()->SetTitleSize(0.75*rDecomp->GetYaxis()->GetTitleSize());
        rDecomp->GetXaxis()->SetLabelSize(0.75*rDecomp->GetYaxis()->GetTitleSize());
        rMC->GetYaxis()->SetRangeUser(.5,1.5);
        rDecomp->GetYaxis()->SetRangeUser(.5,1.5);
        RatioLegend();
        Preliminary();
        NeutrinoMode(hcNames[HC]);
        p1->Draw();
        p2->Draw("same");

        gPad->Print(outDir+name+forfile+".pdf");
      }
    }
    nueDataFile->Close();
    nueMCFile->Close();
    nueMCDecompFile->Close();

  }

  return;
}

void SetAxis(TH1* h, bool visible, int color, int style, bool scale) {
  if(scale) h->Scale(.001);
  h->SetLineStyle(style);
  h->SetLineColor(color);
  if(!visible){
    h->SetTickLength(0.);
    h->SetLabelSize(0.);
  }
  double norm = DerivePowScale(h);
}

TString POTstr(double pot) {
  return TString(TString::Format("10^{3} Events / %.3g#times10^{20} POT", pot/1e20)); 
}

double DerivePowScale(TH1* h) {
  // "Next smallest power of 1000 than maximum bin" - from plot_validation_datamc.py
  return 3*(int(log10(h->GetMaximum())))/3;
}

TH1* MakeRatio(TH1* num, TH1* denom, int color) {
  TH1* ret = (TH1*)num->Clone();
  ret->Divide(denom);

  ret->SetMarkerStyle(kFullCircle);
  ret->SetMarkerColor(color);
  ret->SetLineColor(color);
  ret->GetYaxis()->SetRangeUser(.5, 1.5);
  ret->GetYaxis()->SetTitle("Data / MC");
  ret->GetXaxis()->SetNdivisions(1900, kFALSE);
  //  ret->GetXaxis()->SetTickLength(1.);
  return ret;
}

void Zero() {
  TGraph* g = new TGraph();
  g->SetPoint(0, -1e10,0.003);
  g->SetPoint(1, 1e10, 0.003);
  g->SetLineWidth(2);
  g->SetLineStyle(1);
  g->Draw("same");

}

void OneLine(){
  TGraph* g = new TGraph();
  g->SetPoint(0, -1e10,1);
  g->SetPoint(1, 1e10, 1);
  g->SetLineWidth(2);
  g->SetLineStyle(7);
  g->Draw("same");
}

void Legend(TH1* h_to_clone){
  std::vector <std::pair<int, TString>> entries = {{kBlack, "ND data"},
                                                   {kRed, "Total MC"},
                                                   {kAzure, "NC"},
                                                   {kGreen+2, "#nu_{#mu} CC"},
                                                   {kMagenta, "#nu_{e} CC"},
                                                   {kGreen-3, "#bar{#nu}_{#mu} CC"},
                                                   {kMagenta-2, "#bar{#nu}_{e} CC"},
                                                   {kBlack, "Uncorr. MC"},
                                                   {17, "Uncorr. MC"},
                                                   {1, "Decomp MC"}};
                                                       

  TLegend* leg = new TLegend(.82, .5, .99, .92);

  leg->SetFillColor(kWhite);
  leg->SetBorderSize(1);
  leg->SetMargin(.22);

  TH1* h = (TH1F*)h_to_clone->Clone();
  leg->SetTextSize(0.5*h->GetXaxis()->GetTitleSize());
  h->SetLineStyle(1);
  h->SetLineColor(14);
  h->SetMarkerStyle(kFullCircle);
  
  leg->AddEntry(h, entries[0].second, "lep");
  for(int i = 1; i < (int) entries.size()-2; i++){
    h->SetLineColor(entries[i].first);
    if(i == 7) h->SetLineStyle(7); // Uncorr MC
    else h->SetLineStyle(1);
    if(i == 8){
      h->SetMarkerColor(entries[i].first); // Ratio to Uncorr MC
      h->SetMarkerStyle(kFullCircle);
      leg->AddEntry(h->Clone(), entries[i].second, "lep");
      continue;
    }
    if(i == 9) {
      h->SetMarkerColor(entries[i].first);
      h->SetMarkerStyle(kFullCircle);
      leg->AddEntry(h->Clone(), entries[i].second, "lep");
      continue;
    }
    leg->AddEntry(h->Clone(), entries[i].second, "l");                   
  }
  leg->SetTextSize(1.3*leg->GetTextSize());
  leg->Draw();

}

void RatioLegend(){
  TLegend* leg = new TLegend(.12, .13, .47, .18);
  //  leg->SetFillStyle(0);
  leg->SetNColumns(2);
  leg->SetBorderSize(1);
  TH1* h = new TH1F("","", 1, 0, 1);
  h->SetMarkerColor(17);
  h->SetLineColor(17);
  h->SetMarkerStyle(kFullCircle);
  leg->AddEntry(h->Clone(), "Uncorr. MC", "lep");
  h->SetMarkerColor(1);
  h->SetLineColor(1);
  leg->AddEntry(h->Clone(), "Decomp MC", "lep");
  leg->Draw();
}

void NeutrinoMode(TString mode){
  TString txt = mode == "fhc"? "Neutrino Mode" : "Antineutrino Mode";
  TLatex* lb = new TLatex(.1, .95, txt);
  lb->SetNDC();
  lb->SetTextAlign(13);
  lb->SetTextColor(kGray+3);
  lb->SetTextSize(2/45.);
  lb->Draw();

}