nue_pid_effs_paper_numu_neweff.C
Go to the documentation of this file.
1 #ifdef __CINT__
3 {
4  std::cout << "Sorry, you must run in compiled mode" << std::endl;
5 }
6 #else
7 
9 #include "OscLib/OscCalcDumb.h"
10 
11 #include "CAFAna/Cuts/Cuts.h"
12 #include "CAFAna/Analysis/Style.h"
13 #include "CAFAna/Cuts/NueCutsFirstAna.h"
14 #include "CAFAna/Cuts/NumuCuts.h"
17 #include "CAFAna/Core/Utilities.h"
18 #include "CAFAna/Vars/Vars.h"
20 
21 #include "TArrow.h"
22 #include "TCanvas.h"
23 #include "TH1.h"
24 #include "TH2.h"
25 #include "TGraph.h"
26 #include "TLegend.h"
27 
28 #include <cmath>
29 #include <fstream>
30 #include <iomanip>
31 #include <iostream>
32 
33 #include "Utilities/func/MathUtil.h"
34 
35 using namespace ana;
36 const double kPOTMACRO = 6.754e20;
37 
38 // 250MeV bins
40 const Binning kModeBinning = Binning::Simple(4, -0.5, 3.5);
41 const Binning kYBinning = Binning::Simple(50, 0, 1);
42 
43 /***************************************************************************/
44 
46 {
47  name = "plots/assessPaperNuMu/"+prefix+"_"+name;
48  if(!suffix.empty()) name += "_"+suffix;
49 
50  gPad->Print((name+".eps").c_str());
51  gPad->Print((name+".pdf").c_str());
52  gPad->Print((name+".png").c_str());
53 }
54 
55 const Cut kNumuFAFDContain(
56  //{ "mc.nu.iscc", "mc.nu.pdg", "mc.nu.E",
57  // "sel.contain.missE", "mc.nnu", "energy.numusimp.trkccE","sel.contain.numucontain","slc.nhit","sel.remid.pid","trk.ncosmic"},
58  [](const caf::SRProxy* sr)
59  {
60  if(sr->mc.nnu == 0) return false;
61  assert(sr->mc.nnu == 1);
62  return (
63  (
64  (sr->mc.nu[0].E < 5 && (sr->mc.nu[0].iscc==1 && abs(sr->mc.nu[0].pdg)==14)) ||
65  (sr->mc.nu[0].iscc==0) ||
66  (abs(sr->mc.nu[0].pdg)!=14)
67  ) &&
68  (
69  ((sr->sel.contain.missE/sr->mc.nu[0].E) < 0.01 && (sr->mc.nu[0].iscc==1 && abs(sr->mc.nu[0].pdg)==14)) ||
70  (sr->mc.nu[0].iscc==0) ||
71  (abs(sr->mc.nu[0].pdg)!=14)
72  ) &&
73  sr->slc.nhit>=20 &&
74  sr->sel.remid.pid>0 &&
75  sr->trk.cosmic.ntracks>0 &&
76  sr->mc.nnu>0 &&
77  sr->energy.numu.trkccE < 5 &&
79  );
80  }
81  );
82 
83 
85  //{ "mc.nu.iscc", "mc.nu.pdg", "mc.nu.E",
86  // "sel.contain.missE", "mc.nnu", "energy.numusimp.trkccE","sel.contain.numucontain","slc.nhit","sel.remid.pid","trk.ncosmic"},
87  [](const caf::SRProxy* sr)
88  {
89  if(sr->mc.nnu == 0) return false;
90  assert(sr->mc.nnu == 1);
91  return (
92  (
93  (sr->mc.nu[0].E < 5 && (sr->mc.nu[0].iscc==1 && abs(sr->mc.nu[0].pdg)==14)) ||
94  (sr->mc.nu[0].iscc==0) ||
95  (abs(sr->mc.nu[0].pdg)!=14)
96  ) &&
97  (
98  ((sr->sel.contain.missE/sr->mc.nu[0].E) < 0.01 && (sr->mc.nu[0].iscc==1 && abs(sr->mc.nu[0].pdg)==14)) ||
99  (sr->mc.nu[0].iscc==0) ||
100  (abs(sr->mc.nu[0].pdg)!=14)
101  ) &&
102  sr->mc.nnu>0 &&
103  sr->energy.numu.trkccE < 5
104  );
105  }
106  );
107 
108 // True E, *y for NCs
109 const Var kTrueEVis(//{"mc.nu.E", "mc.nu.y", "mc.nu.iscc"},
110  [](const caf::SRProxy* sr)
111  {
112  if(sr->mc.nu.empty()) return 0.;
113  double ret = sr->mc.nu[0].E;
114  if(!sr->mc.nu[0].iscc) ret *= sr->mc.nu[0].y;
115  return ret;
116  });
117 
118 /*
119 const Var kMode(//{"mc.nu.mode"},
120  [](const caf::SRProxy* sr)
121  {
122  int tmp = 0;
123  if(sr->mc.nu.empty()) tmp;
124  tmp = sr->mc.nu[0].mode;
125  return tmp;
126  });
127 
128 */
129 
130 const Var kY(//{"mc.nu.y"},
131  [](const caf::SRProxy* sr)
132  {
133  float tmp = 0.f;
134  if(sr->mc.nu.empty()) return tmp;
135  tmp = sr->mc.nu[0].y;
136  return tmp;
137  });
138 
139 /***************************************************************************/
140 // Nominal osc params
142 {
143  calc->SetL(810);
144  calc->SetRho(2.75);
145  calc->SetDmsq21(7.6e-5);
146  calc->SetDmsq32(2.35e-3);
147  calc->SetTh12(asin(sqrt(.87))/2); // PDG
148  calc->SetTh13(asin(sqrt(.095))/2);
149  calc->SetTh23(TMath::Pi()/4);
150  calc->SetdCP(0);
151 }
152 
153 /***************************************************************************/
154 
155 // Draw cuts on plots
156 void Arrows(TH1* sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow = false)
157 {
158  sigId->GetYaxis()->SetRangeUser(0,1.05*sigId->GetMaximum());
159 
160  double maxy = sigId->GetMaximum()/2;
161  if(noarrow)
162  maxy = maxy*2;
163 
164  TLine* lin = new TLine(bestCutSoB, 0, bestCutSoB, maxy);
165  lin->SetLineWidth(3);
166  lin->SetLineColor(kGreen+2);
167  lin->Draw();
168  if(!noarrow){
169  TArrow* arr = new TArrow(bestCutSoB-0.003, maxy, bestCutSoB+.03, maxy, .02, "|>");
170  arr->SetLineWidth(3);
171  arr->SetLineColor(kGreen+2);
172  arr->SetFillColor(kGreen+2);
173  arr->Draw();
174  maxy *= 1.5;
175  }
176 
177 
178 
179  lin = new TLine(bestCutSoSB, 0, bestCutSoSB, maxy);
180  lin->SetLineWidth(3);
181  lin->SetLineStyle(2);
182  lin->SetLineColor(kGreen+3);
183  lin->Draw();
184 
185  if(!noarrow){
186  TArrow* arr = new TArrow(bestCutSoSB-0.003, maxy, bestCutSoSB+.04, maxy, .02, "|>");
187  arr->SetLineWidth(3);
188  arr->SetLineColor(kGreen+3);
189  arr->SetFillColor(kGreen+3);
190  arr->Draw();
191  }
192 
193 
194  lin = new TLine(faCut, 0, faCut, maxy);
195  lin->SetLineWidth(3);
196  lin->SetLineColor(kRed);
197  //lin->Draw();
198 
199  if(!noarrow){
200  TArrow* arr = new TArrow(faCut, maxy, faCut+.04, maxy, .02, "|>");
201  arr->SetLineWidth(3);
202  arr->SetLineColor(kRed);
203  arr->SetFillColor(kRed);
204  //arr->Draw();
205  }
206 }
207 
208 /***************************************************************************/
209 // Finds best cuts
210 void FOMPlot(TH1* sig, TH1* bkg, TH1* bkgCC, TH1* bkgBeam,
211  TH1* sig_denom, double* bestCuts, double faCut,
212  std::string suffix, std::string eff_suffix)
213 {
214  std::cout << suffix << " " << eff_suffix << ":" << std::endl;
215  std::cout << " &\t Cut &\t $s/\\sqrt{b}$ &\t $s/\\sqrt{s+b}$ &\t $\\nu_\\mu$ sig &\t Tot bkg &\t NC &\t Appeared $\\nu_e$ &\t Beam $\\nu_e$ &\t Sig eff &\t Purity\\\\"
216  << std::endl;
217 
218  TH1F* eff = new TH1F(("eff_"+suffix+"_"+eff_suffix).c_str(),
219  "Efficiency vs Cut Value; ;Efficiency (%)",
220  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
221  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
222 
223  TH1F* effxpur = new TH1F(("effxpur_"+suffix+"_"+eff_suffix).c_str(),
224  "Efficiency*Purity vs Cut Value; ;Efficiency*Purity (%)",
225  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
226  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
227 
228  TH1F* pur = new TH1F(("pur_"+suffix+"_"+eff_suffix).c_str(),
229  "Purity vs Cut Value; ;Purity (%)",
230  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
231  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
232  TH1F* fomsob = new TH1F(("sob_"+suffix+"_"+eff_suffix).c_str(),
233  "s/#sqrt{b} vs Cut Value; ;s/#sqrt{b}",
234  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
235  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
236  TH1F* fomsosb = new TH1F(("sosb_"+suffix+"_"+eff_suffix).c_str(),
237  "s/#sqrt{s+b} vs Cut Value; ;s/#sqrt{s+b}",
238  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
239  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
240 
241  eff->GetXaxis()->CenterTitle();
242  eff->GetYaxis()->CenterTitle();
243  pur->GetXaxis()->CenterTitle();
244  pur->GetYaxis()->CenterTitle();
245  fomsob->GetXaxis()->CenterTitle();
246  fomsob->GetYaxis()->CenterTitle();
247  fomsosb->GetXaxis()->CenterTitle();
248  fomsosb->GetYaxis()->CenterTitle();
249 
250  eff->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
251  pur->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
252 
253  fomsob->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
254  fomsosb->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
255  int ngpoints=0;
256  for(int n = 0; n < sig->GetNbinsX()+2; ++n){
257  if(bkg->Integral(n, -1) == 0) break;
258  ngpoints++;
259  }
260  TGraph* effpurGraph=new TGraph (ngpoints);
261  for(int metric = 0; metric < 3; ++metric){
262  double bestSoB = -1;
263  double bestSoSB = -1;
264  double bestSig = 0;
265  double bestBkg = 0;
266  double bestBkgCC = 0;
267  double bestBeam = 0;
268  double totsig = sig_denom->Integral(0, -1);
269  double fomeff, fompur;
270 
271  for(int n = 0; n < sig->GetNbinsX()+2; ++n){
272  const double nsig = sig->Integral(n, -1);
273  const double nbkg = bkg->Integral(n, -1);
274  const double nbkgcc = bkgCC->Integral(n, -1);
275  const double nbkgbeam = bkgBeam->Integral(n, -1);
276 
277  if(nbkg == 0) break;
278 
279  double neff = nsig*100/(totsig);
280  double npur = nsig*100/(nsig+nbkg);
281 
282  const double fomSoB = nsig/sqrt(nbkg);
283  const double fomSoSB = nsig/sqrt(nsig+nbkg);
284 
285  if(metric == 0){
286  // only fill these hists once.
287  pur->Fill(sig->GetXaxis()->GetBinCenter(n), npur);
288  eff->Fill(sig->GetXaxis()->GetBinCenter(n), neff);
289  effxpur->Fill(sig->GetXaxis()->GetBinCenter(n), 100*(neff/100)*(npur/100));
290  effpurGraph->SetPoint(n,neff,npur);
291  fomsob->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoB);
292  fomsosb->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoSB);
293  }
294 
295  if((metric == 0 && fomSoB > bestSoB) ||
296  (metric == 1 && fomSoSB > bestSoSB) ||
297  (metric == 2 && fabs(sig->GetXaxis()->GetBinLowEdge(n) - faCut) < .01)){
298  bestSoB = fomSoB;
299  bestSoSB = fomSoSB;
300  bestSig = nsig;
301  bestBkg = nbkg;
302  bestBkgCC = nbkgcc;
303  bestBeam = nbkgbeam;
304  bestCuts[metric] = sig->GetXaxis()->GetBinLowEdge(n);
305  fomeff = neff;
306  fompur = npur;
307  }
308  }
309 
310  if(metric == 0) std::cout << "$s/\\sqrt{b}$ opt";
311  if(metric == 1) std::cout << "$s/\\sqrt{s+b}$ opt";
312  if(metric == 2) std::cout << "FA cut";
313 
314  std::cout << std::fixed << std::setprecision(2);
315  std::cout << " &\t " << ((metric == 2) ? faCut : bestCuts[metric])
316  << " &\t " << bestSoB << " &\t " << bestSoSB
317  << " &\t " << bestSig << " &\t " << bestBkg
318  << " &\t " << bestBkg-bestBkgCC-bestBeam
319  << " &\t " << bestBkgCC << " &\t " << bestBeam
320  << " &\t " << fomeff << "\\% &\t " << fompur << "\\%\\\\"<< std::endl;
321 
322  } // end for metric
323 
324  new TCanvas;
325  eff->Draw("l");
326  Arrows(eff, bestCuts[0], bestCuts[1], faCut, true);
327  Print("eff", suffix+"_"+eff_suffix);
328 
329  new TCanvas;
330  pur->Draw("l");
331  Arrows(pur, bestCuts[0], bestCuts[1], faCut, true);
332  Print("pur", suffix);
333 
334  new TCanvas;
335  eff->SetTitle("");
336  eff->GetYaxis()->SetTitle("Percentage");
337  eff->SetLineColor(kRed);
338  pur->SetLineColor(kBlue);
339  effxpur->SetLineColor(kGreen+2);
340  eff->SetMaximum(110);
341 
342  eff->Draw("c");
343  pur->Draw("csame");
344  effxpur->Draw("csame");
345  //Arrows(pur, bestCuts[0], bestCuts[1], faCut, true);
346 
347  TLegend* legeffpur = new TLegend(.125, .12, .5, .32);
348  legeffpur->AddEntry(eff, "Cumulative Efficiency", "l");
349  legeffpur->AddEntry(pur, "Cumulative Purity", "l");
350  legeffpur->AddEntry(effxpur, "Cumulative Efficiency #times Purity", "l");
351  legeffpur->SetFillStyle(0);
352  legeffpur->Draw();
353  Print("effpur", suffix+"_"+eff_suffix);
354 
355  new TCanvas;
356  TH2F* hpx = new TH2F("hpx","hpx",100,0,100,100,0,100); // axis range hpx->SetStats(kFALSE); // no statistics
357  hpx->SetStats(kFALSE);
358  hpx->SetTitle("");
359  hpx->GetXaxis()->SetTitle("Efficiency (%)");
360  hpx->GetYaxis()->SetTitle("Purity (%)");
361  hpx->GetXaxis()->CenterTitle();
362  hpx->GetYaxis()->CenterTitle();
363  hpx->Draw();
364  effpurGraph->SetLineWidth(2);
365  effpurGraph->Draw("C");
366  Print("effpurgraph", suffix+"_"+eff_suffix);
367 
368 
369  new TCanvas;
370  fomsob->Draw("l");
371  Arrows(fomsob, bestCuts[0], bestCuts[1], faCut, true);
372  Print("fomsob", suffix);
373 
374  new TCanvas;
375  fomsosb->Draw("l");
376  Arrows(fomsosb, bestCuts[0], bestCuts[1], faCut, true);
377  Print("fomsosb", suffix);
378 }
379 
380 /***************************************************************************/
381 
383  TH1*& hSig, TH1*& hNC, TH1*& hCC, TH1*& hBeam, TH1*&hTotBkg)
384 {
385  hSig = pred->PredictComponent(calc,
387  Current::kCC,
389  hSig->SetLineColor(kNumuBackgroundColor);
390 
391  hNC = pred->PredictComponent(calc,
393  Current::kNC,
395  hNC->SetLineColor(kNCBackgroundColor);
396 
397  hCC = pred->PredictComponent(calc,
399  Current::kCC,
401  hCC->SetLineColor(kNueSignalColor);
402 
403  hBeam = pred->PredictComponent(calc,
405  Current::kCC,
407  hBeam->SetLineColor(kBeamNueBackgroundColor);
408 
409  hTotBkg = (TH1*)hNC->Clone(UniqueName().c_str());
410  hTotBkg->Add(hCC);
411  hTotBkg->Add(hBeam);
412 
413  hTotBkg->SetLineColor(kTotalMCColor);
414 }
415 
416 /***************************************************************************/
417 
419 {
420  TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
421  GetSpectra(pred, calc, hSig, hNC, hCC, hBeam, hTotBkg);
422 
423  new TCanvas;
424  hSig->SetTitle("");//title.c_str());
425  hSig->GetYaxis()->SetTitle("Events / 18#times10^{20} POT");
426  //hSig->GetYaxis()->SetTitle("Arbitary Units");
427  hSig->GetXaxis()->CenterTitle();
428  hSig->GetYaxis()->CenterTitle();
429  hSig->Draw("hist");
430  if(maxy) hSig->GetYaxis()->SetRangeUser(0, maxy);
431  hNC->Draw("hist same");
432  hCC->Draw("hist same");
433  hBeam->Draw("hist same");
434  TLegend* leg2 = new TLegend(.125, .55, .5, .85);
435  leg2->AddEntry(hCC, "Appeared #nu_{e}", "l");
436  leg2->AddEntry(hSig, "Survived #nu_{#mu}", "l");
437  leg2->AddEntry(hNC, "NC background", "l");
438  leg2->AddEntry(hBeam, "Beam #nu_{e} background", "l");
439  leg2->SetFillStyle(0);
440  leg2->Draw();
441 
442  static bool once = true;
443  if(once){
444  once = false;
445  TVirtualPad* bak = gPad;
446  new TCanvas;
447  TLegend* leg = new TLegend(.1, .1, .9, .9);
448  leg->AddEntry(hSig, "#nu_{#mu} signal", "l");
449  leg->AddEntry(hNC, "NC background", "l");
450  leg->AddEntry(hCC, "Appeared #nu_{e} background", "l");
451  leg->AddEntry(hBeam, "Beam #nu_{e} background", "l");
452  leg->Draw();
453  gPad->Print("plots/assessPaperNuMu/eff_leg.pdf");
454  gPad->Print("plots/assessPaperNuMu/eff_leg.eps");
455  gPad->Print("plots/assessPaperNuMu/eff_leg.png");
456  bak->cd();
457  }
458 }
459 
460 void PrintEffs(TH1* h, const char* title)
461 {
462  std::cout << title;
463  for(int i = 0; i < h->GetNbinsX()+2; ++i){
464  if(h->GetBinCenter(i) > .75 &&
465  h->GetBinCenter(i) < 3.5){
466  std::cout << " &\t " << h->GetBinContent(i);
467  }
468  }
469  std::cout << "\\\\" << std::endl;
470 }
471 
472 /***************************************************************************/
473 
474 void PlotEffs(IPrediction* predNum, IPrediction* predDenom,
476  std::string label, std::string id, double maxy = 0)
477 {
478  TH1 *hSigNum, *hNCNum, *hCCNum, *hBeamNum, *hTotBkgNum;
479  GetSpectra(predNum, calc, hSigNum, hNCNum, hCCNum, hBeamNum, hTotBkgNum);
480 
481  TH1 *hSigDenom, *hNCDenom, *hCCDenom, *hBeamDenom, *hTotBkgDenom;
482  GetSpectra(predDenom, calc, hSigDenom, hNCDenom, hCCDenom, hBeamDenom, hTotBkgDenom);
483 
484  hSigDenom->SetLineStyle(2);
485  hNCDenom->SetLineStyle(2);
486  hBeamDenom->SetLineStyle(2);
487  hTotBkgDenom->SetLineStyle(2);
488  hCCDenom->SetLineStyle(2);
489 
490  new TCanvas();
491 
492  if( label == "mode" ){
493  hNCDenom->GetXaxis()->SetLabelSize(0.06);
494  hNCDenom->GetXaxis()->SetBinLabel(1, "QE");
495  hNCDenom->GetXaxis()->SetBinLabel(2, "Res");
496  hNCDenom->GetXaxis()->SetBinLabel(3, "DIS");
497  hNCDenom->GetXaxis()->SetBinLabel(4, "Coh");
498 
499  hSigNum->GetXaxis()->SetLabelSize(0.06);
500  hSigNum->GetXaxis()->SetBinLabel(1, "QE");
501  hSigNum->GetXaxis()->SetBinLabel(2, "Res");
502  hSigNum->GetXaxis()->SetBinLabel(3, "DIS");
503  hSigNum->GetXaxis()->SetBinLabel(4, "Coh");
504 
505  }
506 
507  hNCDenom->Draw("hist");
508  hCCDenom->Draw("hist same");
509  hBeamDenom->Draw("hist same");
510  hSigDenom->Draw("hist same");
511  Print(label, id, "nocut");
512 
513  new TCanvas();
514  if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
515  //hSigNum->GetYaxis()->SetTitle("Arbitary Units");
516  hSigNum->GetYaxis()->CenterTitle();
517  hSigNum->Draw("hist");
518  hNCNum->Draw("hist same");
519  hBeamNum->Draw("hist same");
520  hCCNum->Draw("hist same");
521  Print(label, id, "cut");
522 
523  new TCanvas;
524 
525  std::cout << "Total efficiency"
526  << " sig: " << 100*hSigNum->Integral(0, -1)/hSigDenom->Integral(0, -1)
527  << " NC: " << 100*hNCNum->Integral(0, -1)/hNCDenom->Integral(0, -1)
528  << " CC: " << 100*hCCNum->Integral(0, -1)/hCCDenom->Integral(0, -1)
529  << " Beam: " << 100*hBeamNum->Integral(0, -1)/hBeamDenom->Integral(0, -1) << std::endl;
530 
531  hSigNum->Divide(hSigDenom);
532  hSigNum->Scale(100);
533  hSigNum->SetTitle(title.c_str());
534  hSigNum->GetYaxis()->SetTitle("Efficiency (%)");
535  hSigNum->GetYaxis()->SetRangeUser(0, 100);
536  //if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
537  hSigNum->GetXaxis()->CenterTitle();
538  hSigNum->GetYaxis()->CenterTitle();
539  hSigNum->Draw("hist");
540 
541  PrintEffs(hSigNum, "$\\nu_e$ sig");
542 
543  hNCNum->Divide(hNCDenom);
544  hNCNum->Scale(100);
545  hNCNum->Draw("hist same");
546 
547  PrintEffs(hNCNum, "NC");
548 
549  hCCNum->Divide(hCCDenom);
550  hCCNum->Scale(100);
551  hCCNum->Draw("hist same");
552 
553  PrintEffs(hCCNum, "$\\nu_\\mu$ CC");
554 
555  hBeamNum->Divide(hBeamDenom);
556  hBeamNum->Scale(100);
557  hBeamNum->Draw("hist same");
558 
559  PrintEffs(hBeamNum, "beam $\\nu_e$");
560 }
561 
562 /***************************************************************************/
563 
564 std::string PIDLongName(int pidIdx, bool rhc)
565 {
567  if(pidIdx == -1) ret = "No selection";
568  if(pidIdx == 0) ret = "remID";
569  if(pidIdx == 1) ret = "CVN#mu";
570 
571  if(rhc) ret += " (RHC)"; else ret += " (FHC)";
572 
573  return ret;
574 }
575 
576 /***************************************************************************/
577 
578 std::string PIDFileTag(int pidIdx, bool rhc)
579 {
581  if(pidIdx == -1) ret = "nosel";
582  if(pidIdx == 0) ret = "lem";
583  if(pidIdx == 1) ret = "cvn";
584 
585  if(rhc) ret += "_rhc"; else ret += "_fhc";
586 
587  return ret;
588 }
589 
590 /***************************************************************************/
591 
593 {
594 
596  calc.SetL(810);
597  calc.SetRho(2.84);
598  calc.SetDmsq21(7.53e-5);
599  calc.SetDmsq32(2.44e-3);
600  calc.SetTh12(asin(sqrt(0.846))/2);
601  calc.SetTh13(asin(sqrt(.085))/2);
602  calc.SetTh23(asin(sqrt(0.999))/2);
603  calc.SetdCP(0);
604 
605  for(int rhc = false; rhc <= true; ++rhc){
606  if(rhc) continue; // Don't exist yet for S13-06-18
607 
608  /*
609  const std::string fnameMC =
610  TString::Format("defname: prod_caf_S15-05-22_fd_genie_fhc_nonswap_ideal with stride 1",
611  rhc ? "rhc" : "fhc").Data();
612 
613  const std::string fnameSwap =
614  TString::Format("defname: prod_caf_S15-05-22_fd_genie_fhc_fluxswap_ideal with stride 1",
615  rhc ? "rhc" : "fhc").Data();
616  */
617 
618  //SpectrumLoader loader("defname: prod_caf_R17-03-01-prod3reco.g_fd_genie_nonswap_fhc_nova_v08_full_v1 with limit 3");
619  //SpectrumLoader loaderSwap("defname: prod_caf_R17-03-01-prod3reco.g_fd_genie_fluxswap_fhc_nova_v08_full_v1 with limit 3");
620 
621  //SpectrumLoader loader("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/fardet_genie_nonswap_genierw_fhc_v08_1000_r00014457_s20_c000_S17-05-31_v1_20170330_113817_sim.caf.root");
622  //SpectrumLoader loaderSwap("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/fardet_genie_fluxswap_genierw_fhc_v08_1000_r00014031_s02_c000_S17-05-31_v1_20170330_145648_sim.caf.root");
623 
624  SpectrumLoader loader("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/semfe_nonS_big.root");
625  SpectrumLoader loaderSwap("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/semfe_S_big.root");
626 
627  //SpectrumLoader loader("/pnfs/nova/scratch/users/atsaris/cvnFS_eval_fullInt_afterColl_swap/00014*/*/*.root");
628  //SpectrumLoader loaderSwap("/pnfs/nova/scratch/users/atsaris/cvnFS_eval_fullInt_afterColl_nonSwap/00014*/*/*.root");
629 
630  const int kNumPIDs = 2;
631 
632  IPrediction* preds[kNumPIDs];
633  //kNumuFAFDContain
634  preds[0] = new PredictionNoExtrap(loader, loaderSwap,
635  "remID", Binning::Simple(100, 0.0, 1.0),
637 
638  preds[1] = new PredictionNoExtrap(loader, loaderSwap,
639  "#nu_{#mu} CC Classifier Output", Binning::Simple(100, 0.0, 1.0),
640  //kANGNuMu, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain);
642 
643  PredictionNoExtrap predDenomNoCutRecoE(loader, loaderSwap,
644  "Reco E (GeV)", kRecoEBinning,
645  kCaloE, kNumuFAFDContainLight);//kNumuQuality && kNumuContainFD);
646 
647  IPrediction* predDenomPreselRecoE[kNumPIDs];
648  predDenomPreselRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
649  "Reco E (GeV)",
650  kRecoEBinning,
651  kCaloE,
653 
654  predDenomPreselRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
655  "Reco E (GeV)",
656  kRecoEBinning,
657  kCaloE,
659 
660  PredictionNoExtrap predDenomNoCutTrueE(loader, loaderSwap,
661  "True Visible E (GeV)",
662  kRecoEBinning,
663  kTrueEVis, kNoCut);
664 
665  IPrediction* predDenomPreselTrueE[kNumPIDs];
666  predDenomPreselTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
667  "True Visible E (GeV)",
668  kRecoEBinning,
669  kTrueEVis,
671  predDenomPreselTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
672  "True Visible E (GeV)",
673  kRecoEBinning,
674  kTrueEVis,
676 
677  IPrediction* predsRecoE[kNumPIDs];
678  IPrediction* predsTrueE[kNumPIDs];
679 
680  PredictionNoExtrap predDenomMode(loader, loaderSwap,
681  "", kModeBinning,
682  kMode, kNoCut);
683 
684  IPrediction* predsMode[kNumPIDs];
685 
686  PredictionNoExtrap predDenomY(loader, loaderSwap,
687  "", kYBinning,
688  kY, kNoCut);
689 
690  IPrediction* predsY[kNumPIDs];
691 
692  predsRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
693  "Reco E (GeV)", kRecoEBinning,
694  kCaloE, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaLEMOpt);
695 
696  predsRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
697  "Reco E (GeV)", kRecoEBinning,
698  //kCaloE, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaANGOpt);
700 
701  predsTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
702  "True Visible E (GeV)", kRecoEBinning,
703  kTrueEVis, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaLEMOpt);
704 
705  predsTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
706  "True Visible E (GeV)", kRecoEBinning,
707  //kTrueEVis, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaANGOpt);
709 
710  predsMode[0] = new PredictionNoExtrap(loader, loaderSwap,
711  "", kModeBinning,
712  kMode, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaLEMOpt);
713 
714  predsMode[1] = new PredictionNoExtrap(loader, loaderSwap,
715  "", kModeBinning,
716  //kMode, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaANGOpt);
718 
719  predsY[0] = new PredictionNoExtrap(loader, loaderSwap,
720  "True y", kYBinning,
721  kY, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaLEMOpt);
722 
723  predsY[1] = new PredictionNoExtrap(loader, loaderSwap,
724  "True y", kYBinning,
725  //kY, kNumuQuality && kNumuContainFD && kNumuCosmicRej && kNumuFAFDContain && kNueFirstAnaANGOpt);
727 
728  loader.Go();
729  loaderSwap.Go();
730 
731  for(int pidIdx = 0; pidIdx < kNumPIDs; ++pidIdx){
732  TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
733  GetSpectra(preds[pidIdx], &calc, hSig, hNC, hCC, hBeam, hTotBkg);
734 
735  double cuts[2];
736  const double faCut = 0.75;
737 
738  TH1* backoscnueNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
740  Current::kCC,
742 
743  TH1* sigdenomNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
745  Current::kCC,
747 
748  TH1* backnueNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
750  Current::kCC,
752 
753  TH1* backncNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
755  Current::kNC,
757 
758 
759  double totalBckOscNuECC=backoscnueNoCut->Integral(0, -1);
760  double totalSig=sigdenomNoCut->Integral(0, -1);
761  double totalBckNuECC=backnueNoCut->Integral(0, -1);
762  double totalBckNC=backncNoCut->Integral(0, -1);
763  double totalBackground=totalBckOscNuECC+totalBckNuECC+totalBckNC;
764  double totalsob=totalSig/sqrt(totalBackground);
765  double totalsosb=totalSig/sqrt(totalSig+totalBackground);
766  double totalefficiency=100;
767  double totalpurity=(totalSig*100)/(totalSig+totalBackground);
768 
769  std::cout << std::fixed << std::setprecision(2);
770  std::cout << " &\t " << totalsob << " &\t " << totalsosb
771  << " &\t " << totalSig << " &\t " << totalBackground
772  << " &\t " << totalBckNC
773  << " &\t " << totalBckOscNuECC << " &\t " << totalBckNuECC
774  << " &\t " << totalefficiency << "\\% &\t " << totalpurity << "\\%\\\\"<< std::endl;
775  // & 2.47 & 2.14 & 18.31 & 54.93 & 18.31 & 18.31 18.31 & 100.00\% & 75.00\%\\
776 
777  sigdenomNoCut->SetLineColor(kNueSignalColor);
778 
779  TH1* sigdenomPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
781  Current::kCC,
783  sigdenomPresel->SetLineColor(kNueSignalColor);
784 
785  FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomNoCut, cuts, faCut, PIDFileTag(pidIdx, rhc), "nocut");
786  //FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomPresel, cuts, faCut, PIDFileTag(pidIdx, rhc), "presel");
787 
788  PlotSpectra(preds[pidIdx], &calc, PIDLongName(pidIdx, rhc));//,16.0);
789 
790  Arrows(hSig, cuts[0], cuts[1], faCut);
791  Print("pid", PIDFileTag(pidIdx, rhc));
792  } // end for pidIdx
793 
794  PlotSpectra(&predDenomNoCutRecoE, &calc, PIDLongName(-1, rhc), 40);
795  Print("recoE", PIDFileTag(-1, rhc));
796  PlotSpectra(&predDenomNoCutTrueE, &calc, PIDLongName(-1, rhc), 40);
797  Print("trueE", PIDFileTag(-1, rhc));
798 
799  return;
800  } // end for rhc
801 }
802 
803 #endif
const Binning kModeBinning
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
void PrintEffs(TH1 *h, const char *title)
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:99
void Print(std::string prefix, std::string name, std::string suffix="")
caf::Proxy< float > trkccE
Definition: SRProxy.h:195
virtual void SetL(double L)=0
const double kPOTMACRO
enum BeamMode kRed
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetTh13(const T &th13) override
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
const Binning kYBinning
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void SetL(double L) override
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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:148
virtual void SetDmsq21(const T &dmsq21)=0
double maxy
(&#39;beam &#39;)
Definition: IPrediction.h:15
caf::Proxy< caf::SRContain > contain
Definition: SRProxy.h:1251
caf::Proxy< caf::SRNumuEnergy > numu
Definition: SRProxy.h:214
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Cut kNumuFAFDContainLight( [](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return(( (sr->mc.nu[0].E< 5 &&(sr->mc.nu[0].iscc==1 &&abs(sr->mc.nu[0].pdg)==14))|| (sr->mc.nu[0].iscc==0)|| (abs(sr->mc.nu[0].pdg)!=14) )&&( ((sr->sel.contain.missE/sr->mc.nu[0].E)< 0.01 &&(sr->mc.nu[0].iscc==1 &&abs(sr->mc.nu[0].pdg)==14))|| (sr->mc.nu[0].iscc==0)|| (abs(sr->mc.nu[0].pdg)!=14) )&&sr->mc.nnu >0 && sr->energy.numu.trkccE< 5 );})
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
caf::Proxy< float > pid
Definition: SRProxy.h:1136
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 SetTh13(const T &th13)=0
const Cut kNumuCosmicRej([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.535 && sr->slc.nhit< 400);})
Definition: NumuCuts.h:32
void abs(TH1 *hist)
const Cut kNumuContainFD([](const caf::SRProxy *sr){ std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return( sr->slc.ncellsfromedge > 1 &&planestofront > 1 &&planestoback > 1 &&sr->sel.contain.kalfwdcell > 10 &&sr->sel.contain.kalbakcell > 10 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 0);})
Definition: NumuCuts.h:20
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
Float_t tmp
Definition: plot.C:36
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2136
osc::OscCalcDumb calc
const Color_t kNumuBackgroundColor
Definition: Style.h:30
const Var kY([](const caf::SRProxy *sr){float tmp=0.f;if(sr->mc.nu.empty()) return tmp;tmp=sr->mc.nu[0].y;return tmp;})
virtual void SetDmsq32(const T &dmsq32)=0
const Color_t kNueSignalColor
Definition: Style.h:19
const char * label
Charged-current interactions.
Definition: IPrediction.h:39
caf::Proxy< unsigned int > nhit
Definition: SRProxy.h:1315
void nue_pid_effs_paper_numu_neweff()
const Var kRemID
PID
Definition: Vars.cxx:81
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2145
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const Cut kNumuFAFDContain( [](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return( ( (sr->mc.nu[0].E< 5 &&(sr->mc.nu[0].iscc==1 &&abs(sr->mc.nu[0].pdg)==14))|| (sr->mc.nu[0].iscc==0)|| (abs(sr->mc.nu[0].pdg)!=14) )&& ( ((sr->sel.contain.missE/sr->mc.nu[0].E)< 0.01 &&(sr->mc.nu[0].iscc==1 &&abs(sr->mc.nu[0].pdg)==14))|| (sr->mc.nu[0].iscc==0)|| (abs(sr->mc.nu[0].pdg)!=14) )&& sr->slc.nhit >=20 && sr->sel.remid.pid >0 && sr->trk.cosmic.ntracks >0 && sr->mc.nnu >0 && sr->energy.numu.trkccE< 5 && sr->sel.contain.numucontain );})
void Arrows(TH1 *sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow=false)
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
void SetTh23(const T &th23) override
caf::Proxy< caf::SRRemid > remid
Definition: SRProxy.h:1269
loader
Definition: demo0.py:10
void SetDmsq21(const T &dmsq21) override
const XML_Char * prefix
Definition: expat.h:380
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
caf::Proxy< float > missE
Definition: SRProxy.h:827
void SetDmsq32(const T &dmsq32) override
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
void FOMPlot(TH1 *sig, TH1 *bkg, TH1 *bkgCC, TH1 *bkgBeam, TH1 *sig_denom, double *bestCuts, double faCut, std::string suffix, std::string eff_suffix)
virtual void SetRho(double rho)=0
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
std::string PIDFileTag(int pidIdx, bool rhc)
void SetdCP(const T &dCP) override
virtual void SetTh23(const T &th23)=0
caf::Proxy< caf::SRTrackBase > cosmic
Definition: SRProxy.h:1795
void SetTh12(const T &th12) override
void resetCalc(osc::IOscCalcAdjustable *calc)
Neutral-current interactions.
Definition: IPrediction.h:40
assert(nhit_max >=nhit_nbins)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
std::string PIDLongName(int pidIdx, bool rhc)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const Color_t kTotalMCColor
Definition: Style.h:16
Prediction that just uses FD MC, with no extrapolation.
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
const Binning kRecoEBinning
const Cut kNumuQuality
Definition: NumuCuts.h:18
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
enum BeamMode kGreen
const Color_t kNCBackgroundColor
Definition: Style.h:22
enum BeamMode kBlue
const Var kCVNm
PID
Definition: Vars.cxx:39
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
caf::Proxy< bool > numucontain
Definition: SRProxy.h:830
void PlotEffs(IPrediction *predNum, IPrediction *predDenom, osc::IOscCalc *calc, std::string title, std::string label, std::string id, double maxy=0)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void GetSpectra(IPrediction *pred, osc::IOscCalc *calc, TH1 *&hSig, TH1 *&hNC, TH1 *&hCC, TH1 *&hBeam, TH1 *&hTotBkg)
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1730
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
void SetRho(double rho) override
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string