Functions | Variables
nue_pid_effs.C File Reference
#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Cuts/NueCutsFirstAna.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/NueVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "StandardRecord/Proxy/SRProxy.h"
#include "TArrow.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TH1.h"
#include "TH2.h"
#include "TMarker.h"
#include "TGraph.h"
#include "TLegend.h"
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include "Utilities/func/MathUtil.h"

Go to the source code of this file.

Functions

void Simulation ()
 
void Print (std::string prefix, std::string name, std::string suffix="")
 
void resetCalc (osc::IOscCalcAdjustable *calc)
 
void Arrows (TH1 *sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow=false)
 
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)
 
void GetSpectra (IPrediction *pred, osc::IOscCalc *calc, TH1 *&hSig, TH1 *&hNC, TH1 *&hCC, TH1 *&hBeam, TH1 *&hTotBkg)
 
void PlotSpectra (IPrediction *pred, osc::IOscCalc *calc, std::string title, double maxy=0)
 
void PrintEffs (TH1 *h, const char *title)
 
void PlotEffs (IPrediction *predNum, IPrediction *predDenom, osc::IOscCalc *calc, std::string title, std::string label, std::string id, double maxy=0)
 
std::string PIDLongName (int pidIdx, bool rhc)
 
std::string PIDFileTag (int pidIdx, bool rhc)
 
void nue_pid_effs ()
 

Variables

string fFileName ="test.txt"
 
const double kPOTMACRO = 6.754e20
 
const Binning kRecoEBinning = Binning::Simple(16, 0, 4)
 
const Binning kModeBinning = Binning::Simple(4, -0.5, 3.5)
 
const Binning kYBinning = Binning::Simple(50, 0, 1)
 
const Binning kYvertexBinning = Binning::Simple(40, -800, 800)
 
const Binning kXvertexBinning = Binning::Simple(40, -800, 800)
 
const Binning kZvertexBinning = Binning::Simple(40, 0, 6000)
 
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;})
 
const Cut ISMC ([](const caf::SRProxy *sr){if(sr->hdr.ismc) return true;return false;})
 
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 Var kXvertex ([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.x;})
 
const Var kYvertex ([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.y;})
 
const Var kZvertex ([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.z;})
 

Function Documentation

void Arrows ( TH1 *  sigId,
double  bestCutSoB,
double  bestCutSoSB,
double  faCut,
bool  noarrow = false 
)

Definition at line 141 of file nue_pid_effs.C.

References kGreen, kRed, and maxy.

Referenced by FOMPlot().

142 {
143  sigId->GetYaxis()->SetRangeUser(0,1.05*sigId->GetMaximum());
144 
145  double maxy = sigId->GetMaximum()/2;
146  if(noarrow)
147  maxy = maxy*2;
148 
149  TLine* lin = new TLine(bestCutSoB, 0, bestCutSoB, maxy);
150  lin->SetLineWidth(3);
151  lin->SetLineColor(kGreen+2);
152  lin->Draw();
153  if(!noarrow){
154  TArrow* arr = new TArrow(bestCutSoB-0.003, maxy, bestCutSoB+.03, maxy, .02, "|>");
155  arr->SetLineWidth(3);
156  arr->SetLineColor(kGreen+2);
157  arr->SetFillColor(kGreen+2);
158  arr->Draw();
159  maxy *= 1.5;
160  }
161 
162 
163 
164  lin = new TLine(bestCutSoSB, 0, bestCutSoSB, maxy);
165  lin->SetLineWidth(3);
166  lin->SetLineStyle(2);
167  lin->SetLineColor(kGreen+3);
168  lin->Draw();
169 
170  if(!noarrow){
171  TArrow* arr = new TArrow(bestCutSoSB-0.003, maxy, bestCutSoSB+.04, maxy, .02, "|>");
172  arr->SetLineWidth(3);
173  arr->SetLineColor(kGreen+3);
174  arr->SetFillColor(kGreen+3);
175  arr->Draw();
176  }
177 
178 
179  lin = new TLine(faCut, 0, faCut, maxy);
180  lin->SetLineWidth(3);
181  lin->SetLineColor(kRed);
182  //lin->Draw();
183 
184  if(!noarrow){
185  TArrow* arr = new TArrow(faCut, maxy, faCut+.04, maxy, .02, "|>");
186  arr->SetLineWidth(3);
187  arr->SetLineColor(kRed);
188  arr->SetFillColor(kRed);
189  //arr->Draw();
190  }
191 }
kRed
enum BeamMode kRed
maxy
double maxy
Definition: plot_datapredictions.C:15
kGreen
enum BeamMode kGreen
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 
)

Definition at line 195 of file nue_pid_effs.C.

References Arrows(), plot_validation_datamc::c, om::cout, nd_projection_maker::eff, allTimeWatchdog::endl, stan::math::fabs(), MECModelEnuComparisons::i, calib::j, kBlue, kGreen, kRed, Munits::m2, productionTest::metric, getGoodRuns4SAM::n, outFile, Print(), PandAna.Demos.pi0_spectra::pur, and std::sqrt().

Referenced by nue_pid_effs().

