Analyse_GetEfficiency.C
Go to the documentation of this file.
5 
7 #include "CAFAna/Core/Spectrum.h"
9 
17 #include "CAFAna/Cuts/SpillCuts.h"
18 #include "CAFAna/Cuts/TimingCuts.h"
19 
21 
22 #include "TFile.h"
23 #include "TKey.h"
24 #include "TH2.h"
25 #include "TStyle.h"
26 #include "TCanvas.h"
27 #include "TLegend.h"
28 #include "TGraph.h"
29 
30 #include <fstream>
31 
32 
33 using namespace ana;
34 
35 
36 
37 std::pair<TH1*,TH1*> takeRatios(TH1 *h_MCData_Num, TH1 *h_MCData_Denom, TH1 *h_MCMC_Num, TH1 *h_MCMC_Denom)
38 {
39  /*
40  h_MCData_Num ->Rebin(3);
41  h_MCData_Denom->Rebin(3);
42  h_MCMC_Num ->Rebin(3);
43  h_MCMC_Denom ->Rebin(3);
44  */
45 
46  TH1 *h_MCData_Rat = (TH1*)h_MCData_Num->Clone();
47  h_MCData_Rat->Reset();
48  TH1 *h_MCMC_Rat = (TH1*)h_MCMC_Num->Clone();
49  h_MCMC_Rat->Reset();
50 
51  h_MCData_Rat->Divide(h_MCData_Num, h_MCData_Denom, 1, 1, "B");
52  h_MCMC_Rat ->Divide(h_MCMC_Num, h_MCMC_Denom, 1, 1, "B");
53 
54  std::pair<TH1*,TH1*> h_Rats = {h_MCData_Rat, h_MCMC_Rat};
55 
56  return h_Rats;
57 }
58 
59 
60 
61 TCanvas* makeCanvas(TH1 *h_MCData, TH1 *h_MCMC, TString s_Obj)
62 {
63  Helper help;
64 
65  h_MCData->SetLineColor(kBlack);
66  h_MCData->SetLineWidth(2);
67  h_MCData->GetXaxis()->SetTitle("");
68  h_MCData->GetYaxis()->SetTitle("Efficiency");
69  h_MCData->GetYaxis()->SetTitleSize(0.05);
70  h_MCData->GetYaxis()->SetTitleOffset(1.2);
71  h_MCData->GetYaxis()->CenterTitle();
72  h_MCMC->SetLineColor(kRed);
73  h_MCMC->SetLineWidth(2);
74  double ymax = std::max(h_MCMC->GetMaximum(), h_MCData->GetMaximum());
75 
76  TCanvas *c = new TCanvas("c", "c", 800, 800);
77  TPad *p1 = new TPad("p1", "p1", 0, 0, 1, 1);
78  p1->SetFillStyle(0);
79  p1->SetBottomMargin(0.3);
80  p1->SetTopMargin(0.01);
81  p1->SetRightMargin(0.05);
82  p1->SetLeftMargin(0.13);
83  TPad *p2 = new TPad("p2", "p2", 0, 0, 1, 1);
84  p2->SetFillStyle(0);
85  p2->SetTopMargin(0.7);
86  p2->SetRightMargin(0.05);
87  p2->SetBottomMargin(0.1);
88  p2->SetLeftMargin(0.13);
89 
90  p1->cd();
91 
92  h_MCData->Draw();
93  h_MCMC->Draw("SAME");
94  h_MCData->GetYaxis()->SetRangeUser(0, ymax + ymax*0.3);
95 
96  p2->cd();
97 
98  TH1D *h_Rat = (TH1D*)h_MCData->Clone();
99  h_MCData->GetXaxis()->SetLabelSize(0);
100  h_Rat->Divide(h_MCMC);
101  h_Rat->SetAxisRange(0.89,1.11,"Y");
102  h_Rat->SetLineColor(kBlue);
103  h_Rat->GetXaxis()->SetTitle(help.VarLabel((std::string)s_Obj, true).c_str());
104  h_Rat->GetYaxis()->SetTitle("MCData/MCMC");
105  h_Rat->GetYaxis()->SetTitleSize(0.05);
106  h_Rat->GetXaxis()->CenterTitle();
107  h_Rat->GetYaxis()->CenterTitle();
108  h_Rat->Draw();
109 
110  double x[2] = {h_Rat->GetXaxis()->GetXmin(), h_Rat->GetXaxis()->GetXmax()};
111  double y[2] = {1.0, 1.0};
112  TGraph* l_Xaxis = new TGraph(2, x, y);
113  l_Xaxis->SetLineColor(kGray);
114  l_Xaxis->SetLineWidth(3);
115  l_Xaxis->SetLineStyle(2);
116  l_Xaxis->Draw("L SAME");
117 
118  double x1_leg; double x2_leg;
119  s_Obj.Contains("kPDG") ? x1_leg = 0.4 : x1_leg = 0.65;
120  s_Obj.Contains("kPDG") ? x2_leg = 0.6 : x2_leg = 0.89;
121  TLegend *leg = new TLegend(x1_leg, 0.84, x2_leg, 0.98);
122  leg->AddEntry(h_MCData, "MC singles in Data", "L");
123  leg->AddEntry(h_MCMC, "MC singles in MC", "L");
124 
125  c->cd();
126  p1->Draw();
127  p2->Draw();
128  leg->Draw();
129 
130  return c;
131 }
132 
133 
134 double getEffDiffPc(std::map<std::string,std::pair<TH1*,TH1*>> map_HistPair, std::string const &cut, std::string const &s_Denom)
135 {
136  double effDiff = map_HistPair[cut].first ->Integral()/map_HistPair[s_Denom].first->Integral() - map_HistPair[cut].second->Integral()/map_HistPair[s_Denom].second->Integral();
137  effDiff /= map_HistPair[cut].second->Integral()/map_HistPair[s_Denom].second->Integral();
138 
139  if(std::abs(effDiff)<1.e-6) return 0.;
140 
141  return 100.*effDiff;
142 }
143 
144 
145 double getEffDiffErrPc(std::map<std::string,std::pair<TH1*,TH1*>> map_HistPair, std::string const &cut, std::string const &s_Denom)
146 {
147  double effDiffErrPc = 0.;
148 
149  /*OLD METHOD
150  double nSelectedData = map_HistPair[cut].first ->Integral();
151  double nSelectedMC = map_HistPair[cut].second->Integral();
152 
153  double pData = nSelectedData/map_HistPair[s_Denom].first ->Integral();
154  double pMC = nSelectedMC /map_HistPair[s_Denom].second->Integral();
155 
156  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));
157  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));
158 
159  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;
160  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;;
161  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;
162  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;;
163 
164  double sdSelectedData = std::sqrt(varSelectedData);
165  double sdSelectedMC = std::sqrt(varSelectedMC );
166 
167  double sdSelectedFracData = sdSelectedData*100/nSelectedData;
168  double sdSelectedFracMC = sdSelectedMC*100 /nSelectedMC;
169 
170  double varDelta = sdSelectedFracData*sdSelectedFracData + sdSelectedFracMC*sdSelectedFracMC;
171  double sdDelta = std::sqrt(varDelta);
172 
173  effDiffErrPc = sdDelta;
174  */
175 
176  //NEW METHOD
177  double nSelectedData = map_HistPair[cut].first ->Integral();
178  double nSelectedMC = map_HistPair[cut].second->Integral();
179 
180  double pData = nSelectedData/map_HistPair[s_Denom].first ->Integral();
181  double pMC = nSelectedMC /map_HistPair[s_Denom].second->Integral();
182 
183  double pVarData = pData*(1-pData)/map_HistPair[s_Denom].first ->Integral();
184  double pVarMC = pMC *(1-pMC) /map_HistPair[s_Denom].second->Integral();
185 
186  double pSdData = std::sqrt(pVarData);
187  double pSdMC = std::sqrt(pVarMC);
188 
189  double varDelta = pSdData*pSdData + pSdMC*pSdMC;
190  double sdDelta = std::sqrt(varDelta);
191 
192  double sdDeltaFrac = sdDelta/pMC;
193 
194  //MULTIPLY BY 100 TO MAKE A %.
195  effDiffErrPc = sdDeltaFrac*100.;
196 
197  return effDiffErrPc;
198 }
199 
200 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)
201 {
202  Helper help;
203 
204  gStyle->SetOptStat(0);
205 
206  TFile *f_InFileMCData = new TFile(s_InFileMCData.c_str(), "READ");
207  TFile *f_InFileMCMC = new TFile(s_InFileMCMC .c_str(), "READ");
208 
209  std::map<std::string,const Var> map_Vars;
210  std::map<std::string,const NuTruthVar> map_NuVars;
211  scheduleVars(map_Vars, map_NuVars, s_FileAppend);
212 
213  std::vector<std::string> vec_Cuts = help.CutFlow(s_FileAppend, useAnaTrueECut, getDetailed);
214  std::string firstEntry = vec_Cuts.at(0);
215  std::string lastEntry = vec_Cuts.at(vec_Cuts.size()-1);
216  std::map<std::string,std::pair<TH1*,TH1*>> map_HistPair;
217  std::map<std::string,std::pair<TH1*,TH1*>> map_HistPair_AllVars;
218  TIter iter(f_InFileMCData->GetListOfKeys());
219  TKey *key;
220  while((key=(TKey*)iter())){
221  TString s_Obj = ((TObject*)key->ReadObj())->GetName();
222 
223  for(auto var : map_Vars){
224  std::string s_EffVar = var.first;
225  if(!makePlots) s_FileAppend.find("nue")!=std::string::npos ? s_EffVar = "kTrueEben" : s_EffVar = "kTrueE";
226  for(auto cut : vec_Cuts){
227  std::string s_Spec = "dir_-CUT-" + cut + "-VAR-"+ s_EffVar;
228  std::string s_NuSpec = "dir_-CUT-" + cut + "-VAR-"+ s_EffVar + "_Nu";
229 
230  if(((std::string)s_Obj==s_Spec) || ((std::string)s_Obj==s_NuSpec)){
231 
232  std::unique_ptr<Spectrum> spec_MCData = Spectrum::LoadFrom(f_InFileMCData, s_Obj);
233  std::unique_ptr<Spectrum> spec_MCMC = Spectrum::LoadFrom(f_InFileMCMC , s_Obj);
234 
235  TH1 *h_MCData = spec_MCData->ToTH1(spec_MCData->POT());
236  h_MCData->SetName(s_Obj);
237  TH1 *h_MCMC = spec_MCMC ->ToTH1(spec_MCMC ->POT());
238  h_MCMC->SetName(s_Obj);
239  map_HistPair[cut] = {h_MCData, h_MCMC};
240  if((cut==firstEntry) || (cut==lastEntry)){
241  map_HistPair_AllVars[(std::string)s_Obj] = {h_MCData, h_MCMC};
242  }
243 
244  break;
245  }
246  }
247  if(!makePlots) break;
248  }
249  }
250 
251 
252  //SPECIAL PAIRS.
253  std::map<std::string,std::pair<std::string,std::string>> map_SpecialPairs;
254  std::map<std::string,std::vector<TH1*>> map_SpecialHists;
255  map_SpecialPairs["kSpillPOT"] = {"dir_-CUT-" + lastEntry + "-VAR-kSpillPOT","dir_-CUT-" + firstEntry + "-VAR-kSpillPOTInNuTruth_Nu"};
256 
257  for(auto pair : map_SpecialPairs){
258  std::unique_ptr<Spectrum> spec_MCData_Num = Spectrum::LoadFrom(f_InFileMCData, (TString)pair.second.first );
259  std::unique_ptr<Spectrum> spec_MCData_Denom = Spectrum::LoadFrom(f_InFileMCData, (TString)pair.second.second);
260  std::unique_ptr<Spectrum> spec_MCMC_Num = Spectrum::LoadFrom(f_InFileMCMC , (TString)pair.second.first );
261  std::unique_ptr<Spectrum> spec_MCMC_Denom = Spectrum::LoadFrom(f_InFileMCMC , (TString)pair.second.second);
262 
263  std::vector<TH1*> vec_Temp;
264  vec_Temp.push_back(spec_MCData_Num ->ToTH1(spec_MCData_Num ->POT()));
265  vec_Temp.at(vec_Temp.size()-1)->SetName(pair.second.first.c_str());
266  vec_Temp.push_back(spec_MCData_Denom->ToTH1(spec_MCData_Denom->POT()));
267  vec_Temp.at(vec_Temp.size()-1)->SetName(pair.second.second.c_str());
268  vec_Temp.push_back(spec_MCMC_Num ->ToTH1(spec_MCMC_Num ->POT()));
269  vec_Temp.at(vec_Temp.size()-1)->SetName(pair.second.first.c_str());
270  vec_Temp.push_back(spec_MCMC_Denom ->ToTH1(spec_MCMC_Denom ->POT()));
271  vec_Temp.at(vec_Temp.size()-1)->SetName(pair.second.second.c_str());
272 
273  map_SpecialHists[pair.first] = vec_Temp;
274  }
275 
276 
277  if(useAnaTrueECut) s_FileAppend += "_ATE";
278  if(getDetailed) s_FileAppend += "_Detailed";
279  /*
280  s_FileAppend += "_Rebin";
281  */
282  if(!makePlots){
283  std::ofstream f_Out;
284  f_Out.open(s_OutDir + "/CutFlowTable_" + s_FileAppend + ".tex");
285 
286  f_Out << "\\documentclass{article}" << std::endl;
287  f_Out << "\\usepackage[a4paper,margin=1in,landscape]{geometry}"<< std::endl;
288  f_Out << "\\usepackage{multirow}" << std::endl;
289  f_Out << "\\usepackage{array}" << std::endl;
290  f_Out << "\\usepackage{amsmath}" << std::endl;
291  f_Out << "\\newcolumntype{A}{>{\\centering\\arraybackslash} m{3cm}}" << std::endl;
292  f_Out << "\\newcolumntype{B}{>{\\centering\\arraybackslash} m{4.2cm}}" << std::endl;
293  f_Out << std::endl;
294  f_Out << "\\begin{document}" << std::endl;
295  f_Out << "\\begin{center}" << std::endl;
296  f_Out << "\\begin{tabular}{|A|A|A|A|A|B|}" << std::endl;
297  f_Out << "\\hline" << std::endl;
298  f_Out << "\\multirow{2}{*}{}&\\multicolumn{2}{c|}{\\textbf{MC in Data}} & \\multicolumn{2}{c|}{\\textbf{MC in MC}} &" << std::endl;
299  f_Out << "\\multirow{2}{*}{$\\frac{\\left(\\textbf{Data}_\\textbf{eff} - \\textbf{MC}_\\textbf{eff}\\right)}{\\textbf{MC}_\\textbf{eff}}$, \\%}\\\\" << std::endl;
300  f_Out << "\\cline{2-5}" << std::endl;
301  f_Out << "& \\textbf{Events} & \\textbf{Eff, \\%} & \\textbf{Events} & \\textbf{Eff, \\%} & \\\\" << std::endl;
302  f_Out << "\\hline" << std::endl;
303 
304  std::string s_Denom = vec_Cuts.at(0);
305  for(unsigned int cut = 0; cut < vec_Cuts.size(); cut++){
306  f_Out << "\\textbf{" << help.getNiceCutLabel(vec_Cuts.at(cut)) << "} & & & & & \\\\" << std::endl;
307  f_Out << " & " << map_HistPair[vec_Cuts.at(cut)].first->Integral();
308  f_Out << " & " << map_HistPair[vec_Cuts.at(cut)].first->Integral()*100./map_HistPair[s_Denom].first->Integral();
309  f_Out << " & " << map_HistPair[vec_Cuts.at(cut)].second->Integral();
310  f_Out << " & " << map_HistPair[vec_Cuts.at(cut)].second->Integral()*100./map_HistPair[s_Denom].second->Integral();
311  f_Out << " & \\large " << getEffDiffPc(map_HistPair, vec_Cuts.at(cut), s_Denom);
312  if(cut==vec_Cuts.size()-1){
313  f_Out << " $\\pm$ " << getEffDiffErrPc(map_HistPair, vec_Cuts.at(cut), s_Denom);
314  }
315  f_Out << "\\\\" << std::endl;
316  f_Out << "\\hline" << std::endl;
317  }
318 
319  f_Out << "\\end{tabular}" << std::endl;
320  f_Out << "\\end{center}" << std::endl;
321  f_Out << "\\end{document}" << std::endl;
322  f_Out.close();
323  }
324  else{
325  for(auto hist : map_HistPair_AllVars){
326  std::string s_Key = hist.first;
327  if(hist.first.find("CUT-"+firstEntry+"-VAR")!=std::string::npos){
328  std::string s_Denom = hist.first;
329  std::string s_Num = s_Denom;
330  s_Num.replace(s_Denom.find(firstEntry), firstEntry.length(), lastEntry);
331  s_Num = s_Num.substr(0, s_Num.length()-3);
332 
333  TH1 *h_MCData_Num = map_HistPair_AllVars[s_Num].first;
334  TH1 *h_MCData_Denom = hist.second.first;
335  TH1 *h_MCMC_Num = map_HistPair_AllVars[s_Num].second;
336  TH1 *h_MCMC_Denom = hist.second.second;
337 
338  std::pair<TH1*,TH1*> h_Rats = takeRatios(h_MCData_Num, h_MCData_Denom, h_MCMC_Num, h_MCMC_Denom);
339  TCanvas *c = makeCanvas(h_Rats.first, h_Rats.second, s_Num);
340 
341  std::string s_OutName = s_OutDir + "/GetEfficiency_" + help.VarLabel((std::string)s_Key, false) + "_" + s_FileAppend + ".pdf";
342  c->SaveAs(s_OutName.c_str());
343 
344  delete c;
345  }
346  }
347 
348  //SPECIAL PAIRS.
349  for(auto hists : map_SpecialHists){
350  std::pair<TH1*,TH1*> h_Rats = takeRatios(hists.second.at(0), hists.second.at(1), hists.second.at(2), hists.second.at(3));
351  TCanvas *c = makeCanvas(h_Rats.first, h_Rats.second, hists.second.at(1)->GetName());
352 
353  if(((TString)hists.second.at(1)->GetName()).Contains("kSpillPOT")){
354  TBox *box = new TBox(0.955, 0.1, 1.0, 0.2);
355  box->SetFillColor(kWhite);
356  box->Draw();
357  }
358 
359  std::string s_OutName = s_OutDir + "/GetEfficiency_" + help.VarLabel((std::string)hists.second.at(1)->GetName(), false) + "_" + s_FileAppend + ".pdf";
360  c->SaveAs(s_OutName.c_str());
361 
362  delete c;
363  }
364  }
365 
366  std::cout << "ANALYSIS DONE" << std::endl;
367 
368  return;
369 }
370 
371 
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)
T max(const caf::Proxy< T > &a, T b)
std::vector< std::string > CutFlow(std::string const &s_Selection, bool useAnaTrueECut, bool getDetailed)
Definition: Helper.h:284
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double Integral(const Spectrum &s, const double pot, int cvnregion)
double getEffDiffPc(std::map< std::string, std::pair< TH1 *, TH1 * >> map_HistPair, std::string const &cut, std::string const &s_Denom)
void makePlots()
Definition: makePlots.C:17
T sqrt(T number)
Definition: d0nt_math.hpp:156
TString hists[nhists]
Definition: bdt_com.C:3
float abs(float number)
Definition: d0nt_math.hpp:39
Double_t ymax
Definition: plot.C:25
std::string getNiceCutLabel(std::string const &s_Cut)
Definition: Helper.h:548
void scheduleVars(std::map< std::string, const Var > &map_Vars, std::map< std::string, const NuTruthVar > &map_NuVars, std::string s_Selection)
Definition: Scheduling.h:97
std::string GetName(int i)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
std::pair< TH1 *, TH1 * > takeRatios(TH1 *h_MCData_Num, TH1 *h_MCData_Denom, TH1 *h_MCMC_Num, TH1 *h_MCMC_Denom)
std::string VarLabel(std::string s_DirName, bool isNice)
Definition: Helper.h:222
double getEffDiffErrPc(std::map< std::string, std::pair< TH1 *, TH1 * >> map_HistPair, std::string const &cut, std::string const &s_Denom)
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
const Cut cut
Definition: exporter_fd.C:30
TCanvas * makeCanvas(TH1 *h_MCData, TH1 *h_MCMC, TString s_Obj)
Definition: Helper.h:14
Float_t e
Definition: plot.C:35
enum BeamMode kBlue
enum BeamMode string