Analyse_GetEfficiency.C

Go to the documentation of this file.

#include "3FlavorAna/Ana2020/SingleNuOverlay/Scheduling.h"
#include "3FlavorAna/Ana2020/SingleNuOverlay/Helper.h"
#include "3FlavorAna/Cuts/OverlayCuts.h"
#include "3FlavorAna/Vars/OverlayVars.h"

#include "CAFAna/Core/EventList.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"

#include "3FlavorAna/Cuts/BeamNueCuts.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/NumuCuts2020.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"

#include "3FlavorAna/Vars/NumuVars.h"

#include "TFile.h"
#include "TKey.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TGraph.h"

#include <fstream>


using namespace ana;



std::pair<TH1*,TH1*> takeRatios(TH1 *h_MCData_Num, TH1 *h_MCData_Denom, TH1 *h_MCMC_Num, TH1 *h_MCMC_Denom)
{
  /*
  h_MCData_Num  ->Rebin(3);
  h_MCData_Denom->Rebin(3);
  h_MCMC_Num    ->Rebin(3);
  h_MCMC_Denom  ->Rebin(3);
  */

  TH1 *h_MCData_Rat = (TH1*)h_MCData_Num->Clone();
  h_MCData_Rat->Reset();
  TH1 *h_MCMC_Rat   = (TH1*)h_MCMC_Num->Clone();
  h_MCMC_Rat->Reset();

  h_MCData_Rat->Divide(h_MCData_Num, h_MCData_Denom, 1, 1, "B");
  h_MCMC_Rat  ->Divide(h_MCMC_Num,   h_MCMC_Denom,   1, 1, "B");

  std::pair<TH1*,TH1*> h_Rats = {h_MCData_Rat, h_MCMC_Rat};

  return h_Rats;
}



TCanvas* makeCanvas(TH1 *h_MCData, TH1 *h_MCMC, TString s_Obj)
{
  Helper help;

  h_MCData->SetLineColor(kBlack);
  h_MCData->SetLineWidth(2);
  h_MCData->GetXaxis()->SetTitle("");
  h_MCData->GetYaxis()->SetTitle("Efficiency");
  h_MCData->GetYaxis()->SetTitleSize(0.05);
  h_MCData->GetYaxis()->SetTitleOffset(1.2);
  h_MCData->GetYaxis()->CenterTitle();
  h_MCMC->SetLineColor(kRed);
  h_MCMC->SetLineWidth(2);
  double ymax = std::max(h_MCMC->GetMaximum(), h_MCData->GetMaximum());

  TCanvas *c = new TCanvas("c", "c", 800, 800);
  TPad   *p1 = new TPad("p1", "p1", 0, 0, 1, 1);
  p1->SetFillStyle(0);
  p1->SetBottomMargin(0.3);
  p1->SetTopMargin(0.01);
  p1->SetRightMargin(0.05);
  p1->SetLeftMargin(0.13);
  TPad   *p2 = new TPad("p2", "p2", 0, 0, 1, 1);
  p2->SetFillStyle(0);
  p2->SetTopMargin(0.7);
  p2->SetRightMargin(0.05);
  p2->SetBottomMargin(0.1);
  p2->SetLeftMargin(0.13);

  p1->cd();

  h_MCData->Draw();
  h_MCMC->Draw("SAME");
  h_MCData->GetYaxis()->SetRangeUser(0, ymax + ymax*0.3);

  p2->cd();

  TH1D *h_Rat = (TH1D*)h_MCData->Clone();
  h_MCData->GetXaxis()->SetLabelSize(0);
  h_Rat->Divide(h_MCMC);
  h_Rat->SetAxisRange(0.89,1.11,"Y");
  h_Rat->SetLineColor(kBlue);
  h_Rat->GetXaxis()->SetTitle(help.VarLabel((std::string)s_Obj, true).c_str());
  h_Rat->GetYaxis()->SetTitle("MCData/MCMC");
  h_Rat->GetYaxis()->SetTitleSize(0.05);
  h_Rat->GetXaxis()->CenterTitle();
  h_Rat->GetYaxis()->CenterTitle();
  h_Rat->Draw();

  double x[2] = {h_Rat->GetXaxis()->GetXmin(), h_Rat->GetXaxis()->GetXmax()};
  double y[2] = {1.0, 1.0};
  TGraph* l_Xaxis = new TGraph(2, x, y);
  l_Xaxis->SetLineColor(kGray);
  l_Xaxis->SetLineWidth(3);
  l_Xaxis->SetLineStyle(2);
  l_Xaxis->Draw("L SAME");

  double x1_leg; double x2_leg;
  s_Obj.Contains("kPDG") ? x1_leg = 0.4 : x1_leg = 0.65;
  s_Obj.Contains("kPDG") ? x2_leg = 0.6 : x2_leg = 0.89;
  TLegend *leg = new TLegend(x1_leg, 0.84, x2_leg, 0.98);
  leg->AddEntry(h_MCData, "MC singles in Data", "L");
  leg->AddEntry(h_MCMC,   "MC singles in MC",   "L");

  c->cd();
  p1->Draw();
  p2->Draw();
  leg->Draw();

  return c;
}