198 {
199  std::cout << suffix << " " << eff_suffix << ":" << std::endl;
200  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 Beam $\\nu_\\tau$&\t Sig eff &\t Purity\\\\"
201  << std::endl;
202 
203  TH1F* eff = new TH1F(("eff_"+suffix+"_"+eff_suffix).c_str(),
204  "Efficiency vs Cut Value; ;Efficiency (%)",
205  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
206  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
207  TH1F* effxpur = new TH1F(("effxpur_"+suffix+"_"+eff_suffix).c_str(),
208  "Efficiency*Purity vs Cut Value; ;Efficiency*Purity (%)",
209  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
210  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
211  TH1F* pur = new TH1F(("pur_"+suffix+"_"+eff_suffix).c_str(),
212  "Purity vs Cut Value; ;Purity (%)",
213  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
214  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
215  TH1F* fomsob = new TH1F(("sob_"+suffix+"_"+eff_suffix).c_str(),
216  "s/#sqrt{b} vs Cut Value; ;s/#sqrt{b}",
217  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
218  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
219  TH1F* fomsosb = new TH1F(("sosb_"+suffix+"_"+eff_suffix).c_str(),
220  "s/#sqrt{s+b} vs Cut Value; ;s/#sqrt{s+b}",
221  sig->GetNbinsX(), sig->GetXaxis()->GetBinLowEdge(1),
222  sig->GetXaxis()->GetBinUpEdge(sig->GetNbinsX()));
223 
224  eff->GetXaxis()->CenterTitle();
225  eff->SetMaximum(100);
226  eff->GetYaxis()->CenterTitle();
227  pur->GetXaxis()->CenterTitle();
228  pur->GetYaxis()->CenterTitle();
229  fomsob->GetXaxis()->CenterTitle();
230  fomsob->GetYaxis()->CenterTitle();
231  fomsosb->GetXaxis()->CenterTitle();
232  fomsosb->GetYaxis()->CenterTitle();
233 
234  eff->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
235  pur->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
236 
237  fomsob->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
238  fomsosb->GetXaxis()->SetTitle(sig->GetXaxis()->GetTitle());
239  int ngpoints=0;
240  for(int n = 0; n < sig->GetNbinsX()+2; ++n){
241  if(bkg->Integral(n, -1) == 0) break;
242  ngpoints++;
243  }
244  TGraph* effpurGraph=new TGraph (ngpoints);
245  TGraph* rocGraph=new TGraph (ngpoints);
246  double fompur_SoB=0.;
247  double fomeff_SoB=0.;
248  double fompur_SoSB=0.;
249  double fomeff_SoSB=0.;
250  for(int metric = 0; metric < 3; ++metric){
251  double bestSoB = -1;
252  double bestSoSB = -1;
253  double bestSig = 0;
254  double bestBkg = 0;
255  double bestBkgCC = 0;
256  double bestBeam = 0;
257  double totsig = sig_denom->Integral(0, -1);
258  double totbck = bck_denom->Integral(0, -1);
259  double fomeff, fompur;
260 
261  for(int n = 0; n < sig->GetNbinsX()+2; ++n){
262  const double nsig = sig->Integral(n, -1);
263  const double nbkg = bkg->Integral(n, -1);
264  const double nbkgcc = bkgCC->Integral(n, -1);
265  const double nbkgbeam = bkgBeam->Integral(n, -1);
266 
267  if(nbkg == 0) break;
268 
269  double neff = nsig*100/(totsig);
270  double npur = nsig*100/(nsig+nbkg);
271 
272  const double fomSoB = nsig/sqrt(nbkg);
273  const double fomSoSB = nsig/sqrt(nsig+nbkg);
274 
275  if(metric == 0){
276  // only fill these hists once.
277  pur->Fill(sig->GetXaxis()->GetBinCenter(n), npur);
278  eff->Fill(sig->GetXaxis()->GetBinCenter(n), neff);
279  effxpur->Fill(sig->GetXaxis()->GetBinCenter(n), 100*(neff/100)*(npur/100));
280  effpurGraph->SetPoint(n,neff,npur);
281  rocGraph->SetPoint(n,nbkg/totbck,nsig/totsig);
282  fomsob->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoB);
283  fomsosb->Fill(sig->GetXaxis()->GetBinCenter(n), fomSoSB);
284  }
285 
286  if((metric == 0 && fomSoB > bestSoB) ||
287  (metric == 1 && fomSoSB > bestSoSB) ||
288  (metric == 2 && fabs(sig->GetXaxis()->GetBinLowEdge(n) - faCut) < .01)){
289  bestSoB = fomSoB;
290  bestSoSB = fomSoSB;
291  bestSig = nsig;
292  bestBkg = nbkg;
293  bestBkgCC = nbkgcc;
294  bestBeam = nbkgbeam;
295  bestCuts[metric] = sig->GetXaxis()->GetBinLowEdge(n);
296  fomeff = neff;
297  fompur = npur;
298  }
299  }
300 
301  if(metric == 0) std::cout << "$s/\\sqrt{b}$ opt";
302  if(metric == 1) std::cout << "$s/\\sqrt{s+b}$ opt";
303  if(metric == 2) std::cout << "FA cut";
304 
305  std::cout << std::fixed << std::setprecision(2);
306  std::cout << " &\t " << ((metric == 2) ? faCut : bestCuts[metric])
307  << " &\t " << bestSoB << " &\t " << bestSoSB
308  << " &\t " << bestSig << " &\t " << bestBkg
309  << " &\t " << bestBkg-bestBkgCC-bestBeam
310  << " &\t " << bestBkgCC << " &\t " << bestBeam << " &\t "
311  << " &\t " << fomeff << "\\% &\t " << fompur << "\\%\\\\"<< std::endl;
312 
313  if(metric == 0){
314  fomeff_SoB=fomeff;
315  fompur_SoB=fompur;
316  }
317  if(metric == 1){
318  fomeff_SoSB=fomeff;
319  fompur_SoSB=fompur;
320  }
321 
322  } // end for metric
323 
324  /*
325  new TCanvas;
326  eff->Draw("l");
327  Arrows(eff, bestCuts[0], bestCuts[1], faCut, true);
328  Print("eff", suffix+"_"+eff_suffix);
329 
330  new TCanvas;
331  pur->Draw("l");
332  Arrows(pur, bestCuts[0], bestCuts[1], faCut, true);
333  Print("pur", suffix);
334  */
335 
336 
337  new TCanvas;
338  eff->SetTitle("");
339  eff->GetYaxis()->SetTitle("Percentage");
340  eff->SetLineColor(kRed);
341  pur->SetLineColor(kBlue);
342  effxpur->SetLineColor(kGreen+2);
343  eff->SetMaximum(110);
344  eff->Draw("c");
345  pur->Draw("csame");
346  effxpur->Draw("csame");
347 
348  /*
349 
350  TLegend* legeffpur = new TLegend(.125, .12, .5, .32);
351  legeffpur->AddEntry(eff, "Cumulative Efficiency", "l");
352  legeffpur->AddEntry(pur, "Cumulative Purity", "l");
353  legeffpur->AddEntry(effxpur, "Cumulative Efficiency #times Purity", "l");
354  legeffpur->SetFillStyle(0);
355  legeffpur->Draw();
356  Print("effpur", suffix+"_"+eff_suffix);
357  */
358 
359  new TCanvas;
360  TH2F* hpx = new TH2F("hpx","hpx",100,0,100,100,0,100); // axis range hpx->SetStats(kFALSE); // no statistics
361  hpx->SetStats(kFALSE);
362  hpx->SetTitle("");
363  hpx->GetXaxis()->SetTitle("Efficiency (%)");
364  hpx->GetYaxis()->SetTitle("Purity (%)");
365  hpx->GetXaxis()->CenterTitle();
366  hpx->GetYaxis()->CenterTitle();
367  hpx->Draw();
368  for(double i = .5; i < 10; i += .5){
369  TGraph* c = new TGraph;
370  for(double j = .01; j < 101; j += .1){
371  c->SetPoint(c->GetN(), j, (100*i*i)/j);
372  }
373  c->SetLineWidth(2);
374  c->SetLineColor(kGray+1);
375  c->SetLineStyle(7);
376  c->Draw("l same");
377  }
378 
379 
380  TMarker* m1 = new TMarker(fomeff_SoB, fompur_SoB, kFullCircle);
381  m1->SetMarkerSize(1.5*m1->GetMarkerSize());
382  m1->Draw();
383 
384  TMarker* m2 = new TMarker(fomeff_SoSB, fompur_SoSB, kFullTriangleUp);
385  m2->SetMarkerSize(1.5*m2->GetMarkerSize());
386  m2->Draw();
387 
388  effpurGraph->SetLineWidth(2);
389  effpurGraph->Draw("C");
390  Print("effpurgraph", suffix+"_"+eff_suffix);
391 
392  new TCanvas;
393  TH2F* roch = new TH2F("roch","roch",100,0,1,100,0,1); // axis range roch->SetStats(kFALSE); // no statistics
394  roch->SetStats(kFALSE);
395  roch->SetTitle("");
396  roch->GetXaxis()->SetTitle("False Positive Rate");
397  roch->GetYaxis()->SetTitle("True Positive Rate");
398  roch->GetXaxis()->CenterTitle();
399  roch->GetYaxis()->CenterTitle();
400  roch->Draw();
401  TLine* coin = new TLine(0, 0, 1, 1);
402  coin->SetLineStyle(2);
403  coin->Draw();
404  rocGraph->SetLineWidth(2);
405  rocGraph->Draw("C");
406  Print("rocGraph", suffix+"_"+eff_suffix);
407 
408  TFile* outFile=new TFile((suffix+"_"+eff_suffix+".root").c_str(),"RECREATE");
409  rocGraph->Write("rocGraph");
410  effpurGraph->Write("effpurGraph");
411  outFile->Close();
412 
413  new TCanvas;
414  fomsob->Draw("l");
415  Arrows(fomsob, bestCuts[0], bestCuts[1], faCut, true);
416  Print("fomsob", suffix);
417 
418  new TCanvas;
419  fomsosb->Draw("l");
420  Arrows(fomsosb, bestCuts[0], bestCuts[1], faCut, true);
421  Print("fomsosb", suffix);
422 
423 }
kRed
enum BeamMode kRed
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
PandAna.Demos.pi0_spectra.bkg
bkg
Definition: pi0_spectra.py:116
outFile
TFile * outFile
Definition: PlotXSec.C:135
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
PandAna.Demos.pi0_spectra.pur
tuple pur
Definition: pi0_spectra.py:135
calib::j
const double j
Definition: BetheBloch.cxx:29
nd_projection_maker.eff
list eff
Definition: nd_projection_maker.py:3
om::cout
OStream cout
Definition: OStream.cxx:6
Print
void Print(std::string prefix, std::string name, std::string suffix="")
Definition: nue_pid_effs.C:68
Munits::m2
static constexpr Double_t m2
Definition: Munits.h:145
Arrows
void Arrows(TH1 *sigId, double bestCutSoB, double bestCutSoSB, double faCut, bool noarrow=false)
Definition: nue_pid_effs.C:141
kGreen
enum BeamMode kGreen
kBlue
enum BeamMode kBlue
void GetSpectra ( IPrediction *  pred,
osc::IOscCalc *  calc,
TH1 *&  hSig,
TH1 *&  hNC,
TH1 *&  hCC,
TH1 *&  hBeam,
TH1 *&  hTotBkg 
)

Definition at line 427 of file nue_pid_effs.C.

References ana::Flavors::kAll, ana::Flavors::kAllNuMu, ana::kBeamNueBackgroundColor, ana::Sign::kBoth, ana::Current::kCC, ana::Current::kNC, ana::kNCBackgroundColor, ana::kNueSignalColor, ana::Flavors::kNuEToNuE, ana::kNumuBackgroundColor, ana::Flavors::kNuMuToNuE, kPOTMACRO, ana::kTotalMCColor, ana::IPrediction::PredictComponent(), ana::Spectrum::ToTH1(), and ana::UniqueName().

Referenced by ana::NumuCC2p2hBkgdEstimator::Background(), nue_pid_effs(), PlotEffs(), and PlotSpectra().

429 {
430  hSig = pred->PredictComponent(calc,
431  Flavors::kNuMuToNuE,
432  Current::kCC,
433  Sign::kBoth).ToTH1(kPOTMACRO);
434  hSig->SetLineColor(kNueSignalColor);
435 
436  hNC = pred->PredictComponent(calc,
437  Flavors::kAll,
438  Current::kNC,
439  Sign::kBoth).ToTH1(kPOTMACRO);
440  hNC->SetLineColor(kNCBackgroundColor);
441 
442  hCC = pred->PredictComponent(calc,
443  Flavors::kAllNuMu,
444  Current::kCC,
445  Sign::kBoth).ToTH1(kPOTMACRO);
446  hCC->SetLineColor(kNumuBackgroundColor);
447 
448  hBeam = pred->PredictComponent(calc,
449  Flavors::kNuEToNuE,
450  Current::kCC,
451  Sign::kBoth).ToTH1(kPOTMACRO);
452  hBeam->SetLineColor(kBeamNueBackgroundColor);
453 
454 
455  hTotBkg = (TH1*)hNC->Clone(UniqueName().c_str());
456  hTotBkg->Add(hCC);
457  hTotBkg->Add(hBeam);
458 
459  hTotBkg->SetLineColor(kTotalMCColor);
460 }
ana::IPrediction::PredictComponent
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
ana::Flavors::kNuMuToNuE
(&#39; appearance&#39;)
Definition: IPrediction.h:18
ana::Spectrum::ToTH1
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
ana::Flavors::kNuEToNuE
(&#39;beam &#39;)
Definition: IPrediction.h:15
ana::kNumuBackgroundColor
const Color_t kNumuBackgroundColor
Definition: Style.h:30
ana::kNueSignalColor
const Color_t kNueSignalColor
Definition: Style.h:19
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:39
ana::kBeamNueBackgroundColor
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
ana::kTotalMCColor
const Color_t kTotalMCColor
Definition: Style.h:16
ana::Flavors::kAll
All neutrinos, any flavor.
Definition: IPrediction.h:26
ana::kNCBackgroundColor
const Color_t kNCBackgroundColor
Definition: Style.h:22
ana::UniqueName
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
kPOTMACRO
const double kPOTMACRO
Definition: nue_pid_effs.C:44
void nue_pid_effs ( )

Definition at line 662 of file nue_pid_effs.C.

References calc, om::cout, update_sam_good_runs_metadata::cuts, allTimeWatchdog::endl, FOMPlot(), GetSpectra(), ana::SpectrumLoader::Go(), ISMC, ana::Flavors::kAll, ana::kApplySecondAnalysisMask, ana::Sign::kBoth, ana::Current::kCC, ana::kCVNe, ana::kLEM, ana::kMode, ana::Current::kNC, ana::kNoShift, PandAna.var.analysis_vars::kNueEnergy, ana::kNueSignalColor, ana::Flavors::kNuEToNuE, ana::Flavors::kNuMuToNuE, ana::Flavors::kNuMuToNuMu, kPOTMACRO, ana::kTrueEVis(), ana::kTuftsWeightCC, kXvertex, kY, kYvertex, kZvertex, demo0::loader, demo4::loaderSwap, PIDFileTag(), PIDLongName(), ana::PlotSpectra(), ana::PredictionExtrap::PredictComponent(), ana::IPrediction::PredictComponent(), Print(), ana::Binning::Simple(), std::sqrt(), and ana::Spectrum::ToTH1().

663 {
664 
665  //chose you oscillation parameters:
666  osc::OscCalcDumb calc;
667  // osc::OscCalcPMNSOpt calc;
668  // calc.SetL(810);
669  // calc.SetRho(2.84);
670  // calc.SetDmsq21(7.53e-5);
671  // calc.SetDmsq32(2.44e-3);
672  // calc.SetTh12(asin(sqrt(0.846))/2);
673  // calc.SetTh13(asin(sqrt(.085))/2);
674  // calc.SetTh23(asin(sqrt(0.999))/2);
675  // calc.SetdCP(0);
676 
677  for(int rhc = false; rhc <= true; ++rhc){
678  if(rhc) continue; // Don't exist yet for S13-06-18
679 
680  //SpectrumLoader loader("/pnfs/nova/scratch/users/atsaris/cvnFS_eval_fullInt_afterColl_swap/00014*/*/*.root");
681  //SpectrumLoader loaderSwap("/pnfs/nova/scratch/users/atsaris/cvnFS_eval_fullInt_afterColl_nonSwap/00014*/*/*.root");
682 
683  SpectrumLoader loader("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/semfe_nonS_big.root");
684  SpectrumLoader loaderSwap("/nova/ana/users/atsaris/work.files/cvn_full_chain_test_new/NEW_ERA_NEW_hadd_new/semfe_S_big.root");
685 
686  //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");
687  //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");
688 
689  //how many pids are you optimising cuts on
690  const int kNumPIDs = 3;
691 
692  IPrediction* preds[kNumPIDs];
693 
694  //prediction objects for pid districutions to do the fom optimisation
695  preds[0] = new PredictionNoExtrap(loader, loaderSwap,
696  "LEM", Binning::Simple(100, 0.0, 1.),
697  kLEM, kNueSecondAnaFullPresel && ISMC,kNoShift, kTuftsWeightCC);
698 
699  preds[1] = new PredictionNoExtrap(loader, loaderSwap,
700  "LID", Binning::Simple(110, 0.0, 1.1),
701  kElecID, kNueSecondAnaFullPresel && ISMC,kNoShift, kTuftsWeightCC);
702 
703  preds[2] = new PredictionNoExtrap(loader, loaderSwap,
704  "#nu_{e} classifier output", Binning::Simple(100, 0.0, 1.),
705  kCVNe, kNueSecondAnaFullPresel && ISMC,kNoShift, kTuftsWeightCC);
706 
707  //used as the reference total event count pre-presel, choose wisely
708  PredictionNoExtrap predDenomNoCutRecoE(loader, loaderSwap,
709  "Reconstructed energy (GeV)", kRecoEBinning,
710  kNueEnergy, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment && ISMC,kNoShift, kTuftsWeightCC);
711 
712  IPrediction* predDenomPreselRecoE[kNumPIDs];
713  //prediction objects to make efficiency plots for a pre-chosen pid cut
714  predDenomPreselRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
715  "Reconstructed energy (GeV)",
716  kRecoEBinning,
717  kNueEnergy,
718  kNueSecondAnaFullPresel && ISMC,kNoShift, kTuftsWeightCC);
719  predDenomPreselRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
720  "Reconstructed energy (GeV)",
721  kRecoEBinning,
722  kNueEnergy,
723  kNueSecondAnaFullPresel && ISMC,kNoShift, kTuftsWeightCC);
724 
725 
726  predDenomPreselRecoE[2] = new PredictionNoExtrap(loader, loaderSwap,
727  "Reconstructed energy (GeV)",
728  kRecoEBinning,
729  kNueEnergy,
730  kNueSecondAnaFullPresel && ISMC,kNoShift, kTuftsWeightCC);
731 
732  PredictionNoExtrap predDenomNoCutTrueE(loader, loaderSwap,
733  "True visible energy (GeV)",
734  kRecoEBinning,
735  kTrueEVis, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment && ISMC,kNoShift, kTuftsWeightCC);
736 
737  IPrediction* predDenomPreselTrueE[kNumPIDs];
738  predDenomPreselTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
739  "True visible energy (GeV)",
740  kRecoEBinning,
741  kTrueEVis,
742  kNueSecondAnaFullPresel && ISMC,kNoShift, kTuftsWeightCC);
743  predDenomPreselTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
744  "True visible energy (GeV)",
745  kRecoEBinning,
746  kTrueEVis,
747  kNueSecondAnaFullPresel && ISMC,kNoShift, kTuftsWeightCC);
748 
749  predDenomPreselTrueE[2] = new PredictionNoExtrap(loader, loaderSwap,
750  "True visible energy (GeV)",
751  kRecoEBinning,
752  kTrueEVis,
753  kNueSecondAnaFullPresel && ISMC,kNoShift, kTuftsWeightCC);
754 
755  IPrediction* predsRecoE[kNumPIDs];
756  IPrediction* predsTrueE[kNumPIDs];
757 
758  PredictionNoExtrap predDenomMode(loader, loaderSwap,
759  "", kModeBinning,
760  kMode, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kTuftsWeightCC);
761 
762  IPrediction* predsMode[kNumPIDs];
763 
764  PredictionNoExtrap predDenomY(loader, loaderSwap,
765  "", kYBinning,
766  kY, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kTuftsWeightCC);
767 
768  IPrediction* predsY[kNumPIDs];
769 
770  PredictionNoExtrap predDenomXvertex(loader, loaderSwap,
771  "", kXvertexBinning,
772  kXvertex, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kTuftsWeightCC);
773 
774  IPrediction* predsXvertex[kNumPIDs];
775 
776  PredictionNoExtrap predDenomYvertex(loader, loaderSwap,
777  "", kYvertexBinning,
778  kYvertex, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kTuftsWeightCC);
779 
780  IPrediction* predsYvertex[kNumPIDs];
781 
782  PredictionNoExtrap predDenomZvertex(loader, loaderSwap,
783  "", kZvertexBinning,
784  kZvertex, kApplySecondAnalysisMask && kNueSecondAnaDQ && kNueSecondAnaContainment &&kNueSecondAnaPresel && ISMC,kNoShift, kTuftsWeightCC);
785 
786  IPrediction* predsZvertex[kNumPIDs];
787 
788 
789  predsRecoE[0] = new PredictionNoExtrap(loader, loaderSwap,
790  "Reconstructed energy (GeV)", kRecoEBinning,
791  kNueEnergy, kNueSecondAnaFullPresel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kTuftsWeightCC);
792 
793  predsRecoE[1] = new PredictionNoExtrap(loader, loaderSwap,
794  "Reconstructed energy (GeV)", kRecoEBinning,
795  kNueEnergy, kNueSecondAnaFullPresel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kTuftsWeightCC);
796 
797  predsRecoE[2] = new PredictionNoExtrap(loader, loaderSwap,
798  "Reconstructed energy (GeV)", kRecoEBinning,
799  kNueEnergy, kNueSecondAnaFullPresel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kTuftsWeightCC);
800 
801  predsTrueE[0] = new PredictionNoExtrap(loader, loaderSwap,
802  "True visible energy (GeV)", kRecoEBinning,
803  kTrueEVis, kNueSecondAnaFullPresel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kTuftsWeightCC);
804 
805  predsTrueE[1] = new PredictionNoExtrap(loader, loaderSwap,
806  "True visible energy (GeV)", kRecoEBinning,
807  kTrueEVis, kNueSecondAnaFullPresel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kTuftsWeightCC);
808  predsTrueE[2] = new PredictionNoExtrap(loader, loaderSwap,
809  "True visible energy (GeV)", kRecoEBinning,
810  kTrueEVis, kNueSecondAnaFullPresel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kTuftsWeightCC);
811 
812  predsMode[0] = new PredictionNoExtrap(loader, loaderSwap,
813  "", kModeBinning,
814  kMode, kNueSecondAnaFullPresel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kTuftsWeightCC);
815 
816  predsMode[1] = new PredictionNoExtrap(loader, loaderSwap,
817  "", kModeBinning,
818  kMode, kNueSecondAnaFullPresel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kTuftsWeightCC);
819 
820  predsMode[2] = new PredictionNoExtrap(loader, loaderSwap,
821  "", kModeBinning,
822  kMode, kNueSecondAnaFullPresel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kTuftsWeightCC);
823 
824  predsY[0] = new PredictionNoExtrap(loader, loaderSwap,
825  "True y", kYBinning,
826  kY, kNueSecondAnaFullPresel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kTuftsWeightCC);
827 
828  predsY[1] = new PredictionNoExtrap(loader, loaderSwap,
829  "True y", kYBinning,
830  kY, kNueSecondAnaFullPresel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kTuftsWeightCC);
831 
832 
833  predsY[2] = new PredictionNoExtrap(loader, loaderSwap,
834  "True y", kYBinning,
835  kY, kNueSecondAnaFullPresel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kTuftsWeightCC);
836 
837  predsXvertex[0] = new PredictionNoExtrap(loader, loaderSwap,
838  "X Vertex", kXvertexBinning,
839  kXvertex, kNueSecondAnaFullPresel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kTuftsWeightCC);
840 
841  predsXvertex[1] = new PredictionNoExtrap(loader, loaderSwap,
842  "X Vertex", kXvertexBinning,
843  kXvertex, kNueSecondAnaFullPresel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kTuftsWeightCC);
844 
845 
846  predsXvertex[2] = new PredictionNoExtrap(loader, loaderSwap,
847  "X Vertex", kXvertexBinning,
848  kXvertex, kNueSecondAnaFullPresel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kTuftsWeightCC);
849 
850  predsYvertex[0] = new PredictionNoExtrap(loader, loaderSwap,
851  "Y Vertex", kYvertexBinning,
852  kYvertex, kNueSecondAnaFullPresel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kTuftsWeightCC);
853 
854  predsYvertex[1] = new PredictionNoExtrap(loader, loaderSwap,
855  "Y Vertex", kYvertexBinning,
856  kYvertex, kNueSecondAnaFullPresel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kTuftsWeightCC);
857 
858 
859  predsYvertex[2] = new PredictionNoExtrap(loader, loaderSwap,
860  "Y Vertex", kYvertexBinning,
861  kYvertex, kNueSecondAnaFullPresel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kTuftsWeightCC);
862 
863  predsZvertex[0] = new PredictionNoExtrap(loader, loaderSwap,
864  "Z Vertex", kZvertexBinning,
865  kZvertex, kNueSecondAnaFullPresel && kNueSecondAnaLEMSsb && ISMC,kNoShift, kTuftsWeightCC);
866 
867  predsZvertex[1] = new PredictionNoExtrap(loader, loaderSwap,
868  "Z Vertex", kZvertexBinning,
869  kZvertex, kNueSecondAnaFullPresel && kNueSecondAnaLIDSsb && ISMC,kNoShift, kTuftsWeightCC);
870 
871 
872  predsZvertex[2] = new PredictionNoExtrap(loader, loaderSwap,
873  "Z Vertex", kZvertexBinning,
874  kZvertex, kNueSecondAnaFullPresel && kNueSecondAnaCVNeSsb && ISMC,kNoShift, kTuftsWeightCC);
875 
876  loader.Go();
877  loaderSwap.Go();
878 
879  //loop over pids to find FOM cuts and output performance at those points
880  //for(int pidIdx = 0; pidIdx < kNumPIDs; ++pidIdx){
881  for(int pidIdx = 2; pidIdx < kNumPIDs; ++pidIdx){
882  TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
883  GetSpectra(preds[pidIdx], &calc, hSig, hNC, hCC, hBeam, hTotBkg);
884 
885  double cuts[2];
886  const double faCut = 0.75; //<-- set a cut value that the pid will be evaluated at in addition to the FOM optimization points
887  TH1* sigdenomNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
888  Flavors::kNuMuToNuE,
889  Current::kCC,
890  Sign::kBoth).ToTH1(kPOTMACRO);
891 
892  TH1* backnumuNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
893  Flavors::kNuMuToNuMu,
894  Current::kCC,
895  Sign::kBoth).ToTH1(kPOTMACRO);
896 
897 
898  TH1* backnueNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
899  Flavors::kNuEToNuE,
900  Current::kCC,
901  Sign::kBoth).ToTH1(kPOTMACRO);
902 
903  TH1* backncNoCut = predDenomNoCutRecoE.PredictComponent(&calc,
904  Flavors::kAll,
905  Current::kNC,
906  Sign::kBoth).ToTH1(kPOTMACRO);
907 
908 
909  double totalSig=sigdenomNoCut->Integral(0, -1);
910  double totalBckNuCuCC=backnumuNoCut->Integral(0, -1);
911  double totalBckNuECC=backnueNoCut->Integral(0, -1);
912  double totalBckNC=backncNoCut->Integral(0, -1);
913  double totalBackground=totalBckNuCuCC+totalBckNuECC+totalBckNC;
914  double totalsob=totalSig/sqrt(totalBackground);
915  double totalsosb=totalSig/sqrt(totalSig+totalBackground);
916  double totalefficiency=100;
917  double totalpurity=(totalSig*100)/(totalSig+totalBackground);
918 
919  std::cout << std::fixed << std::setprecision(2);
920  std::cout << " &\t " << totalsob << " &\t " << totalsosb
921  << " &\t " << totalSig << " &\t " << totalBackground
922  << " &\t " << totalBckNC
923  << " &\t " << totalBckNuCuCC << " &\t " << totalBckNuECC
924  << " &\t " << totalefficiency << "\\% &\t " << totalpurity << "\\%\\\\"<< std::endl;
925  // & 2.47 & 2.14 & 18.31 & 54.93 & 18.31 & 18.31 18.31 & 100.00\% & 75.00\%\\
926 
927  sigdenomNoCut->SetLineColor(kNueSignalColor);
928 
929  TH1* sigdenomPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
930  Flavors::kNuMuToNuE,
931  Current::kCC,
932  Sign::kBoth).ToTH1(kPOTMACRO);
933  sigdenomPresel->SetLineColor(kNueSignalColor);
934 
935  TH1* backnumuPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
936  Flavors::kNuMuToNuMu,
937  Current::kCC,
938  Sign::kBoth).ToTH1(kPOTMACRO);
939 
940  TH1* backnuePresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
941  Flavors::kNuEToNuE,
942  Current::kCC,
943  Sign::kBoth).ToTH1(kPOTMACRO);
944 
945  TH1* backncPresel = predDenomPreselRecoE[pidIdx]->PredictComponent(&calc,
946  Flavors::kAll,
947  Current::kNC,
948  Sign::kBoth).ToTH1(kPOTMACRO);
949 
950  double totalPreselSig=sigdenomPresel->Integral(0, -1);
951  double totalPreselBckNuCuCC=backnumuPresel->Integral(0, -1);
952  double totalPreselBckNuECC=backnuePresel->Integral(0, -1);
953  double totalPreselBckNC=backncPresel->Integral(0, -1);
954  double totalPreselBackground=totalPreselBckNuCuCC+totalPreselBckNuECC+totalPreselBckNC;
955  double totalPreselsob=totalPreselSig/sqrt(totalPreselBackground);
956  double totalPreselsosb=totalPreselSig/sqrt(totalPreselSig+totalPreselBackground);
957  double totalPreselefficiency=100;
958  double totalPreselpurity=(totalPreselSig*100)/(totalPreselSig+totalPreselBackground);
959 
960  std::cout << std::fixed << std::setprecision(2);
961  std::cout << " &\t " << totalPreselsob << " &\t " << totalPreselsosb
962  << " &\t " << totalPreselSig << " &\t " << totalPreselBackground
963  << " &\t " << totalPreselBckNC
964  << " &\t " << totalPreselBckNuCuCC << " &\t " << totalPreselBckNuECC
965  << " &\t " << totalPreselefficiency << "\\% &\t " << totalPreselpurity << "\\%\\\\"<< std::endl;
966 
967 
968  backnumuPresel->Add(backnuePresel);
969  backnumuPresel->Add(backncPresel);
970 
971  backnumuNoCut->Add(backnueNoCut);
972  backnumuNoCut->Add(backncNoCut);
973 
974  FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomNoCut, backnumuNoCut, cuts, faCut, PIDFileTag(pidIdx, rhc), "nocut");
975  //FOMPlot(hSig, hTotBkg, hCC, hBeam, sigdenomPresel, backnumuPresel, cuts, faCut, PIDFileTag(pidIdx, rhc), "presel");
976 
977  PlotSpectra(preds[pidIdx], &calc, PIDLongName(pidIdx, rhc),9.8);
978 
979  Print("pid", PIDFileTag(pidIdx, rhc));
980  } // end for pidIdx
981 
982  PlotSpectra(&predDenomNoCutRecoE, &calc, PIDLongName(-1, rhc), 10);
983  Print("recoE", PIDFileTag(-1, rhc));
984  PlotSpectra(&predDenomNoCutTrueE, &calc, PIDLongName(-1, rhc), 10);
985  Print("trueE", PIDFileTag(-1, rhc));
986 
987  } // end for rhc
988 }
kXvertex
const Var kXvertex([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.x;})
ana::kMode
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
kRecoEBinning
const Binning kRecoEBinning
Definition: nue_pid_effs.C:47
ana::kLEM
const Var kLEM
PID
Definition: Vars.cxx:26
ana::IPrediction::PredictComponent
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
FOMPlot
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)
Definition: nue_pid_effs.C:195
ana::kCVNe
const Var kCVNe
PID
Definition: Vars.cxx:35
ana::kApplySecondAnalysisMask
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
kModeBinning
const Binning kModeBinning
Definition: nue_pid_effs.C:48
ana::Flavors::kNuMuToNuE
(&#39; appearance&#39;)
Definition: IPrediction.h:18
ana::Spectrum::ToTH1
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
kYBinning
const Binning kYBinning
Definition: nue_pid_effs.C:49
ana::Flavors::kNuEToNuE
(&#39;beam &#39;)
Definition: IPrediction.h:15
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
kYvertex
const Var kYvertex([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.y;})
ISMC
const Cut ISMC([](const caf::SRProxy *sr){if(sr->hdr.ismc) return true;return false;})
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
kY
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;})
ana::kNueSignalColor
const Color_t kNueSignalColor
Definition: Style.h:19
kZvertexBinning
const Binning kZvertexBinning
Definition: nue_pid_effs.C:52
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:39
osc::OscCalcDumb
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
demo0.loader
loader
Definition: demo0.py:10
kTrueEVis
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;})
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:22
om::cout
OStream cout
Definition: OStream.cxx:6
ana::Flavors::kNuMuToNuMu
(&#39; survival&#39;)
Definition: IPrediction.h:19
Print
void Print(std::string prefix, std::string name, std::string suffix="")
Definition: nue_pid_effs.C:68
GetSpectra
void GetSpectra(IPrediction *pred, osc::IOscCalc *calc, TH1 *&hSig, TH1 *&hNC, TH1 *&hCC, TH1 *&hBeam, TH1 *&hTotBkg)
Definition: nue_pid_effs.C:427
ana::kTuftsWeightCC
const Var kTuftsWeightCC
Definition: XsecTunes.h:30
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
kYvertexBinning
const Binning kYvertexBinning
Definition: nue_pid_effs.C:50
PIDLongName
std::string PIDLongName(int pidIdx, bool rhc)
Definition: nue_pid_effs.C:632
PandAna.var.analysis_vars.kNueEnergy
kNueEnergy
Definition: analysis_vars.py:160
ana::IPrediction
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
update_sam_good_runs_metadata.cuts
string cuts
Definition: update_sam_good_runs_metadata.py:10
ana::PredictionNoExtrap
Prediction that just uses FD MC, with no extrapolation.
Definition: PredictionNoExtrap.h:17
ana::Flavors::kAll
All neutrinos, any flavor.
Definition: IPrediction.h:26
demo4.loaderSwap
loaderSwap
Definition: demo4.py:9
PlotSpectra
void PlotSpectra(IPrediction *pred, osc::IOscCalc *calc, std::string title, double maxy=0)
Definition: nue_pid_effs.C:464
kXvertexBinning
const Binning kXvertexBinning
Definition: nue_pid_effs.C:51
PIDFileTag
std::string PIDFileTag(int pidIdx, bool rhc)
Definition: nue_pid_effs.C:647
kZvertex
const Var kZvertex([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.z;})
kPOTMACRO
const double kPOTMACRO
Definition: nue_pid_effs.C:44
std::string PIDFileTag ( int  pidIdx,
bool  rhc 
)

Definition at line 647 of file nue_pid_effs.C.

References runNovaSAM::ret, and string.

Referenced by nue_pid_effs().

648 {
649  std::string ret;
650  if(pidIdx == -1) ret = "nosel";
651  if(pidIdx == 0) ret = "lem";
652  if(pidIdx == 1) ret = "lid";
653  if(pidIdx == 2) ret = "cvn";
654 
655  if(rhc) ret += "_rhc"; else ret += "_fhc";
656 
657  return ret;
658 }
string
enum BeamMode string
std::string PIDLongName ( int  pidIdx,
bool  rhc 
)

Definition at line 632 of file nue_pid_effs.C.

References runNovaSAM::ret, and string.

Referenced by nue_pid_effs().

633 {
634  std::string ret;
635  if(pidIdx == -1) ret = "No selection";
636  if(pidIdx == 0) ret = "LEM";
637  if(pidIdx == 1) ret = "LID";
638  if(pidIdx == 2) ret = "CVNe";
639 
640  if(rhc) ret += " (RHC)"; else ret += " (FHC)";
641 
642  return ret;
643 }
string
enum BeamMode string
void PlotEffs ( IPrediction *  predNum,
IPrediction *  predDenom,
osc::IOscCalc *  calc,
std::string  title,
std::string  label,
std::string  id,
double  maxy = 0 
)

Definition at line 521 of file nue_pid_effs.C.

References om::cout, allTimeWatchdog::endl, GetSpectra(), maxy, Print(), PrintEffs(), x1, submit_syst::x2, y1, and submit_syst::y2.

524 {
525  TH1 *hSigNum, *hNCNum, *hCCNum, *hBeamNum, *hTotBkgNum;
526  GetSpectra(predNum, calc, hSigNum, hNCNum, hCCNum, hBeamNum, hTotBkgNum);
527 
528  TH1 *hSigDenom, *hNCDenom, *hCCDenom, *hBeamDenom, *hTotBkgDenom;
529  GetSpectra(predDenom, calc, hSigDenom, hNCDenom, hCCDenom, hBeamDenom, hTotBkgDenom);
530 
531  hSigDenom->SetLineStyle(2);
532  hNCDenom->SetLineStyle(2);
533  hBeamDenom->SetLineStyle(2);
534  hTotBkgDenom->SetLineStyle(2);
535  hCCDenom->SetLineStyle(2);
536 
537  new TCanvas();
538 
539  if( label == "mode" ){
540  hNCDenom->GetXaxis()->SetLabelSize(0.06);
541  hNCDenom->GetXaxis()->SetBinLabel(1, "QE");
542  hNCDenom->GetXaxis()->SetBinLabel(2, "Res");
543  hNCDenom->GetXaxis()->SetBinLabel(3, "DIS");
544  hNCDenom->GetXaxis()->SetBinLabel(4, "Coh");
545 
546  hSigNum->GetXaxis()->SetLabelSize(0.06);
547  hSigNum->GetXaxis()->SetBinLabel(1, "QE");
548  hSigNum->GetXaxis()->SetBinLabel(2, "Res");
549  hSigNum->GetXaxis()->SetBinLabel(3, "DIS");
550  hSigNum->GetXaxis()->SetBinLabel(4, "Coh");
551 
552  }
553 
554  hNCDenom->Draw("hist");
555  hCCDenom->Draw("hist same");
556  hBeamDenom->Draw("hist same");
557  hSigDenom->Draw("hist same");
558  Print(label, id, "nocut");
559 
560  new TCanvas();
561  if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
562  //hSigNum->GetYaxis()->SetTitle("Arbitary Units");
563  hSigNum->GetYaxis()->CenterTitle();
564  hSigNum->Draw("hist");
565  hNCNum->Draw("hist same");
566  hBeamNum->Draw("hist same");
567  hCCNum->Draw("hist same");
568  Print(label, id, "cut");
569 
570  new TCanvas;
571 
572  std::cout << "Total efficiency"
573  << " sig: " << 100*hSigNum->Integral(0, -1)/hSigDenom->Integral(0, -1)
574  << " NC: " << 100*hNCNum->Integral(0, -1)/hNCDenom->Integral(0, -1)
575  << " CC: " << 100*hCCNum->Integral(0, -1)/hCCDenom->Integral(0, -1)
576  << " Beam: " << 100*hBeamNum->Integral(0, -1)/hBeamDenom->Integral(0, -1) << std::endl;
577 
578  hSigNum->Divide(hSigDenom);
579  hSigNum->Scale(100);
580  hSigNum->SetTitle("");
581  hSigNum->GetYaxis()->SetTitle("Efficiency (%)");
582  hSigNum->GetYaxis()->SetRangeUser(0, 100);
583  //if(maxy) hSigNum->GetYaxis()->SetRangeUser(0, maxy);
584  hSigNum->GetXaxis()->CenterTitle();
585  hSigNum->GetYaxis()->CenterTitle();
586  hSigNum->Draw("hist");
587 
588  PrintEffs(hSigNum, "$\\nu_e$ sig");
589 
590  hNCNum->Divide(hNCDenom);
591  hNCNum->Scale(100);
592  hNCNum->Draw("hist same");
593 
594  PrintEffs(hNCNum, "NC");
595 
596  hCCNum->Divide(hCCDenom);
597  hCCNum->Scale(100);
598  hCCNum->Draw("hist same");
599 
600  PrintEffs(hCCNum, "$\\nu_\\mu$ CC");
601 
602  hBeamNum->Divide(hBeamDenom);
603  hBeamNum->Scale(100);
604  hBeamNum->Draw("hist same");
605 
606  PrintEffs(hBeamNum, "beam $\\nu_e$");
607 
608  double x1=125;
609  double y1=.6;
610  double x2=.425;
611  double y2=.85;
612 
613  if(title=="y"){
614  x1=0.6;
615  y1=.6;
616  x2=.9;
617  y2=.85;
618  }
619 
620  TLegend* leg2 = new TLegend(.125, .6, .425, .85);
621  leg2->AddEntry(hSigNum, "Signal", "l");
622  leg2->AddEntry(hCCNum, "#nu_{#mu} CC", "l");
623  leg2->AddEntry(hNCNum, "NC", "l");
624  leg2->AddEntry(hBeamNum, "Beam #nu_{e}", "l");
625  leg2->SetFillStyle(0);
626  leg2->Draw();
627 
628 }
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
maxy
double maxy
Definition: plot_datapredictions.C:15
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
submit_syst.x2
x2
Definition: submit_syst.py:250
submit_syst.y2
y2
Definition: submit_syst.py:250
PrintEffs
void PrintEffs(TH1 *h, const char *title)
Definition: nue_pid_effs.C:507
label
const char * label
Definition: cana_time_peak.C:34
om::cout
OStream cout
Definition: OStream.cxx:6
Print
void Print(std::string prefix, std::string name, std::string suffix="")
Definition: nue_pid_effs.C:68
GetSpectra
void GetSpectra(IPrediction *pred, osc::IOscCalc *calc, TH1 *&hSig, TH1 *&hNC, TH1 *&hCC, TH1 *&hBeam, TH1 *&hTotBkg)
Definition: nue_pid_effs.C:427
void PlotSpectra ( IPrediction *  pred,
osc::IOscCalc *  calc,
std::string  title,
double  maxy = 0 
)

Definition at line 464 of file nue_pid_effs.C.

References genie::utils::style::Format(), GetSpectra(), MECModelEnuComparisons::leg, and maxy.

465 {
466  TH1 *hSig, *hNC, *hCC, *hBeam, *hTotBkg;
467  GetSpectra(pred, calc, hSig, hNC, hCC, hBeam, hTotBkg);
468 
469  new TCanvas;
470  hSig->SetTitle("");//title.c_str());
471  double kRealPot=6e20;
472  hSig->GetYaxis()->SetTitle(TString::Format("Events / %.02lf#times10^{20} POT", kRealPot/1e20));
473  //hSig->GetYaxis()->SetTitle("Arbitary Units");
474  hSig->GetXaxis()->CenterTitle();
475  hSig->GetYaxis()->CenterTitle();
476  hSig->Draw("hist");
477  if(maxy) hSig->GetYaxis()->SetRangeUser(0, maxy);
478  hNC->Draw("hist same");
479  hCC->Draw("hist same");
480  hBeam->Draw("hist same");
481  TLegend* leg2 = new TLegend(.125, .55, .5, .85);
482  leg2->AddEntry(hSig, "Appeared #nu_{e}", "l");
483  leg2->AddEntry(hCC, "Survived #nu_{#mu}", "l");
484  leg2->AddEntry(hNC, "NC background", "l");
485  leg2->AddEntry(hBeam, "Beam #nu_{e} background", "l");
486  leg2->SetFillStyle(0);
487  leg2->Draw();
488 
489  static bool once = true;
490  if(once){
491  once = false;
492  TVirtualPad* bak = gPad;
493  new TCanvas;
494  TLegend* leg = new TLegend(.1, .1, .9, .9);
495  leg->AddEntry(hSig, "#nu_{e} signal", "l");
496  leg->AddEntry(hNC, "NC background", "l");
497  leg->AddEntry(hCC, "#nu_{#mu} background", "l");
498  leg->AddEntry(hBeam, "Beam #nu_{e} background", "l");
499  leg->Draw();
500  gPad->Print("plots/SAbless/eff_leg.pdf");
501  gPad->Print("plots/SAbless/eff_leg.eps");
502  gPad->Print("plots/SAbless/eff_leg.png");
503  bak->cd();
504  }
505 }
maxy
double maxy
Definition: plot_datapredictions.C:15
GetSpectra
void GetSpectra(IPrediction *pred, osc::IOscCalc *calc, TH1 *&hSig, TH1 *&hNC, TH1 *&hCC, TH1 *&hBeam, TH1 *&hTotBkg)
Definition: nue_pid_effs.C:427
genie::utils::style::Format
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void Print ( std::string  prefix,
std::string  name,
std::string  suffix = "" 
)

Definition at line 68 of file nue_pid_effs.C.

References caf::Proxy< caf::SRVertexBranch >::elastic, caf::Proxy< caf::StandardRecord >::hdr, ISMC, caf::Proxy< caf::SRHeader >::ismc, ana::kTrueEVis(), kXvertex, kY, kYvertex, kZvertex, caf::Proxy< caf::StandardRecord >::mc, caf::Proxy< caf::SRTruthBranch >::nu, runNovaSAM::ret, Simulation(), sr, tmp, caf::Proxy< caf::SRElastic >::vtx, caf::Proxy< caf::StandardRecord >::vtx, caf::Proxy< caf::SRVector3D >::x, caf::Proxy< caf::SRVector3D >::y, and caf::Proxy< caf::SRVector3D >::z.

Referenced by genie::geometry::FidPolyhedron::clear(), genie::geometry::FidCylinder::FidCylinder(), genie::geometry::FidSphere::FidSphere(), FOMPlot(), genie::geometry::PlaneParam::IsValid(), main(), nue_pid_effs(), genie::flux::operator<<(), plot_ehade_quantbound(), plot_quantbound_twosamples(), plot_quantile_boundaries_2020(), plot_quantile_boundaries_twosamples_2020(), PlotEffs(), and genie::geometry::FidShape::~FidShape().

69 {
70  name = "plots/SAbless/"+prefix+"_"+name;
71  if(!suffix.empty()) name += "_"+suffix;
72  Simulation();
73  gPad->Print((name+".eps").c_str());
74  gPad->Print((name+".pdf").c_str());
75  gPad->Print((name+".png").c_str());
76 }
name
const XML_Char * name
Definition: expat.h:151
Simulation
void Simulation()
Definition: nue_pid_effs.C:57
prefix
const XML_Char * prefix
Definition: expat.h:380
void PrintEffs ( TH1 *  h,
const char *  title 
)

Definition at line 507 of file nue_pid_effs.C.

References om::cout, allTimeWatchdog::endl, MECModelEnuComparisons::i, and plotROC::title.

Referenced by PlotEffs().

508 {
509  std::cout << title;
510  for(int i = 0; i < h->GetNbinsX()+2; ++i){
511  if(h->GetBinCenter(i) > .75 &&
512  h->GetBinCenter(i) < 3.5){
513  std::cout << " &\t " << h->GetBinContent(i);
514  }
515  }
516  std::cout << "\\\\" << std::endl;
517 }
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
om::cout
OStream cout
Definition: OStream.cxx:6
void resetCalc ( osc::IOscCalcAdjustable *  calc)

Definition at line 126 of file nue_pid_effs.C.

References std::asin(), e, osc::_IOscCalcAdjustable< T >::SetdCP(), osc::_IOscCalcAdjustable< T >::SetDmsq21(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), osc::_IOscCalcAdjustable< T >::SetL(), osc::_IOscCalcAdjustable< T >::SetRho(), osc::_IOscCalcAdjustable< T >::SetTh12(), osc::_IOscCalcAdjustable< T >::SetTh13(), osc::_IOscCalcAdjustable< T >::SetTh23(), and std::sqrt().

127 {
128  calc->SetL(810);
129  calc->SetRho(2.75);
130  calc->SetDmsq21(7.6e-5);
131  calc->SetDmsq32(2.35e-3);
132  calc->SetTh12(asin(sqrt(.87))/2); // PDG
133  calc->SetTh13(asin(sqrt(.095))/2);
134  calc->SetTh23(TMath::Pi()/4);
135  calc->SetdCP(0);
136 }
osc::_IOscCalcAdjustable::SetL
virtual void SetL(double L)=0
osc::_IOscCalcAdjustable::SetDmsq21
virtual void SetDmsq21(const T &dmsq21)=0
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
osc::_IOscCalcAdjustable::SetTh13
virtual void SetTh13(const T &th13)=0
osc::_IOscCalcAdjustable::SetDmsq32
virtual void SetDmsq32(const T &dmsq32)=0
osc::_IOscCalcAdjustable::SetRho
virtual void SetRho(double rho)=0
osc::_IOscCalcAdjustable::SetTh23
virtual void SetTh23(const T &th23)=0
e
Float_t e
Definition: plot.C:35
osc::_IOscCalcAdjustable::SetTh12
virtual void SetTh12(const T &th12)=0
osc::_IOscCalcAdjustable::SetdCP
virtual void SetdCP(const T &dCP)=0
std::asin
T asin(T number)
Definition: d0nt_math.hpp:60
void Simulation ( )

Definition at line 57 of file nue_pid_effs.C.

References prelim.

Referenced by Print().

58 {
59  TLatex* prelim = new TLatex(.9, .95, "NOvA Simulation");
60  prelim->SetTextColor(kGray+1);
61  prelim->SetNDC();
62  prelim->SetTextSize(2/30.);
63  prelim->SetTextAlign(32);
64  prelim->Draw();
65 }
prelim
TLatex * prelim
Definition: Xsec_final.C:133

Variable Documentation

string fFileName ="test.txt"

Definition at line 43 of file nue_pid_effs.C.

const Cut ISMC([](const caf::SRProxy *sr){if(sr->hdr.ismc) return true;return false;})

Referenced by nue_pid_effs(), and Print().

const Binning kModeBinning = Binning::Simple(4, -0.5, 3.5)

Definition at line 48 of file nue_pid_effs.C.

const double kPOTMACRO = 6.754e20

Definition at line 44 of file nue_pid_effs.C.

Referenced by GetSpectra(), and nue_pid_effs().

const Binning kRecoEBinning = Binning::Simple(16, 0, 4)

Definition at line 47 of file nue_pid_effs.C.

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;})
const Var kXvertex([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.x;})

