plot_nuCrystalBall.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void plot_nuCrystalBall(std::string fname, const TString opt)
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"
18 #include "CAFAna/Core/Utilities.h"
19 #include "CAFAna/Vars/Vars.h"
20 #include "CAFAna/Vars/NueVars.h"
22 #include "CAFAna/Cuts/NueCutsSecondAna.h"
24 
25 #include "TArrow.h"
26 #include "TCanvas.h"
27 #include "TLatex.h"
28 #include "TH1.h"
29 #include "TH2.h"
30 #include "TMarker.h"
31 #include "TGraph.h"
32 
33 #include "TLegend.h"
34 
35 #include <cmath>
36 #include <fstream>
37 #include <iomanip>
38 #include <iostream>
39 
40 #include "Utilities/func/MathUtil.h"
41 
42 using namespace ana;
43 
44 string fFileName="test.txt";
45 double kPOTMACRO = 8.85e20;
46 
47 // 250MeV bins
49 const Binning kModeBinning = Binning::Simple(4, -0.5, 3.5);
50 const Binning kYBinning = Binning::Simple(50, 0, 1);
51 const Binning kYvertexBinning = Binning::Simple(40, -800, 800);
52 const Binning kXvertexBinning = Binning::Simple(40, -800, 800);
54 
55 /***************************************************************************/
56 
57 // Put a "NOvA Simulation" tag in the corner
58 void Simulation()
59 {
60  TLatex* prelim = new TLatex(.9, .95, "NOvA Simulation");
61  prelim->SetTextColor(kGray+1);
62  prelim->SetNDC();
63  prelim->SetTextSize(2/30.);
64  prelim->SetTextAlign(32);
65  prelim->Draw();
66 }
67 
68 // print the current canvas ina variety of formats
70 {
71  name = "plots/"+prefix+"_"+name;
72  if(!suffix.empty()) name += "_"+suffix;
73  Simulation();
74  //gPad->Print((name+".eps").c_str());
75  //gPad->Print((name+".pdf").c_str());
76  //gPad->Print((name+".png").c_str());
77 }
78 
79 /***************************************************************************/
80 // Nominal osc params
82 {
83  calc->SetL(810);
84  calc->SetRho(2.75);
85  calc->SetDmsq21(7.6e-5);
86  calc->SetDmsq32(2.35e-3);
87  calc->SetTh12(asin(sqrt(.87))/2); // PDG
88  calc->SetTh13(asin(sqrt(.095))/2);
89  calc->SetTh23(TMath::Pi()/4);
90  calc->SetdCP(0);
91 }
92 
93 /***************************************************************************/
94 
95 // Draw cuts on plots
96 void Arrows(TH1* sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow = false)
97 {
98  sigId->GetYaxis()->SetRangeUser(0,1.05*sigId->GetMaximum());
99 
100  double maxy = sigId->GetMaximum()/2;
101  if(noarrow)
102  maxy = maxy*2;
103 
104  TLine* lin = new TLine(bestCutSoB, 0, bestCutSoB, maxy);
105  lin->SetLineWidth(3);
106  lin->SetLineColor(kGreen+2);
107  lin->Draw();
108  if(!noarrow){
109  TArrow* arr = new TArrow(bestCutSoB-0.003, maxy, bestCutSoB+.03, maxy, .02, "|>");
110  arr->SetLineWidth(3);
111  arr->SetLineColor(kGreen+2);
112  arr->SetFillColor(kGreen+2);
113  arr->Draw();
114  maxy *= 1.5;
115  }
116 
117  lin = new TLine(bestCutSoSB, 0, bestCutSoSB, maxy);
118  lin->SetLineWidth(3);
119  lin->SetLineStyle(2);
120  lin->SetLineColor(kGreen+3);
121  lin->Draw();
122 
123  if(!noarrow){
124  TArrow* arr = new TArrow(bestCutSoSB-0.003, maxy, bestCutSoSB+.04, maxy, .02, "|>");
125  arr->SetLineWidth(3);
126  arr->SetLineColor(kGreen+3);
127  arr->SetFillColor(kGreen+3);
128  arr->Draw();
129  }
130 
131  lin = new TLine(faCut, 0, faCut, maxy);
132  lin->SetLineWidth(3);
133  lin->SetLineColor(kRed);
134  lin->Draw();
135 
136  if(!noarrow){
137  TArrow* arr = new TArrow(faCut, maxy, faCut+.04, maxy, .02, "|>");
138  arr->SetLineWidth(3);
139  arr->SetLineColor(kRed);
140  arr->SetFillColor(kRed);
141  arr->Draw();
142  }
143 }
144 
145 /***************************************************************************/
146 // Finds best cuts
147 // void FOMPlot(TH1* sig, TH1* bkg, TH1* bkgCC, TH1* bkgBeam,
148 // TH1* sig_denom, TH1* bck_denom,double* bestCuts, double faCut,
149 // std::string suffix, std::string eff_suffix)
150 void FOMPlot(TH1* nue, TH1* nc, TH1* cc, TH1* bkgBeam, TH1* bkg,
151  double* bestCuts, double faCut,
152  std::string suffix, std::string eff_suffix,
153  TString opt)
154 
155 {
156  std::cout << suffix << " " << eff_suffix << ":" << std::endl;
157  if(opt == "nue"){
158  // 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 $\\nu_\\tau$&\t Sig eff &\t Purity\\\\"
159  std::cout << " &\t Cut &\t $s/\\sqrt{b}$ &\t $s/\\sqrt{s+b}$ &\t $\\nu_e$ sig &\t Tot bkg\\\\"
160  << std::endl;
161  }
162  if(opt == "numu"){
163  //std::cout << " &\t Cut &\t $s/\\sqrt{b}$ &\t $s/\\sqrt{s+b}$ &\t $\\nu_\\mu$ sig &\t Tot bkg &\t NC &\t $\\nu_e$ CC &\t Beam $\\nu_e$ &\t $\\nu_\\tau$&\t Sig eff &\t Purity\\\\"
164  std::cout << " &\t Cut &\t $s/\\sqrt{b}$ &\t $s/\\sqrt{s+b}$ &\t $\\nu_\\mu$ sig &\t Tot bkg\\\\"
165  << std::endl;
166  }
167  if(opt == "nus") {
168  //std::cout << " &\t Cut &\t $s/\\sqrt{b}$ &\t $s/\\sqrt{s+b}$ &\t NC &\t Tot bkg &\t $\\nu_e$ &\t $\\nu_\\mu$ CC &\t Beam $\\nu_e$ &\t $\\nu_\\tau$&\t Sig eff &\t Purity\\\\"
169  std::cout << " &\t Cut &\t $s/\\sqrt{b}$ &\t $s/\\sqrt{s+b}$ &\t NC &\t Tot bkg\\\\"
170  << std::endl;
171  }
172 
173  TH1* sig;
174  if(opt == "nue") sig = (TH1*)nue->Clone();
175  if(opt == "numu") sig = (TH1*)cc->Clone();
176  if(opt == "nus") sig = (TH1*)nc->Clone();
177 
178 
179  TH1F* eff = new TH1F(("eff_"+suffix+"_"+eff_suffix).c_str(),
180  "Efficiency vs Cut Value; ;Efficiency (%)",
181  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
182  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
183  TH1F* effxpur = new TH1F(("effxpur_"+suffix+"_"+eff_suffix).c_str(),
184  "Efficiency*Purity vs Cut Value; ;Efficiency*Purity (%)",
185  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
186  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
187  TH1F* pur = new TH1F(("pur_"+suffix+"_"+eff_suffix).c_str(),
188  "Purity vs Cut Value; ;Purity (%)",
189  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
190  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
191  TH1F* fomsob = new TH1F(("sob_"+suffix+"_"+eff_suffix).c_str(),
192  "s/#sqrt{b} vs Cut Value; ;s/#sqrt{b}",
193  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
194  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
195  TH1F* fomsosb = new TH1F(("sosb_"+suffix+"_"+eff_suffix).c_str(),
196  "s/#sqrt{s+b} vs Cut Value; ;s/#sqrt{s+b}",
197  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
198  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
199 
200  eff->GetXaxis()->CenterTitle();
201  eff->SetMaximum(100);
202  eff->GetYaxis()->CenterTitle();
203  pur->GetXaxis()->CenterTitle();
204  pur->GetYaxis()->CenterTitle();
205  fomsob->GetXaxis()->CenterTitle();
206  fomsob->GetYaxis()->CenterTitle();
207  fomsosb->GetXaxis()->CenterTitle();
208  fomsosb->GetYaxis()->CenterTitle();
209 
210  eff->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
211  pur->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
212 
213  fomsob->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
214  fomsosb->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
215  int ngpoints=0;
216  for(int n = 0; n < sig->GetNbinsX()+2; ++n){
217  if(bkg->Integral(n, -1) == 0) break;
218  ngpoints++;
219  }
220  TGraph* effpurGraph=new TGraph (ngpoints);
221  TGraph* rocGraph=new TGraph (ngpoints);
222  double fompur_SoB=0.;
223  double fomeff_SoB=0.;
224  double fompur_SoSB=0.;
225  double fomeff_SoSB=0.;
226  for(int metric = 0; metric < 3; ++metric){
227  double bestSoB = -1;
228  double bestSoSB = -1;
229  double bestSig = 0;
230  double bestBkg = 0;
231  double bestBkgCC = 0;
232  double bestBeam = 0;
233  //double totsig = sig_denom->Integral(0, -1);
234  //double totbck = bck_denom->Integral(0, -1);
235  double fomeff, fompur;
236 
237  for(int n = 0; n < sig->GetNbinsX()+2; ++n){
238  const double nsig = sig->Integral(n, -1);
239  const double nbkg = bkg->Integral(n, -1);
240  // const double nbkgcc = bkgCC->Integral(n, -1);
241  // const double nbkgbeam = bkgBeam->Integral(n, -1);
242 
243  if(nbkg == 0) break;
244 
245  //double neff = nsig*100/(totsig);
246  double npur = nsig*100/(nsig+nbkg);
247 
248  const double fomSoB = nsig/sqrt(nbkg);
249  const double fomSoSB = nsig/sqrt(nsig+nbkg);
250 
251  if(metric == 0){
252  // only fill these hists once.
253  pur->Fill(sig->GetXaxis()->GetBinCenter(n), npur);
254  //eff->Fill(sig->GetXaxis()->GetBinCenter(n), neff);
255  //effxpur->Fill(sig->GetXaxis()->GetBinCenter(n), 100*(neff/100)*(npur/100));
256  //effpurGraph->SetPoint(n,neff,npur);
257  //rocGraph->SetPoint(n,nbkg/totbck,nsig/totsig);
258  fomsob->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoB);
259  fomsosb->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoSB);
260  }
261 
262  if((metric == 0 && fomSoB > bestSoB) ||
263  (metric == 1 && fomSoSB > bestSoSB) ||
264  (metric == 2 && fabs(sig->GetXaxis()->GetBinLowEdge(n) - faCut) < .01)){
265  bestSoB = fomSoB;
266  bestSoSB = fomSoSB;
267  bestSig = nsig;
268  bestBkg = nbkg;
269  // bestBkgCC = nbkgcc;
270  // bestBeam = nbkgbeam;
271  bestCuts[metric] = sig->GetXaxis()->GetBinLowEdge(n);
272  //fomeff = neff;
273  //fompur = npur;
274  }
275  }
276 
277  if(metric == 0) std::cout << "$s/\\sqrt{b}$ opt";
278  if(metric == 1) std::cout << "$s/\\sqrt{s+b}$ opt";
279  if(metric == 2) std::cout << "FA cut";
280 
281  std::cout << std::fixed << std::setprecision(2);
282  std::cout << " &\t " << ((metric == 2) ? faCut : bestCuts[metric])
283  << " &\t " << bestSoB << " &\t " << bestSoSB
284  << " &\t " << bestSig << " &\t " << bestBkg << " &\t "
285  // << " &\t " << bestBkg-bestBkgCC-bestBeam
286  // << " &\t " << bestBkgCC << " &\t " << bestBeam << " &\t "
287  << std::endl;
288  //<< " &\t " << fomeff << "\\% &\t " << fompur << "\\%\\\\"<< std::endl;
289 
290  // if(metric == 0){
291  // fomeff_SoB=fomeff;
292  // fompur_SoB=fompur;
293  // }
294  // if(metric == 1){
295  // fomeff_SoSB=fomeff;
296  // fompur_SoSB=fompur;
297  // }
298 
299  } // end for metric
300 
301  /*
302  new TCanvas;
303  eff->Draw("l");
304  Arrows(eff, bestCuts[0], bestCuts[1], faCut, true);
305  Print("eff", suffix+"_"+eff_suffix);
306 
307  new TCanvas;
308  pur->Draw("l");
309  Arrows(pur, bestCuts[0], bestCuts[1], faCut, true);
310  Print("pur", suffix);
311  */
312 
313 
314  // new TCanvas;
315  // eff->SetTitle("");
316  // eff->GetYaxis()->SetTitle("Percentage");
317  // eff->SetLineColor(kRed);
318  // pur->SetLineColor(kBlue);
319  // effxpur->SetLineColor(kGreen+2);
320  // eff->SetMaximum(110);
321  // eff->Draw("c");
322  // pur->Draw("csame");
323  // effxpur->Draw("csame");
324 
325  /*
326 
327  TLegend* legeffpur = new TLegend(.125, .12, .5, .32);
328  legeffpur->AddEntry(eff, "Cumulative Efficiency", "l");
329  legeffpur->AddEntry(pur, "Cumulative Purity", "l");
330  legeffpur->AddEntry(effxpur, "Cumulative Efficiency #times Purity", "l");
331  legeffpur->SetFillStyle(0);
332  legeffpur->Draw();
333  Print("effpur", suffix+"_"+eff_suffix);
334  */
335 
336  // new TCanvas;
337  // TH2F* hpx = new TH2F("hpx","hpx",100,0,100,100,0,100); // axis range hpx->SetStats(kFALSE); // no statistics
338  // hpx->SetStats(kFALSE);
339  // hpx->SetTitle("");
340  // hpx->GetXaxis()->SetTitle("Efficiency (%)");
341  // hpx->GetYaxis()->SetTitle("Purity (%)");
342  // hpx->GetXaxis()->CenterTitle();
343  // hpx->GetYaxis()->CenterTitle();
344  // hpx->Draw();
345  // for(double i = .5; i < 10; i += .5){
346  // TGraph* c = new TGraph;
347  // for(double j = .01; j < 101; j += .1){
348  // c->SetPoint(c->GetN(), j, (100*i*i)/j);
349  // }
350  // c->SetLineWidth(2);
351  // c->SetLineColor(kGray+1);
352  // c->SetLineStyle(7);
353  // c->Draw("l same");
354  // }
355 
356 
357  // TMarker* m1 = new TMarker(fomeff_SoB, fompur_SoB, kFullCircle);
358  // m1->SetMarkerSize(1.5*m1->GetMarkerSize());
359  // m1->Draw();
360 
361  // TMarker* m2 = new TMarker(fomeff_SoSB, fompur_SoSB, kFullTriangleUp);
362  // m2->SetMarkerSize(1.5*m2->GetMarkerSize());
363  // m2->Draw();
364 
365  // effpurGraph->SetLineWidth(2);
366  // effpurGraph->Draw("C");
367  // Print("effpurgraph", suffix+"_"+eff_suffix);
368 
369  // new TCanvas;
370  // TH2F* roch = new TH2F("roch","roch",100,0,1,100,0,1); // axis range roch->SetStats(kFALSE); // no statistics
371  // roch->SetStats(kFALSE);
372  // roch->SetTitle("");
373  // roch->GetXaxis()->SetTitle("False Positive Rate");
374  // roch->GetYaxis()->SetTitle("True Positive Rate");
375  // roch->GetXaxis()->CenterTitle();
376  // roch->GetYaxis()->CenterTitle();
377  // roch->Draw();
378  // TLine* coin = new TLine(0, 0, 1, 1);
379  // coin->SetLineStyle(2);
380  // coin->Draw();
381  // rocGraph->SetLineWidth(2);
382  // rocGraph->Draw("C");
383  // Print("rocGraph", suffix+"_"+eff_suffix);
384 
385  // TFile* outFile=new TFile((suffix+"_"+eff_suffix+".root").c_str(),"RECREATE");
386  // rocGraph->Write("rocGraph");
387  // effpurGraph->Write("effpurGraph");
388  // outFile->Close();
389 
390  new TCanvas;
391  fomsob->Draw("l");
392  Arrows(fomsob, bestCuts[0], bestCuts[1], faCut, true);
393  Print("fomsob", suffix);
394 
395  new TCanvas;
396  fomsosb->Draw("l");
397  Arrows(fomsosb, bestCuts[0], bestCuts[1], faCut, true);
398  Print("fomsosb", suffix);
399 
400 }
401 
402 /***************************************************************************/
403 
405  TH1*& hNue, TH1*& hNC, TH1*& hCC, TH1*& hBeam, TH1*&hTotBkg,
406  TString opt, bool rhc=false)
407 {
408 
409  std::cout << "RHC 2: " << rhc << std::endl;
410  if(rhc) kPOTMACRO = 8.1e20;
411  std::cout << "POT: " << kPOTMACRO << std::endl;
412  hNue = pred->PredictComponent(calc,
414  Current::kCC,
416  hNue->SetLineColor(kNueSignalColor);
417 
418  hNC = pred->PredictComponent(calc,
420  Current::kNC,
422  hNC->SetLineColor(kNCBackgroundColor);
423 
424  hCC = pred->PredictComponent(calc,
426  Current::kCC,
428  hCC->SetLineColor(kNumuBackgroundColor);
429 
430  hBeam = pred->PredictComponent(calc,
432  Current::kCC,
434  hBeam->SetLineColor(kBeamNueBackgroundColor);
435 
436  if(opt == "nue"){
437  hTotBkg = (TH1*)hNC->Clone(UniqueName().c_str());
438  hTotBkg->Add(hCC);
439  hTotBkg->Add(hBeam);
440  }
441  if(opt == "numu"){
442  hTotBkg = (TH1*)hNC->Clone(UniqueName().c_str());
443  hTotBkg->Add(hNue);
444  hTotBkg->Add(hBeam);
445  }
446  if(opt == "nus"){
447  hTotBkg = (TH1*)hNue->Clone(UniqueName().c_str());
448  hTotBkg->Add(hCC);
449  hTotBkg->Add(hBeam);
450  }
451 
452  hTotBkg->SetLineColor(kTotalMCColor);
453 
454 }
455 
456 /***************************************************************************/
457 
458 void PlotSpectra(IPrediction* pred, osc::IOscCalc* calc, std::string title, TString opt, double maxy = 0, bool isrhc=false)
459 {
460  TH1 *hNue, *hNC, *hCC, *hBeam, *hTotBkg, *hSig;
461  std::cout << "RHC 3 : " << isrhc << std::endl;
462  GetSpectra(pred, calc, hNue, hNC, hCC, hBeam, hTotBkg, opt, isrhc);
463 
464  if(opt == "nue") hSig = (TH1*)hNue->Clone();
465  if(opt == "numu") hSig = (TH1*)hCC->Clone();
466  if(opt == "nus") hSig = (TH1*)hNC->Clone();
467  new TCanvas;
468  hSig->SetTitle("");//title.c_str());
469  double kRealPot=8.85e20;
470  hSig->GetYaxis()->SetTitle(TString::Format("Events / %.02lf#times10^{20} POT", kRealPot/1e20));
471  //hSig->GetYaxis()->SetTitle("Arbitary Units");
472  hSig->GetXaxis()->CenterTitle();
473  hSig->GetYaxis()->CenterTitle();
474  hSig->Draw("hist");
475  if(maxy) hSig->GetYaxis()->SetRangeUser(0, maxy);
476  hNC->Draw("hist same");
477  hCC->Draw("hist same");
478  hBeam->Draw("hist same");
479  TLegend* leg2 = new TLegend(.125, .55, .5, .85);
480  if(opt == "nue"){
481  leg2->AddEntry(hSig,"#nu_{e} signal", "l");
482  leg2->AddEntry(hCC, "#nu_{#mu} bkgd.", "l");
483  leg2->AddEntry(hNC, "NC background", "l");
484  }
485  if(opt == "numu"){
486  leg2->AddEntry(hSig, "#nu_{#mu} signal", "l");
487  leg2->AddEntry(hNue, "#nu_{e} bkgd.", "l");
488  leg2->AddEntry(hNC, "NC background", "l");
489  }
490  if(opt == "nus"){
491  leg2->AddEntry(hSig, "NC signal", "l");
492  leg2->AddEntry(hCC, "#nu_{#mu} bkgd.", "l");
493  leg2->AddEntry(hNue, "#nu_{e} bkgd.", "l");
494  }
495  leg2->AddEntry(hBeam, "Beam #nu_{e} background", "l");
496  leg2->SetFillStyle(0);
497  leg2->Draw();
498 
499  TFile *f = new TFile("plots/spectra_"+opt+title+".root", "update");
500  hSig->Write("signal");
501  hTotBkg->Write("background");
502  f->Close();
503 
504  //gPad->Print("plots/spectra"+opt+title+".root");
505 
506 
507  // static bool once = true;
508  // if(once){
509  // once = false;
510  // TVirtualPad* bak = gPad;
511  // new TCanvas;
512  // TLegend* leg = new TLegend(.1, .1, .9, .9);
513  // leg->AddEntry(hSig, "#nu_{e} signal", "l");
514  // leg->AddEntry(hNC, "NC background", "l");
515  // leg->AddEntry(hCC, "#nu_{#mu} background", "l");
516  // leg->AddEntry(hBeam, "Beam #nu_{e} background", "l");
517  // leg->Draw();
518  // gPad->Print("plots/eff_leg.pdf");
519  // gPad->Print("plots/eff_leg.eps");
520  // gPad->Print("plots/eff_leg.png");
521  // bak->cd();
522  // }
523 }
524 
525 void PrintEffs(TH1* h, const char* title)
526 {
527  std::cout << title;
528  for(int i = 0; i < h->GetNbinsX()+2; ++i){
529  if(h->GetBinCenter(i) > .75 &&
530  h->GetBinCenter(i) < 3.5){
531  std::cout << " &\t " << h->GetBinContent(i);
532  }
533  }
534  std::cout << "\\\\" << std::endl;
535 }
536 
537 /***************************************************************************/
538 
539 void PlotEffs(IPrediction* predNum, IPrediction* predDenom,
541  std::string label, std::string id, TString opt, double maxy = 0)
542 {
543  TH1 *hSigNum, *hNCNum, *hCCNum, *hBeamNum, *hTotBkgNum;
544  GetSpectra(predNum, calc, hSigNum, hNCNum, hCCNum, hBeamNum, hTotBkgNum, opt);
545 
546  TH1 *hSigDenom, *hNCDenom, *hCCDenom, *hBeamDenom, *hTotBkgDenom;
547  GetSpectra(predDenom, calc, hSigDenom, hNCDenom, hCCDenom, hBeamDenom, hTotBkgDenom, opt);
548 
549  hSigDenom->SetLineStyle(2);
550  hNCDenom->SetLineStyle(2);
551  hBeamDenom->SetLineStyle(2);
552  hTotBkgDenom->SetLineStyle(2);
553  hCCDenom->SetLineStyle(2);
554 
555  new TCanvas();
556 
557  if( label == "mode" ){
558  hNCDenom->GetXaxis()->SetLabelSize(0.06);
559  hNCDenom->GetXaxis()->SetBinLabel(1, "QE");
560  hNCDenom->GetXaxis()->SetBinLabel(2, "Res");
561  hNCDenom->GetXaxis()->SetBinLabel(3, "DIS");
562  hNCDenom->GetXaxis()->SetBinLabel(4, "Coh");
563 
564  hSigNum->GetXaxis()->SetLabelSize(0.06);
565  hSigNum->GetXaxis()->SetBinLabel(1, "QE");
566  hSigNum->GetXaxis()->SetBinLabel(2, "Res");
567  hSigNum->GetXaxis()->SetBinLabel(3, "DIS");
568  hSigNum->GetXaxis()->SetBinLabel(4, "Coh");
569 
570  }
571 
572  hNCDenom->Draw("hist");
573  hCCDenom->Draw("hist same");
574  hBeamDenom->Draw("hist same");
575  hSigDenom->Draw("hist same");
576  Print(label, id, "nocut");
577 
578  new TCanvas();
579  if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
580  //hSigNum->GetYaxis()->SetTitle("Arbitary Units");
581  hSigNum->GetYaxis()->CenterTitle();
582  hSigNum->Draw("hist");
583  hNCNum->Draw("hist same");
584  hBeamNum->Draw("hist same");
585  hCCNum->Draw("hist same");
586  Print(label, id, "cut");
587 
588  new TCanvas;
589 
590  std::cout << "Total efficiency"
591  << " sig: " << 100*hSigNum->Integral(0, -1)/hSigDenom->Integral(0, -1)
592  << " NC: " << 100*hNCNum->Integral(0, -1)/hNCDenom->Integral(0, -1)
593  << " CC: " << 100*hCCNum->Integral(0, -1)/hCCDenom->Integral(0, -1)
594  << " Beam: " << 100*hBeamNum->Integral(0, -1)/hBeamDenom->Integral(0, -1) << std::endl;
595 
596  hSigNum->Divide(hSigDenom);
597  hSigNum->Scale(100);
598  hSigNum->SetTitle("");
599  hSigNum->GetYaxis()->SetTitle("Efficiency (%)");
600  hSigNum->GetYaxis()->SetRangeUser(0, 100);
601  //if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
602  hSigNum->GetXaxis()->CenterTitle();
603  hSigNum->GetYaxis()->CenterTitle();
604  hSigNum->Draw("hist");
605 
606  PrintEffs(hSigNum, "$\\nu_e$ sig");
607 
608  hNCNum->Divide(hNCDenom);
609  hNCNum->Scale(100);
610  hNCNum->Draw("hist same");
611 
612  PrintEffs(hNCNum, "NC");
613 
614  hCCNum->Divide(hCCDenom);
615  hCCNum->Scale(100);
616  hCCNum->Draw("hist same");
617 
618  PrintEffs(hCCNum, "$\\nu_\\mu$ CC");
619 
620  hBeamNum->Divide(hBeamDenom);
621  hBeamNum->Scale(100);
622  hBeamNum->Draw("hist same");
623 
624  PrintEffs(hBeamNum, "beam $\\nu_e$");
625 
626  double x1=125;
627  double y1=.6;
628  double x2=.425;
629  double y2=.85;
630 
631  if(title=="y"){
632  x1=0.6;
633  y1=.6;
634  x2=.9;
635  y2=.85;
636  }
637 
638  TLegend* leg2 = new TLegend(.125, .6, .425, .85);
639  leg2->AddEntry(hSigNum, "Signal", "l");
640  leg2->AddEntry(hCCNum, "#nu_{#mu} CC", "l");
641  leg2->AddEntry(hNCNum, "NC", "l");
642  leg2->AddEntry(hBeamNum, "Beam #nu_{e}", "l");
643  leg2->SetFillStyle(0);
644  leg2->Draw();
645 
646 }
647 
648 /***************************************************************************/
649 
650 std::string PIDLongName(int pidIdx, bool rhc)
651 {
653  if(pidIdx == -1) ret = "No selection";
654 
655  if(pidIdx == 0) ret = "CVNeShortSimple";
656  if(pidIdx == 1) ret = "CVNeELU";
657  if(pidIdx == 2) ret = "CVNe2017";
658  if(pidIdx == 3) ret = "ReMID";
659 
660  if(rhc) ret += "_RHC"; else ret += "_FHC";
661 
662  return ret;
663 }
664 
665 /***************************************************************************/
666 
667 std::string PIDFileTag(int pidIdx, bool rhc)
668 {
670  if(pidIdx == -1) ret = "nosel";
671  if(pidIdx == 0) ret = "shortsimple";
672  if(pidIdx == 1) ret = "elu";
673  if(pidIdx == 2) ret = "classic";
674  if(pidIdx == 3) ret = "remid";
675 
676  if(rhc) ret += "_rhc"; else ret += "_fhc";
677 
678  return ret;
679 }
680 
681 /***************************************************************************/
682 
684 {
685 
686  TFile fin(fname.c_str());
687  if(fin.IsZombie()){
688  std::cerr << "File: " << fname << " does not exist (or is Zombie). Double-check input filename." << std::endl;
689  exit(0);
690  }
691 
692  //chose you oscillation parameters:
694 
695  std::vector<TString> pidNames;
696  std::vector<TString> cutNames;
697 
698  if(opt == "nue"){
699  pidNames.push_back("CVNeShortSimple");
700  pidNames.push_back("CVNeELU");
701  pidNames.push_back("CVNe2017");
702 
703  //cutNames.push_back("Nue2017FDMinusCVN");
704  cutNames.push_back("Nue2017CorePresel");
705  }
706 
707  if(opt == "numu"){
708  pidNames.push_back("CVNmuShortSimple");
709  pidNames.push_back("CVNmuELU");
710  pidNames.push_back("CVNmu2017");
711  pidNames.push_back("RemID");
712 
713  cutNames.push_back("Numu2017MinusCVN");
714  }
715 
716  if(opt == "nus"){
717  pidNames.push_back("CVNsShortSimple");
718  pidNames.push_back("CVNsELU");
719  pidNames.push_back("CVNs2017");
720 
721  cutNames.push_back("Nus2017MinusCVN");
722  }
723 
724  struct PredDef
725  {
726  const IPrediction *pred;
727  const TString cutName;
728  const TString varName;
729  const TString fhcName;
730  };
731 
732  std::vector<PredDef> predDefs;
733 
734  for(int rhc = false; rhc <= true; ++rhc){
735 
736  TString polarity = "fhc";
737  if(rhc) // Just test fhc first
738  polarity = "rhc";
739 
740 
741  TString xp = "pred_nxp";
742 
743  //how many pids are you optimising cuts on
744  const int kNumPIDs = 3;
745 
746  for(int cutIdx = 0; cutIdx < 1; ++cutIdx) {
747 
748  for(int pidIdx = 0; pidIdx < kNumPIDs; ++pidIdx) {
749 
750  IPrediction* pred =
751  LoadFromFile<IPrediction>(fname,
752  (xp + "_" +
753  cutNames[cutIdx] + "_" +
754  pidNames[pidIdx] + "_" +
755  polarity).Data()).release();
756 
757  //loop over pids to find FOM cuts and output performance at those points
758  //for(int pidIdx = 0; pidIdx < kNumPIDs; ++pidIdx){
759  TH1 *hNue, *hNC, *hCC, *hBeam, *hTotBkg;
760  GetSpectra(pred, &calc, hNue, hNC, hCC, hBeam, hTotBkg, opt);
761 
762  double cuts[2];
763  double faCut = 0.;
764  if(opt == "nue") faCut = 0.75;
765  if(opt == "numu") faCut = 0.5;
766  if(opt == "nus") faCut = 0.2;
767  //const double faCut = 0.; //<- set a cut value that the pid will be evaluated at in addition to the FOM optimization points
768  // TH1* sigdenomNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
769  // Flavors::kNuMuToNuE,
770  // Current::kCC,
771  // Sign::kBoth).ToTH1(kPOTMACRO);
772 
773  // TH1* backnumuNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
774  // Flavors::kNuMuToNuMu,
775  // Current::kCC,
776  // Sign::kBoth).ToTH1(kPOTMACRO);
777 
778 
779  // TH1* backnueNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
780  // Flavors::kNuEToNuE,
781  // Current::kCC,
782  // Sign::kBoth).ToTH1(kPOTMACRO);
783 
784  // TH1* backncNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
785  // Flavors::kAll,
786  // Current::kNC,
787  // Sign::kBoth).ToTH1(kPOTMACRO);
788 
789 
790  // double totalSig=sigdenomNoCut->Integral(0, -1);
791  // double totalBckNuCuCC=backnumuNoCut->Integral(0, -1);
792  // double totalBckNuECC=backnueNoCut->Integral(0, -1);
793  // double totalBckNC=backncNoCut->Integral(0, -1);
794  // double totalBackground=totalBckNuCuCC+totalBckNuECC+totalBckNC;
795  // double totalsob=totalSig/sqrt(totalBackground);
796  // double totalsosb=totalSig/sqrt(totalSig+totalBackground);
797  // double totalefficiency=100;
798  // double totalpurity=(totalSig*100)/(totalSig+totalBackground);
799 
800  // std::cout << std::fixed << std::setprecision(2);
801  // std::cout << " &\t " << totalsob << " &\t " << totalsosb
802  // << " &\t " << totalSig << " &\t " << totalBackground
803  // << " &\t " << totalBckNC
804  // << " &\t " << totalBckNuCuCC << " &\t " << totalBckNuECC
805  // << " &\t " << totalefficiency << "\\% &\t " << totalpurity << "\\%\\\\"<< std::endl;
806  // // & 2.47 & 2.14 & 18.31 & 54.93 & 18.31 & 18.31 18.31 & 100.00\% & 75.00\%\ \
807 
808  // sigdenomNoCut->SetLineColor(kNueSignalColor);
809 
810  // TH1* sigdenomPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
811  // Flavors::kNuMuToNuE,
812  // Current::kCC,
813  // Sign::kBoth).ToTH1(kPOTMACRO);
814  // sigdenomPresel->SetLineColor(kNueSignalColor);
815 
816  // TH1* backnumuPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
817  // Flavors::kNuMuToNuMu,
818  // Current::kCC,
819  // Sign::kBoth).ToTH1(kPOTMACRO);
820 
821  // TH1* backnuePresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
822  // Flavors::kNuEToNuE,
823  // Current::kCC,
824  // Sign::kBoth).ToTH1(kPOTMACRO);
825 
826  // TH1* backncPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
827  // Flavors::kAll,
828  // Current::kNC,
829  // Sign::kBoth).ToTH1(kPOTMACRO);
830 
831  // double totalPreselSig=sigdenomPresel->Integral(0, -1);
832  // double totalPreselBckNuCuCC=backnumuPresel->Integral(0, -1);
833  // double totalPreselBckNuECC=backnuePresel->Integral(0, -1);
834  // double totalPreselBckNC=backncPresel->Integral(0, -1);
835  // double totalPreselBackground=totalPreselBckNuCuCC+totalPreselBckNuECC+totalPreselBckNC;
836  // double totalPreselsob=totalPreselSig/sqrt(totalPreselBackground);
837  // double totalPreselsosb=totalPreselSig/sqrt(totalPreselSig+totalPreselBackground);
838  // double totalPreselefficiency=100;
839  // double totalPreselpurity=(totalPreselSig*100)/(totalPreselSig+totalPreselBackground);
840 
841  // std::cout << std::fixed << std::setprecision(2);
842  // std::cout << " &\t " << totalPreselsob << " &\t " << totalPreselsosb
843  // << " &\t " << totalPreselSig << " &\t " << totalPreselBackground
844  // << " &\t " << totalPreselBckNC
845  // << " &\t " << totalPreselBckNuCuCC << " &\t " << totalPreselBckNuECC
846  // << " &\t " << totalPreselefficiency << "\\% &\t " << totalPreselpurity << "\\%\\\\"<< std::endl;
847 
848 
849  // backnumuPresel->Add(backnuePresel);
850  // backnumuPresel->Add(backncPresel);
851 
852  // backnumuNoCut->Add(backnueNoCut);
853  // backnumuNoCut->Add(backncNoCut);
854 
855  FOMPlot(hNue, hNC, hCC, hBeam, hTotBkg, cuts, faCut, PIDFileTag(pidIdx, rhc), "nocut", opt);
856  //FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomPresel, backnumuPresel, cuts, faCut, PIDFileTag(pidIdx, rhc), "presel");
857  std::cout << "RHC: " << rhc << std::endl;
858  PlotSpectra(pred, &calc, PIDLongName(pidIdx, rhc), opt, 0, rhc);
859 
860  //Print("pid", PIDFileTag(pidIdx, rhc));
861  } // end for pidIdx
862  }
863 
864  // PlotSpectra(&predDenomNoCutRecoE, &calc, PIDLongName(-1, rhc), 10);
865  // Print("recoE", PIDFileTag(-1, rhc));
866  // PlotSpectra(&predDenomNoCutTrueE, &calc, PIDLongName(-1, rhc), 10);
867  // Print("trueE", PIDFileTag(-1, rhc));
868 
869  } // end for rhc
870 }
871 
872 #endif
const Binning kModeBinning
TString fin
Definition: Style.C:24
const XML_Char * name
Definition: expat.h:151
const TString cutName
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
string fFileName
const Binning kYBinning
void PlotEffs(IPrediction *predNum, IPrediction *predDenom, osc::IOscCalc *calc, std::string title, std::string label, std::string id, TString opt, double maxy=0)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const TString polarity
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#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
Float_t x1[n_points_granero]
Definition: compare.C:5
Struct to hold prediction information.
(&#39;beam &#39;)
Definition: IPrediction.h:15
void Print(std::string prefix, std::string name, std::string suffix="")
T sqrt(T number)
Definition: d0nt_math.hpp:156
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
double kPOTMACRO
virtual void SetTh13(const T &th13)=0
OStream cerr
Definition: OStream.cxx:7
void resetCalc(osc::IOscCalcAdjustable *calc)
osc::OscCalcDumb calc
const Color_t kNumuBackgroundColor
Definition: Style.h:30
virtual void SetDmsq32(const T &dmsq32)=0
std::string PIDFileTag(int pidIdx, bool rhc)
const Color_t kNueSignalColor
Definition: Style.h:19
const char * label
Charged-current interactions.
Definition: IPrediction.h:39
const Binning kXvertexBinning
void GetSpectra(IPrediction *pred, osc::IOscCalc *calc, TH1 *&hNue, TH1 *&hNC, TH1 *&hCC, TH1 *&hBeam, TH1 *&hTotBkg, TString opt, bool rhc=false)
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
void Arrows(TH1 *sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow=false)
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
TLatex * prelim
Definition: Xsec_final.C:133
const Binning kZvertexBinning
const XML_Char * prefix
Definition: expat.h:380
void Simulation()
OStream cout
Definition: OStream.cxx:6
void FOMPlot(TH1 *nue, TH1 *nc, TH1 *cc, TH1 *bkgBeam, TH1 *bkg, double *bestCuts, double faCut, std::string suffix, std::string eff_suffix, TString opt)
virtual void SetRho(double rho)=0
void plot_nuCrystalBall(std::string fname, const TString opt)
void cc()
Definition: test_ana.C:28
virtual void SetTh23(const T &th23)=0
void PrintEffs(TH1 *h, const char *title)
const TString fhcName
exit(0)
Neutral-current interactions.
Definition: IPrediction.h:40
enum BeamMode nc
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
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
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
enum BeamMode kGreen
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:107
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
const TString varName
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
std::string PIDLongName(int pidIdx, bool rhc)
enum BeamMode string