PlotPreSelectionPlots.C
Go to the documentation of this file.
1 // plots output of DrawPreSelectionPlots.C
2 
3 // defines whether the cut removes low values of the variable
4 // or high values of the variable
5 enum CutSide {
6  kLow,
10 };
11 
12 enum FomType {
13  kSOverSB, // signal over sqrt background
14  kSOverSSB, // signal over sqrt (signal plus background)
15  kEfficiency, // efficiency (total true selected / total true)
16  // -- n.b. "total true" is after the Concat cuts, so not quite a true purity
17  kPurity, // purity (total true selected / total selected)
18  kEffTimesPurity // efficiency times purity
19 };
20 
21 void printValues(std::vector<double> theseValues, std::vector<double> theseFOMValues){
22 
23  std::cout << "cut value: ";
24  for (int i = 0; i < theseValues.size(); ++i){
25  std::cout << " " << theseValues.at(i);
26  }
28 
29  std::cout << "FOM value: ";
30  for (int i = 0; i < theseFOMValues.size(); ++i){
31  std::cout << " " << theseFOMValues.at(i);
32  }
33  std::cout << "\n" << std::endl;
34 
35 }
36 
37 // find cut values where FOMs are maximised
38 void getMaximumFOM(TH1D* h){
39 
40  std::vector<double> cutValues;
41  std::vector<double> FOMValues;
42 
43  double maxFOM = h->GetMaximum();
44  for (int i = 1; i < h->GetNbinsX()+1; ++i)
45  {
46  if (h->GetBinContent(i) == maxFOM)
47  {
48  cutValues.push_back(h->GetBinLowEdge(i));
49  FOMValues.push_back(h->GetBinContent(i));
50  }
51  }
52 
53  printValues(cutValues, FOMValues);
54 
55 }
56 
57 // simple function to calculate SOverSB
58 double calcSOverSB(double sig, double bg){
59  if (sig == 0) return 0;
60  else return sig/std::sqrt(bg);
61 }
62 
63 // simple function to calculate SOverSSB
64 double calcSOverSSB(double sig, double bg, double fracUn){
65  if (sig == 0) return 0;
66  else return sig/std::sqrt(sig+bg+std::pow(fracUn,2));
67 }
68 
69 // calculate efficiency
70 double calcEfficiency(double selSig, double totSig){
71  if (totSig == 0) return 0;
72  else return selSig/totSig;
73 }
74 
75 // calculate purity
76 double calcPurity(double selSig, double selTot){
77  if (selTot == 0) return 0;
78  else return selSig/selTot;
79 }
80 
81 double getFOM(FomType thisFomType, double selSig, double selBg, double totSig, double selTot, double fracUn){
82 
83  switch (thisFomType){
84  case kSOverSB:
85  return calcSOverSB(selSig, selBg);
86  case kSOverSSB:
87  return calcSOverSSB(selSig, selBg, fracUn);
88  case kEfficiency:
89  return calcEfficiency(selSig, totSig);
90  case kPurity:
91  return calcPurity(selSig, selTot);
92  case kEffTimesPurity:
93  return calcEfficiency(selSig, totSig) * calcPurity(selSig, selTot);
94  default:
95  throw std::logic_error("thisFomType not accepted");
96  }
97 
98 }
99 
100 double getFracUncertainty(TH1D* h, std::vector<float> fracUncertainties){
101 
102  if (h->GetNbinsX() != fracUncertainties.size()+1){
103  std::cout << "hBins: " << h->GetNbinsX()
104  << " fracUncertainties: " << fracUncertainties.size() << std::endl;
105  }
106 
107  TH1D* hShifted = (TH1D*)h->Clone("hShifted");
108 
109  for (int i = 2; i < h->GetNbinsX()+1; ++i){
110  hShifted->SetBinContent(i,h->GetBinContent(i)*(1+fracUncertainties.at(i-2)));
111  }
112 
113 
114  double totUnc = (hShifted->Integral()/h->Integral())-1;
115 
116  return totUnc;
117 
118 }
119 
120 // function to return the histogram of the Figure of Merit
121 TH1D* getFOMHist(std::string fomName,
122  CutSide thisCutSide,
123  FomType thisFomType,
124  TH1D* all,
125  TH1D* sig,
126  TH1D* trueSig,
127  TH2D* energyCorr = nullptr,
128  std::vector<float> fracUncertainties = {}){
129 
130  TH1D* fomReturnHist = new TH1D(fomName.c_str(),
131  ";;FOM",
132  all->GetNbinsX(),
133  all->GetBinLowEdge(1),
134  all->GetBinLowEdge(all->GetNbinsX()+1));
135 
136  double selSig, totSig, selTot, selBg;
137  if (thisCutSide == CutSide::kBoth)
138  return fomReturnHist;
139  else if (thisCutSide == CutSide::kLow){
140  for (int i = 1; i < fomReturnHist->GetNbinsX()+1; i++){
141  selTot = all->Integral(i, fomReturnHist->GetNbinsX());
142  selSig = sig->Integral(i, fomReturnHist->GetNbinsX());
143  selBg = selTot - selSig;
144  totSig = trueSig->Integral();
145 
146  TH1D* proj = nullptr;
147  if (energyCorr != nullptr){
148  proj = energyCorr->ProjectionX("proj", i, fomReturnHist->GetNbinsX());
149  }
150  double fracUncertainty = (proj == nullptr) ? 0 : getFracUncertainty(proj, fracUncertainties);
151 
152  fomReturnHist->SetBinContent(i, getFOM(thisFomType, selSig, selBg, totSig, selTot, fracUncertainty));
153  }
154  }
155  else if (thisCutSide == CutSide::kHigh){
156  for (int i = fomReturnHist->GetNbinsX()+1; i > 0; i--){
157  selTot = all->Integral(1,i);
158  selSig = sig->Integral(1,i);
159  selBg = selTot - selSig;
160  totSig = trueSig->Integral();
161 
162  TH1D* proj = nullptr;
163  if (energyCorr != nullptr){
164  proj = energyCorr->ProjectionX("proj", 1, i);
165  }
166  double fracUncertainty = (proj == nullptr) ? 0 : getFracUncertainty(proj, fracUncertainties);
167 
168  fomReturnHist->SetBinContent(i, getFOM(thisFomType, selSig, selBg, totSig, selTot, fracUncertainty));
169 
170  }
171  }
172 
173  return fomReturnHist;
174 }
175 
176 // main function
178 
179  TFile* covmxFile = new TFile("/pnfs/nova/persistent/analysis/nux/nus20/cmf_covmx/2020-02-14/covariance_matrix_hist_all.root",
180  "read");
181 
182  TH2D* covmx = (TH2D*)covmxFile->Get("covar/covarianceMatrix");
183  std::vector<float> fracUncertainty;
184 
185  for (int i = 93; i < 93+38; ++i){
186  double mod = 1;
187  if (i > 93+0 && i < 93+20)
188  mod = 1;
189  fracUncertainty.push_back(std::sqrt(covmx->GetBinContent(i,i))*mod);
190  }
191 
192  std::vector< std::string > cutLevels = {
193  "NoCut" ,
194  "Quality",
195  "QualityVtx",
196  "QualityVtxFiducial",
197  "FullSelection"
198  };
199 
200  // first up define some vectors for stuff we're going to loop over
201  std::vector< std::string > histNames1D = {
202  "kHitsPerPlane",
203  "kNSliceHits",
204  "kPartPtp",
205  "kIsElastic",
206  "kNFuzzyProng",
207  "kNContPlanes",
208  "kVtxX",
209  "kVtxY",
210  "kVtxZ",
211  "kDistAllTop",
212  "kDistAllBottom",
213  "kDistAllEast",
214  "kDistAllWest",
215  "kDistAllFront",
216  "kDistAllBack",
217  "kCVNnc",
218  "kCVNnc_looseptp",
219  "kCVNnc_oldpresel",
220  "kCaloE"
221  };
222 
223  std::vector< std::string > histTitles1D = {
224  ";Hits Per Plane;Events;",
225  ";Number Hits in Slice;Events",
226  ";Beam-Transverse Momentum Fraction; Events",
227  ";Number of Elastic Vertices;Events",
228  ";Number of Fuzzy Prongs;Events",
229  ";Number of Continuous Planes;Events",
230  ";Vertex x (cm);Events",
231  ";Vertex y (cm);Events",
232  ";Vertex z (cm);Events",
233  ";Minimum Distance To Top;Events",
234  ";Minimum Distance To Bottom;Events",
235  ";Minimum Distance To East;Events",
236  ";Minimum Distance To West;Events",
237  ";Minimum Distance To Front;Events",
238  ";Minimum Distance To Back;Events",
239  ";CVNnc Score;Events",
240  ";CVNnc Score (looseptp);Events",
241  ";CVNnc Score (oldpresel);Events",
242  ";Calorimetric Energy (GeV);Events"
243  };
244 
245  std::vector< CutSide > histCutSide = {
246  kHigh , // kHitsPerPlane
247  kLow , // kNSliceHits
248  kHigh , // kPartPtp
249  kLow , // kIsElastic
250  kLow , // kNFuzzyProng
251  kLow , // kNContPlanes
252  kBoth , // kVtxX
253  kBoth , // kVtxY
254  kBoth , // kVtxZ
255  kLow , // kDistAllTopHist
256  kLow , // kDistAllBottomHist
257  kLow , // kDistAllEastHist
258  kLow , // kDistAllWestHist
259  kLow , // kDistAllFrontHist
260  kLow , // kDistAllBackHist
261  kLow , // kCVNnc
262  kLow , // kCVNnc_looseptpHist
263  kLow , // kCVNnc_oldpreselHist
264  kNoCut // kCaloE
265  };
266 
267  std::vector< bool > plotFOM = {
268  true,
269  true,
270  true,
271  false,
272  false,
273  false,
274  false,
275  false,
276  false,
277  false,
278  false,
279  false,
280  false,
281  false,
282  false,
283  true,
284  true,
285  true,
286  false
287  };
288 
289  std::vector< std::vector< double > > currentCuts = {
290  {0.0} , // kHitsPerPlane
291  {0.0} , // kNSliceHits
292  {1.0} , // kPartPtp
293  {1} , // kIsElastic
294  {1} , // kNFuzzyProng
295  {2} , // kNContPlanes
296  {-100, 100} , // kVtxX
297  {-100, 100} , // kVtxY
298  {150 , 1000}, // kVtxZ
299  {20} , // kDistAllTopHist
300  {20} , // kDistAllBottomHist
301  {20} , // kDistAllEastHist
302  {20} , // kDistAllWestHist
303  {150} , // kDistAllFrontHist
304  {50} , // kDistAllBackHist
305  {0} , // kCVNnc
306  {0.42} , // kCVNnc_looseptpHist
307  {0} , // kCVNnc_oldpreselHist
308  {0,0} // kCaloE
309  };
310 
311  std::vector< std::string > histNames2D = {
312  "kHitsPerPlaneNSliceHitsHist",
313  "kHitsPerPlanePartPtpHist",
314  "kNSliceHitsPartPtpHist"
315  };
316 
317  std::vector< std::string > histTitles2D = {
318  ";Hits Per Plane; Number of Hits in Slice",
319  ";Hits Per Plane; kPartPtp",
320  ";Number of Hits in Slice; kPartPtp"
321  };
322 
323  double maxVal = 1.5;
324  int textSize = 28;
325  // fist deal with the 1d histograms
326  for (int iH = 0; iH < histNames1D.size(); iH++){
327 
328  TCanvas* c1Ratio = new TCanvas("c1Ratio", "c1Ratio", 600, 600);
329  TPad *topPad = new TPad("topPad", "", 0.005, 0.3, 0.995, 0.995);
330  TPad *bottomPad = new TPad("bottomPad", "", 0.005, 0.005, 0.995, 0.3);
331  topPad ->SetTopMargin(0.1);
332  topPad ->SetBottomMargin(0.005);
333  bottomPad->SetTopMargin(0.005);
334  bottomPad->SetBottomMargin(0.35);
335  bottomPad->SetGridy();
336  topPad ->Draw();
337  bottomPad->Draw();
338 
339  TCanvas* c1 = new TCanvas("c1", "c1", 600, 500);
340  TPad *thisPad = new TPad("topPad", "", 0.005, 0.005, 0.995, 0.995);
341  thisPad ->SetTopMargin(0.1);
342  thisPad ->SetBottomMargin(0.15);
343  thisPad->Draw();
344  thisPad->cd();
345 
346  for (int icl = 0; icl < cutLevels.size(); ++icl){
347  std::cout << std::string("sim"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)) << std::endl;
348  TH1D* histSim = (TH1D*)_file0->Get(("sim"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)).c_str());
349  TH1D* histSimTrueNC = (TH1D*)_file0->Get(("sim"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)+"TrueNC").c_str());
350  TH1D* histSimTrueNCNoCut = (TH1D*)_file0->Get(("sim"+histNames1D.at(iH)+"Hist"+"NoCutTrueNC").c_str());
351  TH1D* histDatNoCut = (TH1D*)_file0->Get(("dat"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)).c_str());
352 
353  // get 2D histograms for variables correlated with energy
354  std::string name = histNames1D.at(iH);
355  std::string simEnergyString = "sim"+name+"CaloEHist"+cutLevels.at(icl);
356  TH2D* histSimEnergy = (TH2D*)_file0->Get(simEnergyString.c_str());
357 
358  double scalef = 1;
359  if (histDatNoCut != nullptr)
360  scalef = (double)histDatNoCut->Integral()/histSim->Integral();
361 
362  histSim -> Sumw2();
363  histSim -> Scale(scalef);
364  histSimTrueNC -> Sumw2();
365  histSimTrueNC -> Scale(scalef);
366  histSimTrueNCNoCut -> Sumw2();
367  histSimTrueNCNoCut -> Scale(scalef);
368 
369  histSim ->SetLineColor(kGray+1);
370  histSimTrueNC->SetFillColor(kAzure+1);
371  histSimTrueNC->SetLineColor(kAzure+1);
372  histSimTrueNC->SetMarkerColor(kAzure+1);
373 
374  histSim ->SetLineWidth(2);
375  histSim ->GetYaxis()->SetRangeUser(0.0001, histSim->GetMaximum()*maxVal);
376  histSim ->SetTitle(histTitles1D.at(iH).c_str());
377  histSim ->GetYaxis()->SetTitleOffset(1.25);
378  histSim ->GetYaxis()->SetMaxDigits(3);
379  histSim ->GetYaxis()->SetNdivisions(515);
380  histSim ->GetXaxis()->SetTitleFont(43);
381  histSim ->GetXaxis()->SetTitleSize(textSize);
382  histSim ->GetXaxis()->CenterTitle();
383  histSim ->GetYaxis()->SetTitleFont(43);
384  histSim ->GetYaxis()->SetTitleSize(textSize);
385  histSim ->GetYaxis()->CenterTitle();
386  histSim ->GetYaxis()->SetLabelFont(43);
387  histSim ->GetYaxis()->SetLabelSize(textSize);
388  histSim ->GetXaxis()->SetLabelFont(43);
389  histSim ->GetXaxis()->SetLabelSize(textSize);
390  bool isDrawData = false;
391  if ((histNames1D.at(iH).find("CVN") == std::string::npos &&
392  histNames1D.at(iH).find("CaloE") == std::string::npos) ||
393  (histNames1D.at(iH).find("CVN") != std::string::npos &&
394  cutLevels.at(icl).find("QualityVtxFiducial") != std::string::npos) ||
395  (histNames1D.at(iH).find("CVN") != std::string::npos &&
396  cutLevels.at(icl).find("FullSelection") != std::string::npos) ||
397  (histNames1D.at(iH).find("CaloE") != std::string::npos &&
398  cutLevels.at(icl).find("FullSelection") != std::string::npos))
399  isDrawData = true;
400 
401  if(isDrawData){
402  c1Ratio->cd();
403  topPad->cd();
404  histSim ->GetYaxis() ->SetTitleOffset(1.3);
405  histSim ->GetYaxis()->SetMaxDigits(3);
406  histSim ->GetYaxis()->SetNdivisions(515);
407 
408  }
409 
410  histSim ->Draw("hist");
411  histSimTrueNC ->Draw("same hist");
412  histSim ->Draw("same hist");
413  if (isDrawData){
414  histDatNoCut ->SetMarkerStyle(20);
415  histDatNoCut ->SetMarkerSize(0.6);
416 
417  histDatNoCut ->GetXaxis()->SetTitleFont(43);
418  histDatNoCut ->GetXaxis()->SetTitleSize(textSize);
419  histDatNoCut ->GetXaxis()->CenterTitle();
420  histDatNoCut ->GetYaxis()->SetTitleFont(43);
421  histDatNoCut ->GetYaxis()->SetTitleSize(textSize);
422  histDatNoCut ->GetYaxis()->CenterTitle();
423  histDatNoCut ->GetYaxis()->SetLabelFont(43);
424  histDatNoCut ->GetYaxis()->SetLabelSize(textSize);
425  histDatNoCut ->GetXaxis()->SetLabelFont(43);
426  histDatNoCut ->GetXaxis()->SetLabelSize(textSize);
427  histDatNoCut ->SetTitle(histTitles1D.at(iH).c_str());
428 
429 
430  histDatNoCut ->Draw("same p");
431  bottomPad->cd();
432  TH1D* ratioHist = (TH1D*)histDatNoCut->Clone("ratioHist");
433  ratioHist->Divide(histSim);
434  ratioHist->GetYaxis()->SetRangeUser(0.7, 1.29);
435  ratioHist->GetYaxis()->SetNdivisions(505);
436  ratioHist->GetYaxis()->SetTitle("Data/Sim");
437  ratioHist->GetYaxis()->SetTitleOffset(1.3);
438  ratioHist->GetXaxis()->SetTitleOffset(3.0);
439  ratioHist->Draw("p");
440  topPad->cd();
441  }
442  TH1D* EffTimesPurHist = (TH1D*)getFOMHist("EffTimesPur"+histNames1D.at(iH)+cutLevels.at(icl),
443  histCutSide.at(iH),
445  histSim,
446  histSimTrueNC,
447  histSimTrueNCNoCut);
448 
449  TH1D* SOverSSBHist = (TH1D*)getFOMHist("SOverSSB"+histNames1D.at(iH)+cutLevels.at(icl),
450  histCutSide.at(iH),
452  histSim,
453  histSimTrueNC,
454  histSimTrueNCNoCut,
455  histSimEnergy,
456  fracUncertainty);
457 
458  TH1D* SOverSBHist = (TH1D*)getFOMHist("SOverSB"+histNames1D.at(iH)+cutLevels.at(icl),
459  histCutSide.at(iH),
461  histSim,
462  histSimTrueNC,
463  histSimTrueNCNoCut);
464 
465  // define new axis on the right hand side for the EffTimesPurHist
466  TLegend* legFOM = new TLegend(0.5, 0.70, 0.8, 0.89);
467  if (plotFOM.at(iH) == true){
468 
469  std::cout
470  << "Printing cut values: "
471  << "\n-- Efficiency x Purity: "
472  << std::endl;
473  getMaximumFOM(EffTimesPurHist);
474  std::cout
475  << "\n-- SOverSSB:"
476  << std::endl;
477  getMaximumFOM(SOverSSBHist);
478  std::cout
479  << "\n-- SOverSB:"
480  << std::endl;
481  getMaximumFOM(SOverSBHist);
482 
483  double scaleval;
484  if (isDrawData)
485  scaleval = histDatNoCut->GetMaximum();
486  else
487  scaleval = histSim->GetMaximum()/maxVal;
488 
489  EffTimesPurHist->GetYaxis()->SetRangeUser(0,EffTimesPurHist->GetMaximum()*maxVal);
490  EffTimesPurHist->Scale(scaleval/EffTimesPurHist->GetMaximum());
491  EffTimesPurHist->SetLineColor(kGreen-3);
492  EffTimesPurHist->SetLineWidth(2);
493  EffTimesPurHist->Draw("same hist");
494 
495  SOverSSBHist->Scale(1./SOverSSBHist->GetMaximum() * scaleval);
496  SOverSSBHist->SetLineColor(kViolet+1);
497  SOverSSBHist->SetLineWidth(2);
498  SOverSSBHist->Draw("same hist");
499 
500  SOverSBHist->Scale(0.8/SOverSBHist->GetMaximum() * scaleval);
501  SOverSBHist->SetLineColor(kOrange+1);
502  SOverSBHist->SetLineWidth(2);
503  SOverSBHist->Draw("same hist");
504 
505  legFOM->SetFillStyle(0);
506  legFOM->SetLineWidth(0);
507  legFOM->SetTextFont(43);
508  legFOM->SetTextSize(textSize*0.75);
509  legFOM->AddEntry(EffTimesPurHist, "Eff #times Purity", "l");
510  legFOM->AddEntry(SOverSBHist , "S/#sqrt{B}", "l");
511  legFOM->AddEntry(SOverSSBHist , "S/#sqrt{S+B+#sigma^{2}_{syst}}", "l");
512 
513  gPad->RedrawAxis();
514 
515  //draw an axis on the right side
516  TGaxis *axis = new TGaxis(EffTimesPurHist->GetBinLowEdge(EffTimesPurHist->GetNbinsX()+1), 0.,
517  EffTimesPurHist->GetBinLowEdge(EffTimesPurHist->GetNbinsX()+1), histSim->GetMaximum(),
518  0.0001,maxVal,510,"+L");
519 
520  axis->SetTitleFont(43);
521  axis->SetTitleSize(textSize);
522  axis->SetLabelFont(43);
523  axis->SetLabelSize(textSize);
524  axis->SetTitle("Figure Of Merit");
525  axis->SetTitleOffset(1.4);
526  axis->CenterTitle();
527  axis->Draw();
528  }
529  // draw boxes
530  TBox* box = new TBox();
531  TBox* box2 = new TBox(-5000, -5000, -5000, -5000);
532  box->SetY1(0);
533  box->SetY2(histSim->GetMaximum());
534  if (histCutSide.at(iH) == kLow){
535  box->SetX1(histSim->GetBinLowEdge(1));
536  box->SetX2(currentCuts.at(iH).at(0));
537  }
538  if (histCutSide.at(iH) == kHigh){
539  box->SetX1(currentCuts.at(iH).at(0));
540  box->SetX2(histSim->GetBinLowEdge(histSim->GetNbinsX()+1));
541  }
542  if (histCutSide.at(iH) == kBoth){
543  box->SetX1(histSim->GetBinLowEdge(1));
544  box->SetX2(currentCuts.at(iH).at(0));
545  box2->SetX1(currentCuts.at(iH).at(1));
546  box2->SetX2(histSim->GetBinLowEdge(histSim->GetNbinsX()+1));
547  box2->SetLineWidth(0);
548  box2->SetFillColorAlpha(kGray, 0.4);
549  box2->SetY1(0);
550  box2->SetY2(histSim->GetMaximum());
551  box2->Draw("same");
552  }
553  box->SetLineWidth(0);
554  box->SetFillColorAlpha(kGray, 0.4);
555  if (currentCuts.at(iH).at(0) != 0 || currentCuts.at(iH).size() > 1)
556  box->Draw("same");
557 
558  //gPad->Update();
559  gPad->RedrawAxis();
560  gPad->Modified();
561 
562  if (plotFOM.at(iH) == true)
563  legFOM->Draw("same");
564 
565  TLegend* leg = new TLegend(0.2, 0.70, 0.5, 0.89);
566  leg->AddEntry(histSim , "All Events" , "l" );
567  leg->AddEntry(histSimTrueNC, "True NC Events", "f" );
568  if (currentCuts.at(iH).at(0) != 0 || currentCuts.at(iH).size() > 1)
569  leg->AddEntry(box , "Proposed Cuts" , "lf");
570  if (isDrawData)
571  leg->AddEntry(histDatNoCut, "Data", "p");
572  leg->SetLineWidth(0);
573  leg->SetFillStyle(0);
574  leg->SetTextFont(43);
575  leg->SetTextSize(textSize*0.75);
576  leg->Draw("same");
577 
578  if (!isDrawData){
579  c1->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+".png").c_str());
580  c1->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+".pdf").c_str());
581  }
582  else{
583  c1Ratio->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+".png").c_str());
584  c1Ratio->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+".pdf").c_str());
585  }
586  thisPad->cd();
587 
588  // get out of detector plots
589  TH1D* histSimTrueVtxOutOfDet = (TH1D*)_file0->Get(("sim"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)+"TrueVtxOutOfDet" ).c_str());
590  if (histSimTrueVtxOutOfDet == nullptr) continue;
591  histSimTrueVtxOutOfDet -> Sumw2();
592  histSimTrueVtxOutOfDet -> Scale(scalef);
593  histSimTrueVtxOutOfDet -> SetLineColor(kRed+1);
594  histSimTrueVtxOutOfDet -> SetFillColor(kRed+1);
595  histSimTrueVtxOutOfDet -> SetLineWidth(2);
596  histSim -> Draw("hist");
597  histSimTrueVtxOutOfDet -> Draw("same hist");
598  histSim -> Draw("same hist");
599 
600  if (currentCuts.at(iH).at(0) != 0 || currentCuts.at(iH).size() > 1)
601  box->Draw("same");
602  if (histCutSide.at(iH) == kBoth)
603  box2->Draw("same");
604 
605  //gPad->Update();
606  gPad->RedrawAxis();
607  gPad->Modified();
608 
609  TLegend *leg2 = new TLegend(0.2, 0.7, 0.5, 0.89);
610  leg2->AddEntry(histSim, "All Events", "l");
611  leg2->AddEntry(histSimTrueVtxOutOfDet, "Rock Events", "f");
612  leg2->AddEntry(box , "Proposed Cuts" , "lf");
613  leg2->SetLineWidth(0);
614  leg2->SetFillStyle(0);
615  leg2->SetTextFont(43);
616  leg2->SetTextSize(textSize*0.75);
617 
618  leg2->Draw("same");
619 
620  c1->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+"containment.png").c_str());
621  c1->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+"containment.pdf").c_str());
622  }
623  }
624 
625  for (int iH = 0; iH < histNames2D.size(); iH++){
626 
627  TCanvas* c1 = new TCanvas("c1", "c1", 600, 400);
628  TPad *thisPad = new TPad("topPad", "", 0.005, 0.005, 0.995, 0.995);
629  thisPad ->SetTopMargin(0.1);
630  thisPad ->SetBottomMargin(0.15);
631  thisPad->Draw();
632 
633 
634  c1->SetRightMargin(0.18);
635  c1->cd();
636  thisPad->cd();
637 
638  for (int icl = 0; icl < cutLevels.size(); ++icl){
639  TH2D* histSim = (TH2D*)_file0->Get(("sim"+histNames2D.at(iH)+cutLevels.at(icl)).c_str());
640  TH2D* histSimTrueNC = (TH2D*)_file0->Get(("sim"+histNames2D.at(iH)+cutLevels.at(icl)+"TrueNC").c_str());
641 
642  histSim -> SetTitle(("Background"+histTitles2D.at(iH)).c_str());
643  histSimTrueNC -> SetTitle(("True NC in Detector"+histTitles2D.at(iH)).c_str());
644 
645  gStyle -> SetTitleFont(43);
646  gStyle -> SetTitleFontSize(0.06);
647 
648  histSim ->GetXaxis()->SetTitleFont(43);
649  histSim ->GetXaxis()->SetTitleSize(textSize);
650  histSim ->GetXaxis()->CenterTitle();
651  histSim ->GetYaxis()->SetTitleFont(43);
652  histSim ->GetYaxis()->SetTitleSize(textSize);
653  histSim ->GetYaxis()->CenterTitle();
654  histSimTrueNC ->GetXaxis()->SetTitleFont(43);
655  histSimTrueNC ->GetXaxis()->SetTitleSize(textSize);
656  histSimTrueNC ->GetXaxis()->CenterTitle();
657  histSimTrueNC ->GetYaxis()->SetTitleFont(43);
658  histSimTrueNC ->GetYaxis()->SetTitleSize(textSize);
659  histSimTrueNC ->GetYaxis()->CenterTitle();
660 
661  gPad->Modified();
662  gPad->Update();
663 
664  histSim->Add(histSimTrueNC, -1);
665  histSim->Draw("colz");
666  c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+".png").c_str());
667  c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+".pdf").c_str());
668 
669  histSimTrueNC->Draw("colz");
670  c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+"TrueNC.png").c_str());
671  c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+"TrueNC.pdf").c_str());
672 
673  histSim->Draw("box");
674  histSimTrueNC->SetLineColor(kAzure+1);
675  histSimTrueNC->Draw("box same");
676 
677  TLegend* leg = new TLegend(0.6, 0.75, 0.75, 0.85);
678  leg->AddEntry(histSim, "All Events");
679  leg->AddEntry(histSimTrueNC, "True NC Events");
680  leg->SetLineWidth(0);
681  leg->SetFillStyle(0);
682  leg->SetTextFont(43);
683  leg->SetTextSize(textSize*0.75);
684  leg->Draw("same");
685 
686  c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+"together.png").c_str());
687  c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+"together.pdf").c_str());
688  }
689  }
690 }
void PlotPreSelectionPlots()
const XML_Char * name
Definition: expat.h:151
tree Draw("slc.nhit")
enum BeamMode kOrange
bin1_2sigma SetFillColor(3)
enum BeamMode kRed
TH1D * SOverSSBHist
Definition: PlotSpectra.h:65
T sqrt(T number)
Definition: d0nt_math.hpp:156
double calcSOverSSB(double sig, double bg, double fracUn)
constexpr T pow(T x)
Definition: pow.h:75
TPad * topPad
Definition: PlotSpectra.h:74
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
double calcSOverSB(double sig, double bg)
void getMaximumFOM(TH1D *h)
TH1D * EffTimesPurHist
total data
Definition: PlotSpectra.h:63
hmean SetLineWidth(2)
hmean SetLineColor(4)
TPad * thisPad
Definition: PlotSpectra.h:78
double calcPurity(double selSig, double selTot)
OStream cout
Definition: OStream.cxx:6
TH1D * getFOMHist(std::string fomName, CutSide thisCutSide, FomType thisFomType, TH1D *all, TH1D *sig, TH1D *trueSig, TH2D *energyCorr=nullptr, std::vector< float > fracUncertainties={})
TH1D * SOverSBHist
Definition: PlotSpectra.h:64
double getFracUncertainty(TH1D *h, std::vector< float > fracUncertainties)
void printValues(std::vector< double > theseValues, std::vector< double > theseFOMValues)
enum BeamMode kViolet
Float_t proj
Definition: plot.C:35
c1
Definition: demo5.py:24
double calcEfficiency(double selSig, double totSig)
FomType
Definition: FOMUtilities.h:9
enum BeamMode kGreen
simulatedPeak Scale(1/simulationNormalisationFactor)
double getFOM(FomType thisFomType, double selSig, double selBg, double totSig, double selTot, double fracUn)
enum BeamMode string