Referenced by nue_pid_effs(), and Print().

const Binning kXvertexBinning = Binning::Simple(40, -800, 800)

Definition at line 51 of file nue_pid_effs.C.

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;})

Referenced by mono::Monopole::analyze(), nnbar::RecoAnalysis::analyze(), htk::HoughTrack::analyze(), change_highlighted_cell(), dozooming(), draw_background(), draw_hit(), draw_hits(), htk::Track3D::dzdt(), mono::Track3D::dzdt(), nnbar::NNbarUtilities::energyBalancingVertex(), nnbar::NNbarUtilities::getCoordinateOfHits(), mono::MonopoleTrack::hit_is_on_track(), mouseover(), nue_pid_effs(), mono::operator<<(), ddthelpers::print(), Print(), hough::HoughT::produce(), mono::Track3D::r2_zt(), request_edarea_size(), screen_to_activecell(), screen_to_cell(), screen_to_plane_unbounded(), setboxes(), and mono::Monopole::view().

const Binning kYBinning = Binning::Simple(50, 0, 1)

Definition at line 49 of file nue_pid_effs.C.

const Var kYvertex([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.y;})

Referenced by nue_pid_effs(), and Print().

const Binning kYvertexBinning = Binning::Simple(40, -800, 800)

Definition at line 50 of file nue_pid_effs.C.

const Var kZvertex([](const caf::SRProxy *sr){return sr->vtx.elastic[0].vtx.z;})

Referenced by nue_pid_effs(), and Print().

const Binning kZvertexBinning = Binning::Simple(40, 0, 6000)

Definition at line 52 of file nue_pid_effs.C.

Generated on Mon Mar 8 2021 14:15:25 for NOvA by   doxygen 1.8.11