nue_pid_effs.C
Go to the documentation of this file.
3 
4 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Vars/Vars.h"
11 
12 #include "TArrow.h"
13 #include "TCanvas.h"
14 #include "TH1.h"
15 #include "TLegend.h"
16 
17 #include <cmath>
18 #include <fstream>
19 #include <iomanip>
20 #include <iostream>
21 
22 #include "Utilities/func/MathUtil.h"
23 
24 using namespace ana;
25 
26 const double kPOT = 2.82e20;
27 
28 // 250MeV bins
30 const Binning kModeBinning = Binning::Simple(4, -0.5, 3.5);
31 const Binning kYBinning = Binning::Simple(50, 0, 1);
32 
33 /***************************************************************************/
34 
36 {
37  name = "plots/"+prefix+"_"+name;
38  if(!suffix.empty()) name += "_"+suffix;
39 
40  gPad->Print((name+".eps").c_str());
41  gPad->Print((name+".pdf").c_str());
42  gPad->Print((name+".png").c_str());
43 }
44 
45 // True E, *y for NCs
46 const Var kTrueEVis([](const caf::SRProxy* sr)
47  {
48  if(sr->mc.nu.empty()) return 0.;
49  double ret = sr->mc.nu[0].E;
50  if(!sr->mc.nu[0].iscc) ret *= sr->mc.nu[0].y;
51  return ret;
52  });
53 
54 const Var kMode([](const caf::SRProxy* sr)
55  {
56  if(sr->mc.nu.empty()) return 0;
57  return sr->mc.nu[0].mode;
58  });
59 
60 
61 const Var kY([](const caf::SRProxy* sr)
62  {
63  if(sr->mc.nu.empty()) return 0.f;
64  return sr->mc.nu[0].y;
65  });
66 
67 /***************************************************************************/
68 // Nominal osc params
70 {
71  calc->SetL(810);
72  calc->SetRho(2.75);
73  calc->SetDmsq21(7.6e-5);
74  calc->SetDmsq32(2.35e-3);
75  calc->SetTh12(asin(sqrt(.87))/2); // PDG
76  calc->SetTh13(asin(sqrt(.095))/2);
77  calc->SetTh23(TMath::Pi()/4);
78  calc->SetdCP(0);
79 }
80 
81 /***************************************************************************/
82 
83 // Draw cuts on plots
84 void Arrows(TH1* sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow = false)
85 {
86  sigId->GetYaxis()->SetRangeUser(0,1.05*sigId->GetMaximum());
87 
88  double maxy = sigId->GetMaximum()/2;
89  if(noarrow)
90  maxy = maxy*2;
91 
92  TLine* lin = new TLine(bestCutSoB, 0, bestCutSoB, maxy);
93  lin->SetLineWidth(3);
94  lin->SetLineColor(kGreen+2);
95  lin->Draw();
96  if(!noarrow){
97  TArrow* arr = new TArrow(bestCutSoB-0.003, maxy, bestCutSoB+.03, maxy, .02, "|>");
98  arr->SetLineWidth(3);
99  arr->SetLineColor(kGreen+2);
100  arr->SetFillColor(kGreen+2);
101  arr->Draw();
102  maxy *= 1.5;
103  }
104 
105 
106 
107  lin = new TLine(bestCutSoSB, 0, bestCutSoSB, maxy);
108  lin->SetLineWidth(3);
109  lin->SetLineStyle(2);
110  lin->SetLineColor(kGreen+3);
111  lin->Draw();
112 
113  if(!noarrow){
114  TArrow* arr = new TArrow(bestCutSoSB-0.003, maxy, bestCutSoSB+.04, maxy, .02, "|>");
115  arr->SetLineWidth(3);
116  arr->SetLineColor(kGreen+3);
117  arr->SetFillColor(kGreen+3);
118  arr->Draw();
119  }
120 
121 
122  lin = new TLine(faCut, 0, faCut, maxy);
123  lin->SetLineWidth(3);
124  lin->SetLineColor(kRed);
125  lin->Draw();
126 
127  if(!noarrow){
128  TArrow* arr = new TArrow(faCut, maxy, faCut+.04, maxy, .02, "|>");
129  arr->SetLineWidth(3);
130  arr->SetLineColor(kRed);
131  arr->SetFillColor(kRed);
132  arr->Draw();
133  }
134 }
135 
136 /***************************************************************************/
137 // Finds best cuts
138 void FOMPlot(TH1* sig, TH1* bkg, TH1* bkgCC, TH1* bkgBeam,
139  TH1* sig_denom, double* bestCuts, double faCut,
140  std::string suffix, std::string eff_suffix)
141 {
142  std::cout << suffix << " " << eff_suffix << ":" << std::endl;
143  std::cout << " &\t Cut &\t $s/\\sqrt{b}$ &\t $s/\\sqrt{s+b}$ &\t $\\nu_e$ sig &\t Tot bkg &\t NC &\t $\\nu_\\mu$ CC &\t Beam $\\nu_e$ &\t Sig eff &\t Purity\\\\"
144  << std::endl;
145 
146  TH1F* eff = new TH1F(("eff_"+suffix+"_"+eff_suffix).c_str(),
147  "Efficiency vs Cut Value; ;Efficiency (%)",
148  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
149  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
150  TH1F* pur = new TH1F(("pur_"+suffix+"_"+eff_suffix).c_str(),
151  "Purity vs Cut Value; ;Purity (%)",
152  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
153  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
154  TH1F* fomsob = new TH1F(("sob_"+suffix+"_"+eff_suffix).c_str(),
155  "s/#sqrt{b} vs Cut Value; ;s/#sqrt{b}",
156  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
157  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
158  TH1F* fomsosb = new TH1F(("sosb_"+suffix+"_"+eff_suffix).c_str(),
159  "s/#sqrt{s+b} vs Cut Value; ;s/#sqrt{s+b}",
160  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
161  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
162 
163  eff->GetXaxis()->CenterTitle();
164  eff->GetYaxis()->CenterTitle();
165  pur->GetXaxis()->CenterTitle();
166  pur->GetYaxis()->CenterTitle();
167  fomsob->GetXaxis()->CenterTitle();
168  fomsob->GetYaxis()->CenterTitle();
169  fomsosb->GetXaxis()->CenterTitle();
170  fomsosb->GetYaxis()->CenterTitle();
171 
172  eff->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
173  pur->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
174 
175  fomsob->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
176  fomsosb->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
177 
178  for(int metric = 0; metric < 3; ++metric){
179  double bestSoB = -1;
180  double bestSoSB = -1;
181  double bestSig = 0;
182  double bestBkg = 0;
183  double bestBkgCC = 0;
184  double bestBeam = 0;
185  double totsig = sig_denom->Integral(0, -1);
186  double fomeff, fompur;
187 
188  for(int n = 0; n < sig->GetNbinsX()+2; ++n){
189  const double nsig = sig->Integral(n, -1);
190  const double nbkg = bkg->Integral(n, -1);
191  const double nbkgcc = bkgCC->Integral(n, -1);
192  const double nbkgbeam = bkgBeam->Integral(n, -1);
193 
194  if(nbkg == 0) break;
195 
196  double neff = nsig*100/(totsig);
197  double npur = nsig*100/(nsig+nbkg);
198 
199  const double fomSoB = nsig/sqrt(nbkg);
200  const double fomSoSB = nsig/sqrt(nsig+nbkg);
201 
202  if(metric == 0){
203  // only fill these hists once.
204  pur->Fill(sig->GetXaxis()->GetBinCenter(n), npur);
205  eff->Fill(sig->GetXaxis()->GetBinCenter(n), neff);
206  fomsob->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoB);
207  fomsosb->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoSB);
208  }
209 
210  if((metric == 0 && fomSoB > bestSoB) ||
211  (metric == 1 && fomSoSB > bestSoSB) ||
212  (metric == 2 && fabs(sig->GetXaxis()->GetBinLowEdge(n) - faCut) < .01)){
213  bestSoB = fomSoB;
214  bestSoSB = fomSoSB;
215  bestSig = nsig;
216  bestBkg = nbkg;
217  bestBkgCC = nbkgcc;
218  bestBeam = nbkgbeam;
219  bestCuts[metric] = sig->GetXaxis()->GetBinLowEdge(n);
220  fomeff = neff;
221  fompur = npur;
222  }
223  }
224 
225  if(metric == 0) std::cout << "$s/\\sqrt{b}$ opt";
226  if(metric == 1) std::cout << "$s/\\sqrt{s+b}$ opt";
227  if(metric == 2) std::cout << "FA cut";
228 
229  std::cout << std::fixed << std::setprecision(2);
230  std::cout << " &\t " << ((metric == 2) ? faCut : bestCuts[metric])
231  << " &\t " << bestSoB << " &\t " << bestSoSB
232  << " &\t " << bestSig << " &\t " << bestBkg
233  << " &\t " << bestBkg-bestBkgCC-bestBeam
234  << " &\t " << bestBkgCC << " &\t " << bestBeam
235  << " &\t " << fomeff << "\\% &\t " << fompur << "\\%\\\\"<< std::endl;
236 
237  } // end for metric
238 
239  new TCanvas;
240  eff->Draw("l");
241  Arrows(eff, bestCuts[0], bestCuts[1], faCut, true);
242  Print("eff", suffix+"_"+eff_suffix);
243 
244  new TCanvas;
245  pur->Draw("l");
246  Arrows(pur, bestCuts[0], bestCuts[1], faCut, true);
247  Print("pur", suffix);
248 
249  new TCanvas;
250  fomsob->Draw("l");
251  Arrows(fomsob, bestCuts[0], bestCuts[1], faCut, true);
252  Print("fomsob", suffix);
253 
254  new TCanvas;
255  fomsosb->Draw("l");
256  Arrows(fomsosb, bestCuts[0], bestCuts[1], faCut, true);
257  Print("fomsosb", suffix);
258 }
259 
260 /***************************************************************************/
261 
263  TH1*& hSig, TH1*& hNC, TH1*& hCC, TH1*& hBeam, TH1*&hTotBkg)
264 {
265  hSig = pred->PredictComponent(calc,
267  Current::kCC,
269  hSig->SetLineColor(kNueSignalColor);
270 
271  hNC = pred->PredictComponent(calc,
273  Current::kNC,
275  hNC->SetLineColor(kNCBackgroundColor);
276 
277  hCC = pred->PredictComponent(calc,
279  Current::kCC,
281  hCC->SetLineColor(kNumuBackgroundColor);
282 
283  hBeam = pred->PredictComponent(calc,
285  Current::kCC,
287  hBeam->SetLineColor(kBeamNueBackgroundColor);
288 
289  hTotBkg = (TH1*)hNC->Clone(UniqueName().c_str());
290  hTotBkg->Add(hCC);
291  hTotBkg->Add(hBeam);
292 
293  hTotBkg->SetLineColor(kTotalMCColor);
294 }
295 
296 /***************************************************************************/
297 
299 {
300  TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
301  GetSpectra(pred, calc, hSig, hNC, hCC, hBeam, hTotBkg);
302 
303  new TCanvas;
304  hSig->SetTitle(title.c_str());
305  hSig->GetYaxis()->SetTitle(TString::Format("Events / %.02lf#times10^{20} POT", kPOT/1e20));
306  hSig->GetXaxis()->CenterTitle();
307  hSig->GetYaxis()->CenterTitle();
308  hSig->Draw("hist");
309  if(maxy) hSig->GetYaxis()->SetRangeUser(0, maxy);
310  hNC->Draw("hist same");
311  hCC->Draw("hist same");
312  hBeam->Draw("hist same");
313 
314  static bool once = true;
315  if(once){
316  once = false;
317  TVirtualPad* bak = gPad;
318  new TCanvas;
319  TLegend* leg = new TLegend(.1, .1, .9, .9);
320  leg->AddEntry(hSig, "#nu_{e} signal", "l");
321  leg->AddEntry(hNC, "NC background", "l");
322  leg->AddEntry(hCC, "#nu_{#mu} background", "l");
323  leg->AddEntry(hBeam, "Beam #nu_{e} background", "l");
324  leg->Draw();
325  gPad->Print("plots/eff_leg.pdf");
326  gPad->Print("plots/eff_leg.eps");
327  gPad->Print("plots/eff_leg.png");
328  bak->cd();
329  }
330 }
331 
332 void PrintEffs(TH1* h, const char* title)
333 {
334  std::cout << title;
335  for(int i = 0; i < h->GetNbinsX()+2; ++i){
336  if(h->GetBinCenter(i) > .75 &&
337  h->GetBinCenter(i) < 3.5){
338  std::cout << " &\t " << h->GetBinContent(i);
339  }
340  }
341  std::cout << "\\\\" << std::endl;
342 }
343 
344 /***************************************************************************/
345 
346 void PlotEffs(IPrediction* predNum, IPrediction* predDenom,
348  std::string label, std::string id, double maxy = 0)
349 {
350  TH1 *hSigNum, *hNCNum, *hCCNum, *hBeamNum, *hTotBkgNum;
351  GetSpectra(predNum, calc, hSigNum, hNCNum, hCCNum, hBeamNum, hTotBkgNum);
352 
353  TH1 *hSigDenom, *hNCDenom, *hCCDenom, *hBeamDenom, *hTotBkgDenom;
354  GetSpectra(predDenom, calc, hSigDenom, hNCDenom, hCCDenom, hBeamDenom, hTotBkgDenom);
355 
356  hSigDenom->SetLineStyle(2);
357  hNCDenom->SetLineStyle(2);
358  hBeamDenom->SetLineStyle(2);
359  hTotBkgDenom->SetLineStyle(2);
360  hCCDenom->SetLineStyle(2);
361 
362  new TCanvas();
363 
364  if( label == "mode" ){
365  hNCDenom->GetXaxis()->SetLabelSize(0.06);
366  hNCDenom->GetXaxis()->SetBinLabel(1, "QE");
367  hNCDenom->GetXaxis()->SetBinLabel(2, "Res");
368  hNCDenom->GetXaxis()->SetBinLabel(3, "DIS");
369  hNCDenom->GetXaxis()->SetBinLabel(4, "Coh");
370 
371  hSigNum->GetXaxis()->SetLabelSize(0.06);
372  hSigNum->GetXaxis()->SetBinLabel(1, "QE");
373  hSigNum->GetXaxis()->SetBinLabel(2, "Res");
374  hSigNum->GetXaxis()->SetBinLabel(3, "DIS");
375  hSigNum->GetXaxis()->SetBinLabel(4, "Coh");
376 
377  }
378 
379  hNCDenom->Draw("hist");
380  hCCDenom->Draw("hist same");
381  hBeamDenom->Draw("hist same");
382  hSigDenom->Draw("hist same");
383  Print(label, id, "nocut");
384 
385  new TCanvas();
386  hSigNum->Draw("hist");
387  hNCNum->Draw("hist same");
388  hBeamNum->Draw("hist same");
389  hCCNum->Draw("hist same");
390  Print(label, id, "cut");
391 
392  new TCanvas;
393 
394  std::cout << "Total efficiency"
395  << " sig: " << 100*hSigNum->Integral(0, -1)/hSigDenom->Integral(0, -1)
396  << " NC: " << 100*hNCNum->Integral(0, -1)/hNCDenom->Integral(0, -1)
397  << " CC: " << 100*hCCNum->Integral(0, -1)/hCCDenom->Integral(0, -1)
398  << " Beam: " << 100*hBeamNum->Integral(0, -1)/hBeamDenom->Integral(0, -1) << std::endl;
399 
400  hSigNum->Divide(hSigDenom);
401  hSigNum->Scale(100);
402  hSigNum->SetTitle(title.c_str());
403  hSigNum->GetYaxis()->SetTitle("Efficiency (%)");
404  if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
405  hSigNum->GetXaxis()->CenterTitle();
406  hSigNum->GetYaxis()->CenterTitle();
407  hSigNum->Draw("hist");
408 
409  PrintEffs(hSigNum, "$\\nu_e$ sig");
410 
411  hNCNum->Divide(hNCDenom);
412  hNCNum->Scale(100);
413  hNCNum->Draw("hist same");
414 
415  PrintEffs(hNCNum, "NC");
416 
417  hCCNum->Divide(hCCDenom);
418  hCCNum->Scale(100);
419  hCCNum->Draw("hist same");
420 
421  PrintEffs(hCCNum, "$\\nu_\\mu$ CC");
422 
423  hBeamNum->Divide(hBeamDenom);
424  hBeamNum->Scale(100);
425  hBeamNum->Draw("hist same");
426 
427  PrintEffs(hBeamNum, "beam $\\nu_e$");
428 }
429 
430 /***************************************************************************/
431 
432 std::string PIDLongName(int pidIdx, bool rhc)
433 {
435  if(pidIdx == -1) ret = "No selection";
436  if(pidIdx == 0) ret = "LEM";
437  if(pidIdx == 1) ret = "LID";
438  if(pidIdx == 2) ret = "RVP";
439 
440  if(rhc) ret += " (RHC)"; else ret += " (FHC)";
441 
442  return ret;
443 }
444 
445 /***************************************************************************/
446 
447 std::string PIDFileTag(int pidIdx, bool rhc)
448 {
450  if(pidIdx == -1) ret = "nosel";
451  if(pidIdx == 0) ret = "lem";
452  if(pidIdx == 1) ret = "lid";
453  if(pidIdx == 2) ret = "rvp";
454 
455  if(rhc) ret += "_rhc"; else ret += "_fhc";
456 
457  return ret;
458 }
459 
460 /***************************************************************************/
461 
463 {
464 
466 
467  for(int rhc = false; rhc <= true; ++rhc){
468  if(rhc) continue; // Don't exist yet for S13-06-18
469 
470  const std::string fnameMC =
471  TString::Format("defname: prod_caf_S15-05-22_fd_genie_fhc_nonswap_ideal with stride 1",
472  rhc ? "rhc" : "fhc").Data();
473 
474  const std::string fnameSwap =
475  TString::Format("defname: prod_caf_S15-05-22_fd_genie_fhc_fluxswap_ideal with stride 1",
476  rhc ? "rhc" : "fhc").Data();
477 
478  SpectrumLoader loader(fnameMC);
479  SpectrumLoader loaderSwap(fnameSwap);
480 
481  const int kNumPIDs = 2;
482 
483  IPrediction* preds[kNumPIDs];
484 
485  preds[0] = new PredictionNoExtrap(loader, loaderSwap,
486  "LEM", Binning::Simple(110, 0.0, 1.1),
488  preds[1] = new PredictionNoExtrap(loader, loaderSwap,
489  "LID", Binning::Simple(110, 0.0, 1.1),
490  kElecID, kNueFirstAnaFullLIDPresel);
491 
492  PredictionNoExtrap predDenomNoCutRecoE(loader, loaderSwap,
493  "Reco E (GeV)", kRecoEBinning,
494  kCaloE, kNoCut);
495 
496  IPrediction* predDenomPreselRecoE[kNumPIDs];
497  predDenomPreselRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
498  "Reco E (GeV)",
499  kRecoEBinning,
500  kCaloE,
502  predDenomPreselRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
503  "Reco E (GeV)",
504  kRecoEBinning,
505  kCaloE,
507 
508  PredictionNoExtrap predDenomNoCutTrueE(loader, loaderSwap,
509  "True Visible E (GeV)",
510  kRecoEBinning,
511  kTrueEVis, kNoCut);
512 
513  IPrediction* predDenomPreselTrueE[kNumPIDs];
514  predDenomPreselTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
515  "True Visible E (GeV)",
516  kRecoEBinning,
517  kTrueEVis,
519  predDenomPreselTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
520  "True Visible E (GeV)",
521  kRecoEBinning,
522  kTrueEVis,
524 
525  IPrediction* predsRecoE[kNumPIDs];
526  IPrediction* predsTrueE[kNumPIDs];
527 
528  PredictionNoExtrap predDenomMode(loader, loaderSwap,
529  "", kModeBinning,
530  kMode, kNoCut);
531 
532  IPrediction* predsMode[kNumPIDs];
533 
534  PredictionNoExtrap predDenomY(loader, loaderSwap,
535  "", kYBinning,
536  kY, kNoCut);
537 
538  IPrediction* predsY[kNumPIDs];
539 
540  predsRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
541  "Reco E (GeV)", kRecoEBinning,
543 
544  predsRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
545  "Reco E (GeV)", kRecoEBinning,
547 
548  predsTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
549  "True Visible E (GeV)", kRecoEBinning,
551 
552  predsTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
553  "True Visible E (GeV)", kRecoEBinning,
555 
556  predsMode[0] = new PredictionNoExtrap(loader, loaderSwap,
557  "", kModeBinning,
559 
560  predsMode[1] = new PredictionNoExtrap(loader, loaderSwap,
561  "", kModeBinning,
563 ;
564  predsY[0] = new PredictionNoExtrap(loader, loaderSwap,
565  "True y", kYBinning,
567 
568  predsY[1] = new PredictionNoExtrap(loader, loaderSwap,
569  "True y", kYBinning,
571 
572  loader.Go();
573  loaderSwap.Go();
574 
575  for(int pidIdx = 0; pidIdx < kNumPIDs; ++pidIdx){
576  TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
577  GetSpectra(preds[pidIdx], &calc, hSig, hNC, hCC, hBeam, hTotBkg);
578 
579  double cuts[2];
580  const double faCut = (pidIdx == 0) ? 0.8 : 0.95;
581  TH1* sigdenomNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
583  Current::kCC,
585  sigdenomNoCut->SetLineColor(kNueSignalColor);
586 
587  TH1* sigdenomPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
589  Current::kCC,
591  sigdenomPresel->SetLineColor(kNueSignalColor);
592 
593  FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomNoCut, cuts, faCut, PIDFileTag(pidIdx, rhc), "nocut");
594 
595  FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomPresel, cuts, faCut, PIDFileTag(pidIdx, rhc), "presel");
596 
597  PlotSpectra(preds[pidIdx], &calc, PIDLongName(pidIdx, rhc));
598 
599  Arrows(hSig, cuts[0], cuts[1], faCut);
600 
601  Print("pid", PIDFileTag(pidIdx, rhc));
602  } // end for pidIdx
603 
604  PlotSpectra(&predDenomNoCutRecoE, &calc, PIDLongName(-1, rhc), 40);
605  Print("recoE", PIDFileTag(-1, rhc));
606  PlotSpectra(&predDenomNoCutTrueE, &calc, PIDLongName(-1, rhc), 40);
607  Print("trueE", PIDFileTag(-1, rhc));
608 
609  for(int pidIdx = 0; pidIdx < kNumPIDs; ++pidIdx){
610  std::cout << std::endl << std::endl << PIDLongName(pidIdx, rhc) << ":" << std::endl;
611 
612  std::cout << "Reco E:" << std::endl;
613 
614  PlotSpectra(predsRecoE[pidIdx], &calc,
615  PIDLongName(pidIdx, rhc), rhc ? 6 : 2.5);
616  Print("recoE", PIDFileTag(pidIdx, rhc));
617 
618  PlotEffs(predsRecoE[pidIdx], &predDenomNoCutRecoE, &calc,
619  PIDLongName(pidIdx, rhc), "recoE", PIDFileTag(pidIdx, rhc));
620  Print("recoE", PIDFileTag(pidIdx, rhc), "eff_nocut");
621 
622  PlotEffs(predsRecoE[pidIdx], predDenomPreselRecoE[pidIdx], &calc,
623  PIDLongName(pidIdx, rhc), "recoE", PIDFileTag(pidIdx, rhc));
624  Print("recoE", PIDFileTag(pidIdx, rhc), "eff_presel");
625 
626  PlotEffs(predsMode[pidIdx], &predDenomMode, &calc,
627  PIDLongName(pidIdx, rhc), "mode", PIDFileTag(pidIdx, rhc));
628  Print("mode", PIDFileTag(pidIdx, rhc), "eff");
629 
630  PlotEffs(predsY[pidIdx], &predDenomY, &calc,
631  PIDLongName(pidIdx, rhc), "y", PIDFileTag(pidIdx, rhc));
632  Print("y", PIDFileTag(pidIdx, rhc), "eff");
633 
634 
635 
636  std::cout << std::endl << "True E vis:" << std::endl;
637  PlotSpectra(predsTrueE[pidIdx], &calc,
638  PIDLongName(pidIdx, rhc), rhc ? 6 : 2.5);
639  Print("trueE", PIDFileTag(pidIdx, rhc));
640 
641  PlotEffs(predsTrueE[pidIdx], &predDenomNoCutTrueE, &calc,
642  PIDLongName(pidIdx, rhc), "trueE", PIDFileTag(pidIdx, rhc));
643  Print("trueE", PIDFileTag(pidIdx, rhc), "eff_nocut");
644 
645  PlotEffs(predsTrueE[pidIdx], predDenomPreselTrueE[pidIdx], &calc,
646  PIDLongName(pidIdx, rhc), "trueE", PIDFileTag(pidIdx, rhc));
647  Print("trueE", PIDFileTag(pidIdx, rhc), "eff_presel");
648  }
649  } // end for rhc
650 }
const Binning kModeBinning
const Cut kNueFirstAnaLIDOpt([](const caf::SRProxy *sr){return(sr->sel.lid.ann > 0.95);})
From Tian&#39;s docdb 12673 optimized lid cut.
const XML_Char * name
Definition: expat.h:151
const Var kMode([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?-1:int(sr->mc.nu[0].mode);})
Neutrino interaction mode.
Definition: Vars.h:106
void PlotEffs(IPrediction *predNum, IPrediction *predDenom, osc::IOscCalculator *calc, std::string title, std::string label, std::string id, double maxy=0)
Definition: nue_pid_effs.C:346
const Var kLEM
PID
Definition: Vars.cxx:24
virtual void SetL(double L)=0
osc::OscCalculatorDumb calc
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Oscillation analysis framework, runs over CAF files outside of ART.
const Binning kYBinning
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
(&#39; appearance&#39;)
Definition: IPrediction.h:18
void FOMPlot(TH1 *sig, TH1 *bkg, TH1 *bkgCC, TH1 *bkgBeam, TH1 *sig_denom, double *bestCuts, double faCut, std::string suffix, std::string eff_suffix)
Definition: nue_pid_effs.C:138
VectorProxy< SRNeutrinoProxy > nu
Definition: SRProxy.h:1890
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
double maxy
string fnameSwap
Definition: demo4.py:6
(&#39;beam &#39;)
Definition: IPrediction.h:15
virtual void SetRho(double rho)=0
void PrintEffs(TH1 *h, const char *title)
Definition: nue_pid_effs.C:332
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Var kTrueEVis([](const caf::SRProxy *sr){if(sr->mc.nu.empty()) return 0.;double ret=sr->mc.nu[0].E;if(!sr->mc.nu[0].iscc) ret *=sr->mc.nu[0].y;return ret;})
virtual void SetdCP(const T &dCP)=0
void Arrows(TH1 *sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow=false)
Definition: nue_pid_effs.C:84
Proxy for StandardRecord.
Definition: SRProxy.h:2237
virtual void SetDmsq21(const T &dmsq21)=0
void nue_pid_effs()
Definition: nue_pid_effs.C:462
const Color_t kNumuBackgroundColor
Definition: Style.h:31
const Color_t kNueSignalColor
Definition: Style.h:19
const Cut kNueFirstAnaLEMOpt([](const caf::SRProxy *sr){return(sr->sel.lem.pid > 0.80);})
From Tian&#39;s docdb 12673 optimized lem cut.
const char * label
const Cut kNueFirstAnaFullLIDPresel
From Tian&#39;s docdb 12673.
Charged-current interactions.
Definition: IPrediction.h:39
General interface to any calculator that lets you set the parameters.
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:50
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
void resetCalc(osc::IOscCalculatorAdjustable *calc)
Definition: nue_pid_effs.C:69
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
virtual void Go() override
Load all the registered spectra.
virtual void SetTh23(const T &th23)=0
virtual void SetDmsq32(const T &dmsq32)=0
loader
Definition: demo0.py:10
const XML_Char * prefix
Definition: expat.h:380
std::string PIDFileTag(int pidIdx, bool rhc)
Definition: nue_pid_effs.C:447
OStream cout
Definition: OStream.cxx:6
const Var kY([](const caf::SRProxy *sr){if(sr->mc.nu.empty()) return 0.f;return sr->mc.nu[0].y;})
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
SRTruthBranchProxy mc
Definition: SRProxy.h:2253
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
virtual void SetTh13(const T &th13)=0
void Print(std::string prefix, std::string name, std::string suffix="")
Definition: nue_pid_effs.C:35
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Color_t kTotalMCColor
Definition: Style.h:17
Prediction that just uses FD MC, with no extrapolation.
std::string PIDLongName(int pidIdx, bool rhc)
Definition: nue_pid_effs.C:432
const Binning kRecoEBinning
All neutrinos, any flavor.
Definition: IPrediction.h:26
void PlotSpectra(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string title, std::string det, int POT, std::string option)
Float_t e
Definition: plot.C:35
loaderSwap
Definition: demo4.py:9
const Color_t kNCBackgroundColor
Definition: Style.h:22
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
void GetSpectra(IPrediction *pred, osc::IOscCalculator *calc, TH1 *&hSig, TH1 *&hNC, TH1 *&hCC, TH1 *&hBeam, TH1 *&hTotBkg)
Definition: nue_pid_effs.C:262
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:117
static constexpr Double_t sr
Definition: Munits.h:164
h
Definition: demo3.py:41
T asin(T number)
Definition: d0nt_math.hpp:60
virtual void SetTh12(const T &th12)=0
const Cut kNueFirstAnaFullLEMPresel
From Tian&#39;s docdb 12673.