nue_pid_effs_miniprod.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/NueCutsSecondAna.h"
17 #include "CAFAna/Core/Utilities.h"
18 #include "CAFAna/Vars/Vars.h"
19 #include "CAFAna/Vars/NueVars.h"
20 #include "CAFAna/Cuts/NueCutsSecondAna.h"
23 #include "CAFAna/Cuts/NueCuts2017.h"
26 
27 #include "TArrow.h"
28 #include "TCanvas.h"
29 #include "TLatex.h"
30 #include "TH1.h"
31 #include "TH2.h"
32 #include "TMarker.h"
33 #include "TGraph.h"
34 
35 #include "TLegend.h"
36 
37 #include <cmath>
38 #include <fstream>
39 #include <iomanip>
40 #include <iostream>
41 
42 #include "Utilities/func/MathUtil.h"
43 
44 using namespace ana;
45 
46 string fFileName="test.txt";
47 const double kPOTMACRO = 6.754e20;
48 
49 // New CVN CAF variables
50 /// \f$ \nu_e \f$ PID
51 const Var kCVNSS([](const caf::SRProxy* sr){return sr->sel.cvnProd3Train.nueid;});
52 
53 /// \f$ \nu_e \f$ PID
54 const Var kCVNELU([](const caf::SRProxy* sr){return sr->sel.cvn.nueid;});
55 
56 /// \f$ \nu_e \f$ PID
57 const Var kCVN2017([](const caf::SRProxy* sr){return sr->sel.cvn2017.nueid;});
58 
59 // 250MeV bins
61 const Binning kModeBinning = Binning::Simple(4, -0.5, 3.5);
62 const Binning kYBinning = Binning::Simple(50, 0, 1);
63 const Binning kYvertexBinning = Binning::Simple(40, -800, 800);
64 const Binning kXvertexBinning = Binning::Simple(40, -800, 800);
66 
67 // Output to text file
68 std::ofstream output("output.tex");
69 
70 /***************************************************************************/
71 
72 // Put a "NOvA Simulation" tag in the corner
73 void Simulation()
74 {
75  TLatex* prelim = new TLatex(.9, .95, "NOvA Simulation");
76  prelim->SetTextColor(kGray+1);
77  prelim->SetNDC();
78  prelim->SetTextSize(2/30.);
79  prelim->SetTextAlign(32);
80  prelim->Draw();
81 }
82 
83 // print the current canvas ina variety of formats
85 {
86  name = "plots/SAbless/"+prefix+"_"+name;
87  if(!suffix.empty()) name += "_"+suffix;
88  Simulation();
89  gPad->Print((name+".eps").c_str());
90  gPad->Print((name+".pdf").c_str());
91  gPad->Print((name+".png").c_str());
92 }
93 
94 // define some useful variables for PID diagnostics
95 const Var kTrueEVis(//{"mc.nu.E", "mc.nu.y", "mc.nu.iscc"},
96  [](const caf::SRProxy* sr)
97  {
98  if(sr->mc.nu.empty()) return 0.;
99  double ret = sr->mc.nu[0].E;
100  if(!sr->mc.nu[0].iscc) ret *= sr->mc.nu[0].y;
101  return ret;
102  });
103 
104 // define some useful variables for PID diagnostics
105 const Cut ISMC(//{"hdr.ismc"},
106  [](const caf::SRProxy* sr)
107  {
108  if (sr->hdr.ismc) return true;
109  return false;
110  });
111 
112 const Var kY(//{"mc.nu.y"},
113  [](const caf::SRProxy* sr)
114  {
115  float tmp = 0.f;
116  if(sr->mc.nu.empty()) return tmp;
117  tmp = sr->mc.nu[0].y;
118  return tmp;
119  });
120 
121 const Var kXvertex(//{"vtx.elastic.vtx.x"},
122  [](const caf::SRProxy* sr)
123  {
124  return sr->vtx.elastic[0].vtx.x;
125  });
126 
127 const Var kYvertex(//{"vtx.elastic.vtx.y"},
128  [](const caf::SRProxy* sr)
129  {
130  return sr->vtx.elastic[0].vtx.y;
131  });
132 
133 const Var kZvertex(//{"vtx.elastic.vtx.z"},
134  [](const caf::SRProxy* sr)
135  {
136  return sr->vtx.elastic[0].vtx.z;
137  });
138 
139 
140 /***************************************************************************/
141 // Nominal osc params
143 {
144  calc->SetL(810);
145  calc->SetRho(2.75);
146  calc->SetDmsq21(7.6e-5);
147  calc->SetDmsq32(2.35e-3);
148  calc->SetTh12(asin(sqrt(.87))/2); // PDG
149  calc->SetTh13(asin(sqrt(.095))/2);
150  calc->SetTh23(TMath::Pi()/4);
151  calc->SetdCP(0);
152 }
153 
154 /***************************************************************************/
155 
156 // Draw cuts on plots
157 void Arrows(TH1* sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow = false)
158 {
159  sigId->GetYaxis()->SetRangeUser(0,1.05*sigId->GetMaximum());
160 
161  double maxy = sigId->GetMaximum()/2;
162  if(noarrow)
163  maxy = maxy*2;
164 
165  TLine* lin = new TLine(bestCutSoB, 0, bestCutSoB, maxy);
166  lin->SetLineWidth(3);
167  lin->SetLineColor(kGreen+2);
168  lin->Draw();
169  if(!noarrow){
170  TArrow* arr = new TArrow(bestCutSoB-0.003, maxy, bestCutSoB+.03, maxy, .02, "|>");
171  arr->SetLineWidth(3);
172  arr->SetLineColor(kGreen+2);
173  arr->SetFillColor(kGreen+2);
174  arr->Draw();
175  maxy *= 1.5;
176  }
177 
178 
179 
180  lin = new TLine(bestCutSoSB, 0, bestCutSoSB, maxy);
181  lin->SetLineWidth(3);
182  lin->SetLineStyle(2);
183  lin->SetLineColor(kGreen+3);
184  lin->Draw();
185 
186  if(!noarrow){
187  TArrow* arr = new TArrow(bestCutSoSB-0.003, maxy, bestCutSoSB+.04, maxy, .02, "|>");
188  arr->SetLineWidth(3);
189  arr->SetLineColor(kGreen+3);
190  arr->SetFillColor(kGreen+3);
191  arr->Draw();
192  }
193 
194 
195  lin = new TLine(faCut, 0, faCut, maxy);
196  lin->SetLineWidth(3);
197  lin->SetLineColor(kRed);
198  //lin->Draw();
199 
200  if(!noarrow){
201  TArrow* arr = new TArrow(faCut, maxy, faCut+.04, maxy, .02, "|>");
202  arr->SetLineWidth(3);
203  arr->SetLineColor(kRed);
204  arr->SetFillColor(kRed);
205  //arr->Draw();
206  }
207 }
208 
209 /***************************************************************************/
210 // Finds best cuts
211 void FOMPlot(TH1* sig, TH1* bkg, TH1* bkgCC, TH1* bkgBeam,
212  TH1* sig_denom, TH1* bck_denom,double* bestCuts, double faCut,
213  std::string suffix, std::string eff_suffix)
214 {
215  output << suffix << " " << eff_suffix << ":" << std::endl;
216  output << " &\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 Beam $\\nu_\\tau$&\t Sig eff &\t Purity\\\\"
217  << std::endl;
218 
219  TH1F* eff = new TH1F(("eff_"+suffix+"_"+eff_suffix).c_str(),
220  "Efficiency vs Cut Value; ;Efficiency (%)",
221  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
222  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
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  TH1F* pur = new TH1F(("pur_"+suffix+"_"+eff_suffix).c_str(),
228  "Purity vs Cut Value; ;Purity (%)",
229  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
230  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
231  TH1F* fomsob = new TH1F(("sob_"+suffix+"_"+eff_suffix).c_str(),
232  "s/#sqrt{b} vs Cut Value; ;s/#sqrt{b}",
233  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
234  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
235  TH1F* fomsosb = new TH1F(("sosb_"+suffix+"_"+eff_suffix).c_str(),
236  "s/#sqrt{s+b} vs Cut Value; ;s/#sqrt{s+b}",
237  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
238  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
239 
240  eff->GetXaxis()->CenterTitle();
241  eff->SetMaximum(100);
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  TGraph* rocGraph=new TGraph (ngpoints);
262  double fompur_SoB=0.;
263  double fomeff_SoB=0.;
264  double fompur_SoSB=0.;
265  double fomeff_SoSB=0.;
266  for(int metric = 0; metric < 3; ++metric){
267  double bestSoB = -1;
268  double bestSoSB = -1;
269  double bestSig = 0;
270  double bestBkg = 0;
271  double bestBkgCC = 0;
272  double bestBeam = 0;
273  double totsig = sig_denom->Integral(0, -1);
274  double totbck = bck_denom->Integral(0, -1);
275  double fomeff, fompur;
276 
277  for(int n = 0; n < sig->GetNbinsX()+2; ++n){
278  const double nsig = sig->Integral(n, -1);
279  const double nbkg = bkg->Integral(n, -1);
280  const double nbkgcc = bkgCC->Integral(n, -1);
281  const double nbkgbeam = bkgBeam->Integral(n, -1);
282 
283  if(nbkg == 0) break;
284 
285  double neff = nsig*100/(totsig);
286  double npur = nsig*100/(nsig+nbkg);
287 
288  const double fomSoB = nsig/sqrt(nbkg);
289  const double fomSoSB = nsig/sqrt(nsig+nbkg);
290 
291  if(metric == 0){
292  // only fill these hists once.
293  pur->Fill(sig->GetXaxis()->GetBinCenter(n), npur);
294  eff->Fill(sig->GetXaxis()->GetBinCenter(n), neff);
295  effxpur->Fill(sig->GetXaxis()->GetBinCenter(n), 100*(neff/100)*(npur/100));
296  effpurGraph->SetPoint(n,neff,npur);
297  rocGraph->SetPoint(n,nbkg/totbck,nsig/totsig);
298  fomsob->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoB);
299  fomsosb->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoSB);
300  }
301 
302  if((metric == 0 && fomSoB > bestSoB) ||
303  (metric == 1 && fomSoSB > bestSoSB) ||
304  (metric == 2 && fabs(sig->GetXaxis()->GetBinLowEdge(n) - faCut) < .01)){
305  bestSoB = fomSoB;
306  bestSoSB = fomSoSB;
307  bestSig = nsig;
308  bestBkg = nbkg;
309  bestBkgCC = nbkgcc;
310  bestBeam = nbkgbeam;
311  bestCuts[metric] = sig->GetXaxis()->GetBinLowEdge(n);
312  fomeff = neff;
313  fompur = npur;
314  }
315  }
316 
317  if(metric == 0) output << "$s/\\sqrt{b}$ opt";
318  if(metric == 1) output << "$s/\\sqrt{s+b}$ opt";
319  if(metric == 2) output << "FA cut";
320 
321  output << std::fixed << std::setprecision(2);
322  output << " &\t " << ((metric == 2) ? faCut : bestCuts[metric])
323  << " &\t " << bestSoB << " &\t " << bestSoSB
324  << " &\t " << bestSig << " &\t " << bestBkg
325  << " &\t " << bestBkg-bestBkgCC-bestBeam
326  << " &\t " << bestBkgCC << " &\t " << bestBeam << " &\t "
327  << " &\t " << fomeff << "\\% &\t " << fompur << "\\%\\\\"<< std::endl;
328 
329  if(metric == 0){
330  fomeff_SoB=fomeff;
331  fompur_SoB=fompur;
332  }
333  if(metric == 1){
334  fomeff_SoSB=fomeff;
335  fompur_SoSB=fompur;
336  }
337 
338  } // end for metric
339 
340  /*
341  new TCanvas;
342  eff->Draw("l");
343  Arrows(eff, bestCuts[0], bestCuts[1], faCut, true);
344  Print("eff", suffix+"_"+eff_suffix);
345 
346  new TCanvas;
347  pur->Draw("l");
348  Arrows(pur, bestCuts[0], bestCuts[1], faCut, true);
349  Print("pur", suffix);
350  */
351 
352 
353  new TCanvas;
354  eff->SetTitle("");
355  eff->GetYaxis()->SetTitle("Percentage");
356  eff->SetLineColor(kRed);
357  pur->SetLineColor(kBlue);
358  effxpur->SetLineColor(kGreen+2);
359  eff->SetMaximum(110);
360  eff->Draw("c");
361  pur->Draw("csame");
362  effxpur->Draw("csame");
363 
364  /*
365 
366  TLegend* legeffpur = new TLegend(.125, .12, .5, .32);
367  legeffpur->AddEntry(eff, "Cumulative Efficiency", "l");
368  legeffpur->AddEntry(pur, "Cumulative Purity", "l");
369  legeffpur->AddEntry(effxpur, "Cumulative Efficiency #times Purity", "l");
370  legeffpur->SetFillStyle(0);
371  legeffpur->Draw();
372  Print("effpur", suffix+"_"+eff_suffix);
373  */
374 
375  new TCanvas;
376  TH2F* hpx = new TH2F("hpx","hpx",100,0,100,100,0,100); // axis range hpx->SetStats(kFALSE); // no statistics
377  hpx->SetStats(kFALSE);
378  hpx->SetTitle("");
379  hpx->GetXaxis()->SetTitle("Efficiency (%)");
380  hpx->GetYaxis()->SetTitle("Purity (%)");
381  hpx->GetXaxis()->CenterTitle();
382  hpx->GetYaxis()->CenterTitle();
383  hpx->Draw();
384  for(double i = .5; i < 10; i += .5){
385  TGraph* c = new TGraph;
386  for(double j = .01; j < 101; j += .1){
387  c->SetPoint(c->GetN(), j, (100*i*i)/j);
388  }
389  c->SetLineWidth(2);
390  c->SetLineColor(kGray+1);
391  c->SetLineStyle(7);
392  c->Draw("l same");
393  }
394 
395 
396  TMarker* m1 = new TMarker(fomeff_SoB, fompur_SoB, kFullCircle);
397  m1->SetMarkerSize(1.5*m1->GetMarkerSize());
398  m1->Draw();
399 
400  TMarker* m2 = new TMarker(fomeff_SoSB, fompur_SoSB, kFullTriangleUp);
401  m2->SetMarkerSize(1.5*m2->GetMarkerSize());
402  m2->Draw();
403 
404  effpurGraph->SetLineWidth(2);
405  effpurGraph->Draw("C");
406  Print("effpurgraph", suffix+"_"+eff_suffix);
407 
408  new TCanvas;
409  TH2F* roch = new TH2F("roch","roch",100,0,1,100,0,1); // axis range roch->SetStats(kFALSE); // no statistics
410  roch->SetStats(kFALSE);
411  roch->SetTitle("");
412  roch->GetXaxis()->SetTitle("False Positive Rate");
413  roch->GetYaxis()->SetTitle("True Positive Rate");
414  roch->GetXaxis()->CenterTitle();
415  roch->GetYaxis()->CenterTitle();
416  roch->Draw();
417  TLine* coin = new TLine(0, 0, 1, 1);
418  coin->SetLineStyle(2);
419  coin->Draw();
420  rocGraph->SetLineWidth(2);
421  rocGraph->Draw("C");
422  Print("rocGraph", suffix+"_"+eff_suffix);
423 
424  TFile* outFile=new TFile((suffix+"_"+eff_suffix+".root").c_str(),"RECREATE");
425  rocGraph->Write("rocGraph");
426  effpurGraph->Write("effpurGraph");
427  outFile->Close();
428 
429  new TCanvas;
430  fomsob->Draw("l");
431  Arrows(fomsob, bestCuts[0], bestCuts[1], faCut, true);
432  Print("fomsob", suffix);
433 
434  new TCanvas;
435  fomsosb->Draw("l");
436  Arrows(fomsosb, bestCuts[0], bestCuts[1], faCut, true);
437  Print("fomsosb", suffix);
438 
439 }
440 
441 /***************************************************************************/
442 
444  TH1*& hSig, TH1*& hNC, TH1*& hCC, TH1*& hBeam, TH1*&hTotBkg)
445 {
446  hSig = pred->PredictComponent(calc,
448  Current::kCC,
450  hSig->SetLineColor(kNueSignalColor);
451 
452  hNC = pred->PredictComponent(calc,
454  Current::kNC,
456  hNC->SetLineColor(kNCBackgroundColor);
457 
458  hCC = pred->PredictComponent(calc,
460  Current::kCC,
462  hCC->SetLineColor(kNumuBackgroundColor);
463 
464  hBeam = pred->PredictComponent(calc,
466  Current::kCC,
468  hBeam->SetLineColor(kBeamNueBackgroundColor);
469 
470 
471  hTotBkg = (TH1*)hNC->Clone(UniqueName().c_str());
472  hTotBkg->Add(hCC);
473  hTotBkg->Add(hBeam);
474 
475  hTotBkg->SetLineColor(kTotalMCColor);
476 }
477 
478 /***************************************************************************/
479 
481 {
482  TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
483  GetSpectra(pred, calc, hSig, hNC, hCC, hBeam, hTotBkg);
484 
485  new TCanvas;
486  hSig->SetTitle("");//title.c_str());
487  double kRealPot=6e20;
488  hSig->GetYaxis()->SetTitle(TString::Format("Events / %.02lf#times10^{20} POT", kRealPot/1e20));
489  //hSig->GetYaxis()->SetTitle("Arbitary Units");
490  hSig->GetXaxis()->CenterTitle();
491  hSig->GetYaxis()->CenterTitle();
492  hSig->Draw("hist");
493  if(maxy) hSig->GetYaxis()->SetRangeUser(0, maxy);
494  hNC->Draw("hist same");
495  hCC->Draw("hist same");
496  hBeam->Draw("hist same");
497  TLegend* leg2 = new TLegend(.125, .55, .5, .85);
498  leg2->AddEntry(hSig, "Appeared #nu_{e}", "l");
499  leg2->AddEntry(hCC, "Survived #nu_{#mu}", "l");
500  leg2->AddEntry(hNC, "NC background", "l");
501  leg2->AddEntry(hBeam, "Beam #nu_{e} background", "l");
502  leg2->SetFillStyle(0);
503  leg2->Draw();
504 
505  static bool once = true;
506  if(once){
507  once = false;
508  TVirtualPad* bak = gPad;
509  new TCanvas;
510  TLegend* leg = new TLegend(.1, .1, .9, .9);
511  leg->AddEntry(hSig, "#nu_{e} signal", "l");
512  leg->AddEntry(hNC, "NC background", "l");
513  leg->AddEntry(hCC, "#nu_{#mu} background", "l");
514  leg->AddEntry(hBeam, "Beam #nu_{e} background", "l");
515  leg->Draw();
516  gPad->Print("plots/SAbless/eff_leg.pdf");
517  gPad->Print("plots/SAbless/eff_leg.eps");
518  gPad->Print("plots/SAbless/eff_leg.png");
519  bak->cd();
520  }
521 }
522 
523 void PrintEffs(TH1* h, const char* title)
524 {
525  output << title;
526  for(int i = 0; i < h->GetNbinsX()+2; ++i){
527  if(h->GetBinCenter(i) > .75 &&
528  h->GetBinCenter(i) < 3.5){
529  output << " &\t " << h->GetBinContent(i);
530  }
531  }
532  output << "\\\\" << std::endl;
533 }
534 
535 /***************************************************************************/
536 
537 void PlotEffs(IPrediction* predNum, IPrediction* predDenom,
539  std::string label, std::string id, double maxy = 0)
540 {
541  TH1 *hSigNum, *hNCNum, *hCCNum, *hBeamNum, *hTotBkgNum;
542  GetSpectra(predNum, calc, hSigNum, hNCNum, hCCNum, hBeamNum, hTotBkgNum);
543 
544  TH1 *hSigDenom, *hNCDenom, *hCCDenom, *hBeamDenom, *hTotBkgDenom;
545  GetSpectra(predDenom, calc, hSigDenom, hNCDenom, hCCDenom, hBeamDenom, hTotBkgDenom);
546 
547  hSigDenom->SetLineStyle(2);
548  hNCDenom->SetLineStyle(2);
549  hBeamDenom->SetLineStyle(2);
550  hTotBkgDenom->SetLineStyle(2);
551  hCCDenom->SetLineStyle(2);
552 
553  new TCanvas();
554 
555  if( label == "mode" ){
556  hNCDenom->GetXaxis()->SetLabelSize(0.06);
557  hNCDenom->GetXaxis()->SetBinLabel(1, "QE");
558  hNCDenom->GetXaxis()->SetBinLabel(2, "Res");
559  hNCDenom->GetXaxis()->SetBinLabel(3, "DIS");
560  hNCDenom->GetXaxis()->SetBinLabel(4, "Coh");
561 
562  hSigNum->GetXaxis()->SetLabelSize(0.06);
563  hSigNum->GetXaxis()->SetBinLabel(1, "QE");
564  hSigNum->GetXaxis()->SetBinLabel(2, "Res");
565  hSigNum->GetXaxis()->SetBinLabel(3, "DIS");
566  hSigNum->GetXaxis()->SetBinLabel(4, "Coh");
567 
568  }
569 
570  hNCDenom->Draw("hist");
571  hCCDenom->Draw("hist same");
572  hBeamDenom->Draw("hist same");
573  hSigDenom->Draw("hist same");
574  Print(label, id, "nocut");
575 
576  new TCanvas();
577  if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
578  //hSigNum->GetYaxis()->SetTitle("Arbitary Units");
579  hSigNum->GetYaxis()->CenterTitle();
580  hSigNum->Draw("hist");
581  hNCNum->Draw("hist same");
582  hBeamNum->Draw("hist same");
583  hCCNum->Draw("hist same");
584  Print(label, id, "cut");
585 
586  new TCanvas;
587 
588  output << "Total efficiency"
589  << " sig: " << 100*hSigNum->Integral(0, -1)/hSigDenom->Integral(0, -1)
590  << " NC: " << 100*hNCNum->Integral(0, -1)/hNCDenom->Integral(0, -1)
591  << " CC: " << 100*hCCNum->Integral(0, -1)/hCCDenom->Integral(0, -1)
592  << " Beam: " << 100*hBeamNum->Integral(0, -1)/hBeamDenom->Integral(0, -1) << std::endl;
593 
594  hSigNum->Divide(hSigDenom);
595  hSigNum->Scale(100);
596  hSigNum->SetTitle("");
597  hSigNum->GetYaxis()->SetTitle("Efficiency (%)");
598  hSigNum->GetYaxis()->SetRangeUser(0, 100);
599  //if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
600  hSigNum->GetXaxis()->CenterTitle();
601  hSigNum->GetYaxis()->CenterTitle();
602  hSigNum->Draw("hist");
603 
604  PrintEffs(hSigNum, "$\\nu_e$ sig");
605 
606  hNCNum->Divide(hNCDenom);
607  hNCNum->Scale(100);
608  hNCNum->Draw("hist same");
609 
610  PrintEffs(hNCNum, "NC");
611 
612  hCCNum->Divide(hCCDenom);
613  hCCNum->Scale(100);
614  hCCNum->Draw("hist same");
615 
616  PrintEffs(hCCNum, "$\\nu_\\mu$ CC");
617 
618  hBeamNum->Divide(hBeamDenom);
619  hBeamNum->Scale(100);
620  hBeamNum->Draw("hist same");
621 
622  PrintEffs(hBeamNum, "beam $\\nu_e$");
623 
624  double x1=125;
625  double y1=.6;
626  double x2=.425;
627  double y2=.85;
628 
629  if(title=="y"){
630  x1=0.6;
631  y1=.6;
632  x2=.9;
633  y2=.85;
634  }
635 
636  TLegend* leg2 = new TLegend(.125, .6, .425, .85);
637  leg2->AddEntry(hSigNum, "Signal", "l");
638  leg2->AddEntry(hCCNum, "#nu_{#mu} CC", "l");
639  leg2->AddEntry(hNCNum, "NC", "l");
640  leg2->AddEntry(hBeamNum, "Beam #nu_{e}", "l");
641  leg2->SetFillStyle(0);
642  leg2->Draw();
643 
644 }
645 
646 /***************************************************************************/
647 
648 std::string PIDLongName(int pidIdx, bool rhc)
649 {
651  if(pidIdx == -1) ret = "No selection";
652  if(pidIdx == 0) ret = "LEM";
653  if(pidIdx == 1) ret = "LID";
654  if(pidIdx == 2) ret = "CVNe";
655 
656  if(rhc) ret += " (RHC)"; else ret += " (FHC)";
657 
658  return ret;
659 }
660 
661 /***************************************************************************/
662 
663 std::string PIDFileTag(int pidIdx, bool rhc)
664 {
666  if(pidIdx == -1) ret = "nosel";
667  if(pidIdx == 0) ret = "lem";
668  if(pidIdx == 1) ret = "lid";
669  if(pidIdx == 2) ret = "cvn";
670 
671  if(rhc) ret += "_rhc"; else ret += "_fhc";
672 
673  return ret;
674 }
675 
676 /***************************************************************************/
677 
679 {
680 
681  //chose you oscillation parameters:
683  // osc::OscCalcPMNSOpt calc;
684  // calc.SetL(810);
685  // calc.SetRho(2.84);
686  // calc.SetDmsq21(7.53e-5);
687  // calc.SetDmsq32(2.44e-3);
688  // calc.SetTh12(asin(sqrt(0.846))/2);
689  // calc.SetTh13(asin(sqrt(.085))/2);
690  // calc.SetTh23(asin(sqrt(0.999))/2);
691  // calc.SetdCP(0);
692 
693  std::vector<TString> defs{
694  "prod_caf_R17-09-05-prod4recopreview.f_fd_genie_nonswap_fhc_nova_v08_full_v1_addShortSimpleCVN",
695  "prod_caf_R17-09-05-prod4recopreview.f_fd_genie_nonswap_rhc_nova_v08_full_v1_addShortSimpleCVN"};
696  std::vector<TString> swapDefs{
697  "prod_caf_R17-09-05-prod4recopreview.f_fd_genie_fluxswap_fhc_nova_v08_full_v1_addShortSimpleCVN",
698  "prod_caf_R17-09-05-prod4recopreview.f_fd_genie_fluxswap_rhc_nova_v08_full_v1_addShortSimpleCVN"};
699 
700  for(int rhc = false; rhc <= true; ++rhc){
701  //if(rhc) continue; // Don't exist yet for S13-06-18
702 
703  //SpectrumLoader loader("/pnfs/nova/scratch/users/atsaris/cvnFS_eval_fullInt_afterColl_swap/00014*/*/*.root");
704  //SpectrumLoader loaderSwap("/pnfs/nova/scratch/users/atsaris/cvnFS_eval_fullInt_afterColl_nonSwap/00014*/*/*.root");
705 
706  SpectrumLoader loader(defs[rhc].Data());
707  SpectrumLoader loaderSwap(swapDefs[rhc].Data());
708 
709  //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");
710  //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");
711 
712  //how many pids are you optimising cuts on
713  const int kNumPIDs = 3;
714 
715  IPrediction* preds[kNumPIDs];
716 
717  //prediction objects for pid districutions to do the fom optimisation
718  // PIDS:
719  // 1. short Simple
720  // 2. ELU
721  // 3. 2017
722  preds[0] = new PredictionNoExtrap(loader, loaderSwap, "CVN short-simple",
723  Binning::Simple(100, 0.0, 1.),
726 
727  preds[1] = new PredictionNoExtrap(loader, loaderSwap, "CVN ELU",
728  Binning::Simple(110, 0.0, 1.1),
731  kXSecCVWgt2017);
732 
733  preds[2] = new PredictionNoExtrap(loader, loaderSwap, "CVN 2017",
734  Binning::Simple(100, 0.0, 1.),
737  kXSecCVWgt2017);
738 
739  //used as the reference total event count pre-presel, choose wisely
740  PredictionNoExtrap predDenomNoCutRecoE(loader, loaderSwap,
741  "Reconstructed energy (GeV)", kRecoEBinning,
742  kNueEnergy, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment && ISMC, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
743 
744  IPrediction* predDenomPreselRecoE[kNumPIDs];
745  //prediction objects to make efficiency plots for a pre-chosen pid cut
746  predDenomPreselRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
747  "Reconstructed energy (GeV)",
748  kRecoEBinning,
749  kNueEnergy,
750  kNue2017Presel && ISMC, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
751 
752  predDenomPreselRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
753  "Reconstructed energy (GeV)",
754  kRecoEBinning,
755  kNueEnergy,
756  kNue2017Presel && ISMC, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
757 
758  predDenomPreselRecoE[2] = new PredictionNoExtrap(loader, loaderSwap,
759  "Reconstructed energy (GeV)",
760  kRecoEBinning,
761  kNueEnergy,
762  kNue2017Presel && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
763 
764  PredictionNoExtrap predDenomNoCutTrueE(loader, loaderSwap,
765  "True visible energy (GeV)",
766  kRecoEBinning,
767  kTrueEVis, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
768 
769  IPrediction* predDenomPreselTrueE[kNumPIDs];
770  predDenomPreselTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
771  "True visible energy (GeV)",
772  kRecoEBinning,
773  kTrueEVis,
774  kNue2017Presel && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
775  predDenomPreselTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
776  "True visible energy (GeV)",
777  kRecoEBinning,
778  kTrueEVis,
779  kNue2017Presel && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
780 
781  predDenomPreselTrueE[2] = new PredictionNoExtrap(loader, loaderSwap,
782  "True visible energy (GeV)",
783  kRecoEBinning,
784  kTrueEVis,
785  kNue2017Presel && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
786 
787  IPrediction* predsRecoE[kNumPIDs];
788  IPrediction* predsTrueE[kNumPIDs];
789 
790  PredictionNoExtrap predDenomMode(loader, loaderSwap,
791  "", kModeBinning,
792  kMode, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
793 
794  IPrediction* predsMode[kNumPIDs];
795 
796  PredictionNoExtrap predDenomY(loader, loaderSwap,
797  "", kYBinning,
798  kY, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
799 
800  IPrediction* predsY[kNumPIDs];
801 
802  PredictionNoExtrap predDenomXvertex(loader, loaderSwap,
803  "", kXvertexBinning,
804  kXvertex, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
805 
806  IPrediction* predsXvertex[kNumPIDs];
807 
808  PredictionNoExtrap predDenomYvertex(loader, loaderSwap,
809  "", kYvertexBinning,
810  kYvertex, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
811 
812  IPrediction* predsYvertex[kNumPIDs];
813 
814  PredictionNoExtrap predDenomZvertex(loader, loaderSwap,
815  "", kZvertexBinning,
816  kZvertex, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
817 
818  IPrediction* predsZvertex[kNumPIDs];
819 
820 
821  predsRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
822  "Reconstructed energy (GeV)", kRecoEBinning,
823  kNueEnergy, kNue2017Presel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
824 
825  predsRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
826  "Reconstructed energy (GeV)", kRecoEBinning,
827  kNueEnergy, kNue2017Presel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
828 
829  predsRecoE[2] = new PredictionNoExtrap(loader, loaderSwap,
830  "Reconstructed energy (GeV)", kRecoEBinning,
831  kNueEnergy, kNue2017Presel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
832 
833  predsTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
834  "True visible energy (GeV)", kRecoEBinning,
835  kTrueEVis, kNue2017Presel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
836 
837  predsTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
838  "True visible energy (GeV)", kRecoEBinning,
839  kTrueEVis, kNue2017Presel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
840  predsTrueE[2] = new PredictionNoExtrap(loader, loaderSwap,
841  "True visible energy (GeV)", kRecoEBinning,
842  kTrueEVis, kNue2017Presel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
843 
844  predsMode[0] = new PredictionNoExtrap(loader, loaderSwap,
845  "", kModeBinning,
846  kMode, kNue2017Presel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
847 
848  predsMode[1] = new PredictionNoExtrap(loader, loaderSwap,
849  "", kModeBinning,
850  kMode, kNue2017Presel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
851 
852  predsMode[2] = new PredictionNoExtrap(loader, loaderSwap,
853  "", kModeBinning,
854  kMode, kNue2017Presel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
855 
856  predsY[0] = new PredictionNoExtrap(loader, loaderSwap,
857  "True y", kYBinning,
858  kY, kNue2017Presel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
859 
860  predsY[1] = new PredictionNoExtrap(loader, loaderSwap,
861  "True y", kYBinning,
862  kY, kNue2017Presel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
863 
864 
865  predsY[2] = new PredictionNoExtrap(loader, loaderSwap,
866  "True y", kYBinning,
867  kY, kNue2017Presel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
868 
869  predsXvertex[0] = new PredictionNoExtrap(loader, loaderSwap,
870  "X Vertex", kXvertexBinning,
871  kXvertex, kNue2017Presel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
872 
873  predsXvertex[1] = new PredictionNoExtrap(loader, loaderSwap,
874  "X Vertex", kXvertexBinning,
875  kXvertex, kNue2017Presel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
876 
877 
878  predsXvertex[2] = new PredictionNoExtrap(loader, loaderSwap,
879  "X Vertex", kXvertexBinning,
880  kXvertex, kNue2017Presel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
881 
882  predsYvertex[0] = new PredictionNoExtrap(loader, loaderSwap,
883  "Y Vertex", kYvertexBinning,
884  kYvertex, kNue2017Presel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
885 
886  predsYvertex[1] = new PredictionNoExtrap(loader, loaderSwap,
887  "Y Vertex", kYvertexBinning,
888  kYvertex, kNue2017Presel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
889 
890 
891  predsYvertex[2] = new PredictionNoExtrap(loader, loaderSwap,
892  "Y Vertex", kYvertexBinning,
893  kYvertex, kNue2017Presel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
894 
895  predsZvertex[0] = new PredictionNoExtrap(loader, loaderSwap,
896  "Z Vertex", kZvertexBinning,
897  kZvertex, kNue2017Presel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
898 
899  predsZvertex[1] = new PredictionNoExtrap(loader, loaderSwap,
900  "Z Vertex", kZvertexBinning,
901  kZvertex, kNue2017Presel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
902 
903 
904  predsZvertex[2] = new PredictionNoExtrap(loader, loaderSwap,
905  "Z Vertex", kZvertexBinning,
906  kZvertex, kNue2017Presel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
907 
908  loader.Go();
909  loaderSwap.Go();
910 
911  //loop over pids to find FOM cuts and output performance at those points
912  //for(int pidIdx = 0; pidIdx < kNumPIDs; ++pidIdx){
913  for(int pidIdx = 0; pidIdx < kNumPIDs; ++pidIdx){
914  TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
915  GetSpectra(preds[pidIdx], &calc, hSig, hNC, hCC, hBeam, hTotBkg);
916 
917  double cuts[2];
918  const double faCut = 0.75; //<-- set a cut value that the pid will be evaluated at in addition to the FOM optimization points
919  TH1* sigdenomNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
921  Current::kCC,
923 
924  TH1* backnumuNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
926  Current::kCC,
928 
929 
930  TH1* backnueNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
932  Current::kCC,
934 
935  TH1* backncNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
937  Current::kNC,
939 
940 
941  double totalSig=sigdenomNoCut->Integral(0, -1);
942  double totalBckNuCuCC=backnumuNoCut->Integral(0, -1);
943  double totalBckNuECC=backnueNoCut->Integral(0, -1);
944  double totalBckNC=backncNoCut->Integral(0, -1);
945  double totalBackground=totalBckNuCuCC+totalBckNuECC+totalBckNC;
946  double totalsob=totalSig/sqrt(totalBackground);
947  double totalsosb=totalSig/sqrt(totalSig+totalBackground);
948  double totalefficiency=100;
949  double totalpurity=(totalSig*100)/(totalSig+totalBackground);
950 
951  output << std::fixed << std::setprecision(2);
952  output << " &\t " << totalsob << " &\t " << totalsosb
953  << " &\t " << totalSig << " &\t " << totalBackground
954  << " &\t " << totalBckNC
955  << " &\t " << totalBckNuCuCC << " &\t " << totalBckNuECC
956  << " &\t " << totalefficiency << "\\% &\t " << totalpurity << "\\%\\\\"<< std::endl;
957  // & 2.47 & 2.14 & 18.31 & 54.93 & 18.31 & 18.31 18.31 & 100.00\% & 75.00\%\\
958 
959  sigdenomNoCut->SetLineColor(kNueSignalColor);
960 
961  TH1* sigdenomPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
963  Current::kCC,
965  sigdenomPresel->SetLineColor(kNueSignalColor);
966 
967  TH1* backnumuPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
969  Current::kCC,
971 
972  TH1* backnuePresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
974  Current::kCC,
976 
977  TH1* backncPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
979  Current::kNC,
981 
982  double totalPreselSig=sigdenomPresel->Integral(0, -1);
983  double totalPreselBckNuCuCC=backnumuPresel->Integral(0, -1);
984  double totalPreselBckNuECC=backnuePresel->Integral(0, -1);
985  double totalPreselBckNC=backncPresel->Integral(0, -1);
986  double totalPreselBackground=totalPreselBckNuCuCC+totalPreselBckNuECC+totalPreselBckNC;
987  double totalPreselsob=totalPreselSig/sqrt(totalPreselBackground);
988  double totalPreselsosb=totalPreselSig/sqrt(totalPreselSig+totalPreselBackground);
989  double totalPreselefficiency=100;
990  double totalPreselpurity=(totalPreselSig*100)/(totalPreselSig+totalPreselBackground);
991 
992  output << std::fixed << std::setprecision(2);
993  output << " &\t " << totalPreselsob << " &\t " << totalPreselsosb
994  << " &\t " << totalPreselSig << " &\t " << totalPreselBackground
995  << " &\t " << totalPreselBckNC
996  << " &\t " << totalPreselBckNuCuCC << " &\t " << totalPreselBckNuECC
997  << " &\t " << totalPreselefficiency << "\\% &\t " << totalPreselpurity << "\\%\\\\"<< std::endl;
998 
999 
1000  backnumuPresel->Add(backnuePresel);
1001  backnumuPresel->Add(backncPresel);
1002 
1003  backnumuNoCut->Add(backnueNoCut);
1004  backnumuNoCut->Add(backncNoCut);
1005 
1006  FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomNoCut, backnumuNoCut, cuts, faCut, PIDFileTag(pidIdx, rhc), "nocut");
1007  //FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomPresel, backnumuPresel, cuts, faCut, PIDFileTag(pidIdx, rhc), "presel");
1008 
1009  PlotSpectra(preds[pidIdx], &calc, PIDLongName(pidIdx, rhc),9.8);
1010 
1011  Print("pid", PIDFileTag(pidIdx, rhc));
1012  } // end for pidIdx
1013 
1014  PlotSpectra(&predDenomNoCutRecoE, &calc, PIDLongName(-1, rhc), 10);
1015  Print("recoE", PIDFileTag(-1, rhc));
1016  PlotSpectra(&predDenomNoCutTrueE, &calc, PIDLongName(-1, rhc), 10);
1017  Print("trueE", PIDFileTag(-1, rhc));
1018 
1019  } // end for rhc
1020 
1021  output.close();
1022  std::cout << "Wrote to file output.tex\nDone!\n";
1023 }
1024 
1025 #endif
const Binning kModeBinning
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
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 Simulation()
virtual void SetL(double L)=0
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
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
const Var kCVNSS([](const caf::SRProxy *sr){return sr->sel.cvnProd3Train.nueid;})
PID
const Binning kYBinning
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Cut kApplySecondAnalysisMask([](const caf::SRProxy *sr){if(sr->hdr.det!=caf::kFARDET) return true; std::bitset< 14 > binary(sr->hdr.dibmask);std::pair< int, int > planesA=calcFirstLastLivePlane(sr->slc.firstplane, binary);std::pair< int, int > planesB=calcFirstLastLivePlane(sr->slc.lastplane, binary);if((planesA.first!=planesB.first)||(planesA.second!=planesB.second)) return false;return((planesA.second-planesA.first+1)/64 >=4);})
Definition: AnalysisMasks.h:18
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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
const Var kCVNELU([](const caf::SRProxy *sr){return sr->sel.cvn.nueid;})
PID
double maxy
Float_t x1[n_points_granero]
Definition: compare.C:5
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
(&#39;beam &#39;)
Definition: IPrediction.h:15
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
T sqrt(T number)
Definition: d0nt_math.hpp:156
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
void FOMPlot(TH1 *sig, TH1 *bkg, TH1 *bkgCC, TH1 *bkgBeam, TH1 *sig_denom, TH1 *bck_denom, double *bestCuts, double faCut, std::string suffix, std::string eff_suffix)
const Cut kNue2017Presel
Definition: NueCuts2017.h:54
void Print(std::string prefix, std::string name, std::string suffix="")
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 Var kZvertex([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.z;})
Float_t tmp
Definition: plot.C:36
void resetCalc(osc::IOscCalcAdjustable *calc)
const double kPOTMACRO
osc::OscCalcDumb calc
const Color_t kNumuBackgroundColor
Definition: Style.h:30
virtual void SetDmsq32(const T &dmsq32)=0
const Color_t kNueSignalColor
Definition: Style.h:19
void nue_pid_effs_miniprod()
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;})
const Cut ISMC([](const caf::SRProxy *sr){if(sr->hdr.ismc) return true;return false;})
const char * label
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
std::string PIDFileTag(int pidIdx, bool rhc)
Charged-current interactions.
Definition: IPrediction.h:39
const Binning kXvertexBinning
caf::Proxy< float > z
Definition: SRProxy.h:108
const Var kYvertex([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.y;})
TFile * outFile
Definition: PlotXSec.C:135
caf::Proxy< float > x
Definition: SRProxy.h:106
caf::Proxy< caf::SRCVNResult > cvn
Definition: SRProxy.h:1253
caf::Proxy< float > nueid
Definition: SRProxy.h:906
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
void PlotEffs(IPrediction *predNum, IPrediction *predDenom, osc::IOscCalc *calc, std::string title, std::string label, std::string id, double maxy=0)
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
const Binning kYvertexBinning
std::string PIDLongName(int pidIdx, bool rhc)
const Var kXvertex([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.x;})
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const double j
Definition: BetheBloch.cxx:29
const HistDef defs[kNumVars]
Definition: vars.h:15
void GetSpectra(IPrediction *pred, osc::IOscCalc *calc, TH1 *&hSig, TH1 *&hNC, TH1 *&hCC, TH1 *&hBeam, TH1 *&hTotBkg)
loader
Definition: demo0.py:10
TLatex * prelim
Definition: Xsec_final.C:133
const Binning kZvertexBinning
string fFileName
const XML_Char * prefix
Definition: expat.h:380
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
static constexpr Double_t m2
Definition: Munits.h:145
virtual void SetRho(double rho)=0
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
caf::Proxy< float > y
Definition: SRProxy.h:107
virtual void SetTh23(const T &th23)=0
caf::Proxy< bool > ismc
Definition: SRProxy.h:242
std::ofstream output("output.tex")
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
Neutral-current interactions.
Definition: IPrediction.h:40
void PrintEffs(TH1 *h, const char *title)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const Var kCVN2017([](const caf::SRProxy *sr){return sr->sel.cvn2017.nueid;})
PID
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:16
Prediction that just uses FD MC, with no extrapolation.
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
const Binning kRecoEBinning
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
All neutrinos, any flavor.
Definition: IPrediction.h:26
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
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
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void Arrows(TH1 *sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow=false)
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string