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