double getEffDiffPc(std::map<std::string,std::pair<TH1*,TH1*>> map_HistPair, std::string const &cut, std::string const &s_Denom)
{
  double effDiff = map_HistPair[cut].first ->Integral()/map_HistPair[s_Denom].first->Integral() - map_HistPair[cut].second->Integral()/map_HistPair[s_Denom].second->Integral();
  effDiff       /= map_HistPair[cut].second->Integral()/map_HistPair[s_Denom].second->Integral();

  if(std::abs(effDiff)<1.e-6) return 0.;

  return 100.*effDiff;
}


double getEffDiffErrPc(std::map<std::string,std::pair<TH1*,TH1*>> map_HistPair, std::string const &cut, std::string const &s_Denom)
{
  double effDiffErrPc = 0.;

  /*OLD METHOD
  double nSelectedData = map_HistPair[cut].first ->Integral();
  double nSelectedMC   = map_HistPair[cut].second->Integral();

  double pData = nSelectedData/map_HistPair[s_Denom].first ->Integral();
  double pMC   = nSelectedMC  /map_HistPair[s_Denom].second->Integral();

  double varSelectedData  = pData*(1.-pData)*(map_HistPair[s_Denom].first ->GetSumw2()->GetSum() - map_HistPair[s_Denom].first ->GetSumw2()->GetAt(map_HistPair[s_Denom].first ->GetSumw2()->GetSize()-1));
  double varSelectedMC    = pMC*(1-pMC)     *(map_HistPair[s_Denom].second->GetSumw2()->GetSum() - map_HistPair[s_Denom].second->GetSumw2()->GetAt(map_HistPair[s_Denom].second->GetSumw2()->GetSize()-1));

  std::cout << "SSqW in Data: " << map_HistPair[s_Denom].first ->GetSumw2()->GetSum() - map_HistPair[s_Denom].first ->GetSumw2()->GetAt(map_HistPair[s_Denom].first ->GetSumw2()->GetSize()-1) << std::endl;
  std::cout << "SSqW in MC: "   << map_HistPair[s_Denom].second->GetSumw2()->GetSum() - map_HistPair[s_Denom].second->GetSumw2()->GetAt(map_HistPair[s_Denom].second->GetSumw2()->GetSize()-1) << std::endl;;
  std::cout << "SSqW in Data: " << map_HistPair[cut].first ->GetSumw2()->GetSum() - map_HistPair[cut].first ->GetSumw2()->GetAt(map_HistPair[cut].first ->GetSumw2()->GetSize()-1) << std::endl;
  std::cout << "SSqW in MC: "   << map_HistPair[cut].second->GetSumw2()->GetSum() - map_HistPair[cut].second->GetSumw2()->GetAt(map_HistPair[cut].second->GetSumw2()->GetSize()-1) << std::endl;;

  double sdSelectedData = std::sqrt(varSelectedData);
  double sdSelectedMC   = std::sqrt(varSelectedMC  );

  double sdSelectedFracData = sdSelectedData*100/nSelectedData;
  double sdSelectedFracMC   = sdSelectedMC*100  /nSelectedMC;

  double varDelta = sdSelectedFracData*sdSelectedFracData + sdSelectedFracMC*sdSelectedFracMC;
  double sdDelta  = std::sqrt(varDelta);

  effDiffErrPc = sdDelta;
  */

  //NEW METHOD
  double nSelectedData = map_HistPair[cut].first ->Integral();
  double nSelectedMC   = map_HistPair[cut].second->Integral();

  double pData = nSelectedData/map_HistPair[s_Denom].first ->Integral();
  double pMC   = nSelectedMC  /map_HistPair[s_Denom].second->Integral();

  double pVarData = pData*(1-pData)/map_HistPair[s_Denom].first ->Integral();
  double pVarMC   = pMC  *(1-pMC)  /map_HistPair[s_Denom].second->Integral();

  double pSdData = std::sqrt(pVarData);
  double pSdMC   = std::sqrt(pVarMC);

  double varDelta = pSdData*pSdData + pSdMC*pSdMC;
  double sdDelta  = std::sqrt(varDelta);

  double sdDeltaFrac = sdDelta/pMC;

  //MULTIPLY BY 100 TO MAKE A %.
  effDiffErrPc = sdDeltaFrac*100.;

  return effDiffErrPc;
}

void Analyse_GetEfficiency(std::string s_OutDir = "", std::string s_FileAppend = "", std::string s_InFileMCData = "", std::string s_InFileMCMC = "", bool useAnaTrueECut = true, bool makePlots = false, bool getDetailed = false)
{
  Helper help;

  gStyle->SetOptStat(0);

  TFile *f_InFileMCData = new TFile(s_InFileMCData.c_str(), "READ");
  TFile *f_InFileMCMC   = new TFile(s_InFileMCMC  .c_str(), "READ");

  std::map<std::string,const Var>        map_Vars; 
  std::map<std::string,const NuTruthVar> map_NuVars; 
  scheduleVars(map_Vars, map_NuVars, s_FileAppend);
  
  std::vector<std::string> vec_Cuts = help.CutFlow(s_FileAppend, useAnaTrueECut, getDetailed);
  std::string firstEntry = vec_Cuts.at(0);
  std::string lastEntry  = vec_Cuts.at(vec_Cuts.size()-1);
  std::map<std::string,std::pair<TH1*,TH1*>> map_HistPair;
  std::map<std::string,std::pair<TH1*,TH1*>> map_HistPair_AllVars;
  TIter iter(f_InFileMCData->GetListOfKeys());
  TKey *key;
  while((key=(TKey*)iter())){
    TString s_Obj = ((TObject*)key->ReadObj())->GetName();

    for(auto var : map_Vars){
      std::string s_EffVar = var.first;
      if(!makePlots) s_FileAppend.find("nue")!=std::string::npos ? s_EffVar = "kTrueEben" : s_EffVar = "kTrueE";
      for(auto cut : vec_Cuts){
        std::string s_Spec   = "dir_-CUT-" + cut + "-VAR-"+ s_EffVar;
        std::string s_NuSpec = "dir_-CUT-" + cut + "-VAR-"+ s_EffVar + "_Nu";

        if(((std::string)s_Obj==s_Spec) || ((std::string)s_Obj==s_NuSpec)){

          std::unique_ptr<Spectrum> spec_MCData = Spectrum::LoadFrom(f_InFileMCData, s_Obj);
          std::unique_ptr<Spectrum> spec_MCMC   = Spectrum::LoadFrom(f_InFileMCMC  , s_Obj);

          TH1 *h_MCData = spec_MCData->ToTH1(spec_MCData->POT());
          h_MCData->SetName(s_Obj);
          TH1 *h_MCMC   = spec_MCMC  ->ToTH1(spec_MCMC  ->POT());
          h_MCMC->SetName(s_Obj);
          map_HistPair[cut] = {h_MCData, h_MCMC};
          if((cut==firstEntry) || (cut==lastEntry)){
            map_HistPair_AllVars[(std::string)s_Obj] = {h_MCData, h_MCMC};
          }

          break;
        }
      }
      if(!makePlots) break;
    }
  }


  //SPECIAL PAIRS.
  std::map<std::string,std::pair<std::string,std::string>> map_SpecialPairs;
  std::map<std::string,std::vector<TH1*>> map_SpecialHists;
  map_SpecialPairs["kSpillPOT"] = {"dir_-CUT-" + lastEntry + "-VAR-kSpillPOT","dir_-CUT-" + firstEntry + "-VAR-kSpillPOTInNuTruth_Nu"};

  for(auto pair : map_SpecialPairs){
    std::unique_ptr<Spectrum> spec_MCData_Num   = Spectrum::LoadFrom(f_InFileMCData, (TString)pair.second.first );
    std::unique_ptr<Spectrum> spec_MCData_Denom = Spectrum::LoadFrom(f_InFileMCData, (TString)pair.second.second);
    std::unique_ptr<Spectrum> spec_MCMC_Num     = Spectrum::LoadFrom(f_InFileMCMC  , (TString)pair.second.first );
    std::unique_ptr<Spectrum> spec_MCMC_Denom   = Spectrum::LoadFrom(f_InFileMCMC  , (TString)pair.second.second);

    std::vector<TH1*> vec_Temp;
    vec_Temp.push_back(spec_MCData_Num  ->ToTH1(spec_MCData_Num  ->POT()));
    vec_Temp.at(vec_Temp.size()-1)->SetName(pair.second.first.c_str());
    vec_Temp.push_back(spec_MCData_Denom->ToTH1(spec_MCData_Denom->POT()));
    vec_Temp.at(vec_Temp.size()-1)->SetName(pair.second.second.c_str());
    vec_Temp.push_back(spec_MCMC_Num    ->ToTH1(spec_MCMC_Num    ->POT()));
    vec_Temp.at(vec_Temp.size()-1)->SetName(pair.second.first.c_str());
    vec_Temp.push_back(spec_MCMC_Denom  ->ToTH1(spec_MCMC_Denom  ->POT()));
    vec_Temp.at(vec_Temp.size()-1)->SetName(pair.second.second.c_str());

    map_SpecialHists[pair.first] = vec_Temp;
  }


  if(useAnaTrueECut) s_FileAppend += "_ATE";
  if(getDetailed)    s_FileAppend += "_Detailed";
  /*
  s_FileAppend += "_Rebin";
  */
  if(!makePlots){
    std::ofstream f_Out;
    f_Out.open(s_OutDir + "/CutFlowTable_" + s_FileAppend + ".tex");
  
    f_Out << "\\documentclass{article}"      << std::endl;
    f_Out << "\\usepackage[a4paper,margin=1in,landscape]{geometry}"<< std::endl;
    f_Out << "\\usepackage{multirow}"        << std::endl;
    f_Out << "\\usepackage{array}"           << std::endl;
    f_Out << "\\usepackage{amsmath}"         << std::endl;
    f_Out << "\\newcolumntype{A}{>{\\centering\\arraybackslash} m{3cm}}"   << std::endl;
    f_Out << "\\newcolumntype{B}{>{\\centering\\arraybackslash} m{4.2cm}}" << std::endl;
    f_Out << std::endl;
    f_Out << "\\begin{document}"             << std::endl;
    f_Out << "\\begin{center}"               << std::endl;
    f_Out << "\\begin{tabular}{|A|A|A|A|A|B|}" << std::endl;
    f_Out << "\\hline" << std::endl;
    f_Out << "\\multirow{2}{*}{}&\\multicolumn{2}{c|}{\\textbf{MC in Data}} & \\multicolumn{2}{c|}{\\textbf{MC in MC}} &" << std::endl;
    f_Out << "\\multirow{2}{*}{$\\frac{\\left(\\textbf{Data}_\\textbf{eff} - \\textbf{MC}_\\textbf{eff}\\right)}{\\textbf{MC}_\\textbf{eff}}$, \\%}\\\\" << std::endl;
    f_Out << "\\cline{2-5}" << std::endl;
    f_Out << "& \\textbf{Events} & \\textbf{Eff, \\%} & \\textbf{Events} & \\textbf{Eff, \\%} & \\\\" << std::endl;
    f_Out << "\\hline" << std::endl;
    
    std::string s_Denom = vec_Cuts.at(0);
    for(unsigned int cut = 0; cut < vec_Cuts.size(); cut++){
      f_Out << "\\textbf{" << help.getNiceCutLabel(vec_Cuts.at(cut)) << "} & & & & & \\\\" << std::endl;    
      f_Out << " & "  << map_HistPair[vec_Cuts.at(cut)].first->Integral();
      f_Out << " & "  << map_HistPair[vec_Cuts.at(cut)].first->Integral()*100./map_HistPair[s_Denom].first->Integral();
      f_Out << " & "  << map_HistPair[vec_Cuts.at(cut)].second->Integral();
      f_Out << " & "  << map_HistPair[vec_Cuts.at(cut)].second->Integral()*100./map_HistPair[s_Denom].second->Integral();
      f_Out << " & \\large "  << getEffDiffPc(map_HistPair, vec_Cuts.at(cut), s_Denom);
      if(cut==vec_Cuts.size()-1){
        f_Out << " $\\pm$ "  << getEffDiffErrPc(map_HistPair, vec_Cuts.at(cut), s_Denom);
      }
      f_Out << "\\\\" << std::endl;    
      f_Out << "\\hline" << std::endl;
    }
  
    f_Out << "\\end{tabular}"              << std::endl;
    f_Out << "\\end{center}"               << std::endl;
    f_Out << "\\end{document}"             << std::endl;
    f_Out.close();
  }
  else{
    for(auto hist : map_HistPair_AllVars){
      std::string s_Key = hist.first;
      if(hist.first.find("CUT-"+firstEntry+"-VAR")!=std::string::npos){
        std::string s_Denom = hist.first;
        std::string s_Num   = s_Denom;
        s_Num.replace(s_Denom.find(firstEntry), firstEntry.length(), lastEntry);
        s_Num = s_Num.substr(0, s_Num.length()-3);

        TH1 *h_MCData_Num   = map_HistPair_AllVars[s_Num].first;
        TH1 *h_MCData_Denom = hist.second.first;
        TH1 *h_MCMC_Num     = map_HistPair_AllVars[s_Num].second;
        TH1 *h_MCMC_Denom   = hist.second.second;

        std::pair<TH1*,TH1*> h_Rats = takeRatios(h_MCData_Num, h_MCData_Denom, h_MCMC_Num, h_MCMC_Denom);
        TCanvas *c = makeCanvas(h_Rats.first, h_Rats.second, s_Num);

        std::string s_OutName = s_OutDir + "/GetEfficiency_" + help.VarLabel((std::string)s_Key, false) + "_" + s_FileAppend + ".pdf";
        c->SaveAs(s_OutName.c_str());

        delete c;
      }
    }

    //SPECIAL PAIRS.
    for(auto hists : map_SpecialHists){
      std::pair<TH1*,TH1*> h_Rats = takeRatios(hists.second.at(0), hists.second.at(1), hists.second.at(2), hists.second.at(3));
      TCanvas *c = makeCanvas(h_Rats.first, h_Rats.second, hists.second.at(1)->GetName());

      if(((TString)hists.second.at(1)->GetName()).Contains("kSpillPOT")){
        TBox *box = new TBox(0.955, 0.1, 1.0, 0.2);
        box->SetFillColor(kWhite);
        box->Draw();
      }

      std::string s_OutName = s_OutDir + "/GetEfficiency_" + help.VarLabel((std::string)hists.second.at(1)->GetName(), false) + "_" + s_FileAppend + ".pdf";
      c->SaveAs(s_OutName.c_str());

      delete c;
    }
  }

  std::cout << "ANALYSIS DONE" << std::endl;

  return;
}