QuantileCuts2020.cxx
Go to the documentation of this file.
3 #include "CAFAna/Core/Binning.h"
9 
10 #include "TFile.h"
11 #include "TH2.h"
12 
13 #include <iostream>
14 
15 namespace ana
16 {
17 
18  //----------------------------------------------------------------------
19  // 2020 Numu EhadFrac quantile cuts
20 
21  std::vector<Cut> GetNumuEhadFracQuantCuts2020( const bool isRHC, const unsigned int nquants )
22  {
23  if ( nquants <= 1 )
24  {
25  std::cout << "Wrong number of EhadFrac quantiles: " << nquants << std::endl;
26  abort();
27  }
28  std::string beam = isRHC ? "rhc" : "fhc";
29  std::string fname = "$NUMUDATA_LIB_PATH/ana2020/Quantiles/quantiles_" + beam + "_full_numu2020.root";
30  std::cout << "\nGetting EhadFrac quantile hist from file:\n" << fname << std::endl;
31  TFile* inFile = TFile::Open( fname.c_str() );
32  TH2* quant_hist = (TH2*)inFile->FindObjectAny( "FDSpec2D" );
33  quant_hist->SetDirectory(0);
34  inFile->Close();
35  delete inFile;
36  return QuantileCutsFromTH2( quant_hist, kNumuCCOptimisedAxis2020, kHadEFracAxis2020, nquants );
37  }
38 
39  //----------------------------------------------------------------------
40  // Axes for constructing the Pt and Theta quantiles
41 
42  const HistAxis kAxisNumuPt( "Reco Muon #left|#vec{p}_{t}#right| (GeV)" , Binning::Simple( 100, 0., 5. ), kNumuMuonPt );
43  const HistAxis kAxisNuePt ( "Reco Elecron #left|#vec{p}_{t}#right| (GeV)", Binning::Simple( 100, 0., 5. ), kNueElectronPt2020 );
44 
45  const HistAxis kAxisNumuTheta( "Reco #theta_{#mu} (deg)", Binning::Simple( 180, 0., 180. ), kRecoThetaMu );
46  const HistAxis kAxisNueTheta ( "Reco #theta_{e} (deg)" , Binning::Simple( 180, 0., 180. ), kFirstProngTheta );
47 
48  const HistAxis kAxisNueBinsNumuEv( "Reco E_{#nu} (GeV)", Binning::Simple( 10, 0., 5. ), kNumuE2020 );
49  const HistAxis kAxisNueBinsNueEv ( "Reco E_{#nu} (GeV)", Binning::Simple( 10, 0., 5. ), kNueEnergy2020 );
50 
51  //----------------------------------------------------------------------
52  // 2020 Numu Pt quantile cuts
53 
54  std::vector<Cut> GetNumuPtQuantCuts2020( const bool isRHC, const unsigned int ehad_quant, const unsigned int nquants )
55  {
56 
57  if ( ehad_quant < 1 || ehad_quant > 4 )
58  {
59  std::cout << "Wrong Ehad quant: " << ehad_quant << std::endl;
60  abort();
61  }
62  if ( nquants <= 1 )
63  {
64  std::cout << "Wrong number of Numu Pt quantiles: " << nquants << std::endl;
65  abort();
66  }
67 
68  std::string beam = isRHC ? "rhc" : "fhc";
69  std::string fname = "$NUMUDATA_LIB_PATH/ana2020/Quantiles/pt_quantile_hists_2020_nd_" + beam + ".root";
70  std::cout << "\nRetrieving Numu Pt quant hists from file:\n" << fname << std::endl;
71  TFile* inFile = TFile::Open( fname.c_str() );
72  if ( inFile->IsZombie() )
73  {
74  std::cout << "Problem with file:\n" << fname << std::endl;
75  abort();
76  }
77 
78  TH2D* quant_hist = (TH2D*)inFile->Get( Form( "RecoPt_vs_RecoEv_NumuBins_Q%d", ehad_quant ) );
79  quant_hist->SetDirectory(0);
80  inFile->Close();
81  delete inFile;
82 
83  std::cout << Form( "\nReturning Numu Pt quant cuts for %d quantiles in EhadFrac quantile %d\n", nquants, ehad_quant ) << std::endl;
84  return QuantileCutsFromTH2( quant_hist, kNumuCCOptimisedAxis2020, kAxisNumuPt, nquants );
85  }
86 
87  //----------------------------------------------------------------------
88  // 2020 Nue Pt quantile cuts
89 
90  std::vector<Cut> GetNueQuantCuts2020( const bool isRHC, const caf::Det_t det, const unsigned int nquants, const ExtrapVar var )
91  {
92 
93  if ( nquants <= 1 )
94  {
95  std::cout << "Wrong number of Nue extrap quantiles: " << nquants << std::endl;
96  abort();
97  }
98 
99  std::string varName;
101  if ( var == kExtrapPt )
102  {
103  varName = "Pt";
104  histName = "RecoPt_vs_RecoEv_NueBins";
105  }
106  else if ( var == kExtrapTheta )
107  {
108  varName = "Theta";
109  histName = "RecoThMu_vs_RecoEv_NueBins";
110  }
111  else
112  {
113  std::cout << "Nue extrap variable must be kExtrapPt or kExtrapTheta" << std::endl;
114  abort();
115  }
116 
117  std::string beam = isRHC ? "rhc" : "fhc";
118  std::string fname = "$NUMUDATA_LIB_PATH/ana2020/Quantiles/pt_quantile_hists_2020_nd_" + beam + ".root";
119  std::cout << "\nRetrieving Nue quantile hists from file:\n" << fname << std::endl;
120  TFile* inFile = TFile::Open( fname.c_str() );
121  if ( inFile->IsZombie() )
122  {
123  std::cout << "Problem with file:\n" << fname << std::endl;
124  abort();
125  }
126 
127  std::cout << histName << std::endl;
128  TH2D* quant_hist = (TH2D*)inFile->Get( histName.c_str() );
129  quant_hist->SetDirectory(0);
130  inFile->Close();
131  delete inFile;
132 
133  std::vector<Cut> quant_cuts;
134 
135  if ( det == caf::kNEARDET )
136  {
137  std::cout << Form( "\nReturning cuts for %d ND Nue %s quantiles\n", nquants, varName.c_str() ) << std::endl;
138  if ( var == kExtrapPt ) quant_cuts = QuantileCutsFromTH2( quant_hist, kAxisNueBinsNumuEv, kAxisNumuPt , nquants );
139  else if ( var == kExtrapTheta ) quant_cuts = QuantileCutsFromTH2( quant_hist, kAxisNueBinsNumuEv, kAxisNumuTheta, nquants );
140  }
141  else if ( det == caf::kFARDET )
142  {
143  std::cout << Form( "\nReturning cuts for %d FD Nue %s quantiles\n", nquants, varName.c_str() ) << std::endl;
144  if ( var == kExtrapPt ) quant_cuts = QuantileCutsFromTH2( quant_hist, kAxisNueBinsNueEv, kAxisNuePt , nquants );
145  else if ( var == kExtrapTheta ) quant_cuts = QuantileCutsFromTH2( quant_hist, kAxisNueBinsNueEv, kAxisNueTheta, nquants );
146  }
147  else
148  {
149  std::cout << "Detector must be caf::kNEARDET or caf::kFARDET" << std::endl;
150  abort();
151  }
152 
153  return quant_cuts;
154 
155  }
156 
157 }
Near Detector underground.
Definition: SREnums.h:10
Det_t
Which NOvA detector?
Definition: SREnums.h:7
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kNumuMuonPt([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){double pmu=sqrt(util::sqr(kMuE(sr))-util::sqr(0.1056));double Zbeam=kCosNumi(sr);double pt=pmu *sqrt(1-Zbeam *Zbeam);return(float) pt;}return-5.f;})
Definition: NumuVars.h:47
const HistAxis kAxisNueBinsNumuEv("Reco E_{#nu} (GeV)", Binning::Simple(10, 0., 5.), kNumuE2020)
const HistAxis kNumuCCOptimisedAxis2020("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kNumuE2020)
Definition: HistAxes.h:25
const Var kNumuE2020
Definition: NumuEFxs.h:217
ifstream inFile
Definition: AnaPlotMaker.h:34
std::vector< Cut > GetNumuPtQuantCuts2020(const bool isRHC, const unsigned int ehad_quant, const unsigned int nquants)
const HistAxis kAxisNuePt("Reco Elecron #left|#vec{p}_{t}#right| (GeV)", Binning::Simple(100, 0., 5.), kNueElectronPt2020)
const HistAxis kHadEFracAxis2020("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kNumuHadEFrac2020)
Definition: HistAxes.h:31
const Var kFirstProngTheta([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1000.f;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1000.f;float ct=kFirstProngCosTheta(sr);if(fabs(ct) > 1) return-1000.f;return float(acos(ct)*180.0/3.1416);})
Definition: NueVarsExtra.h:63
std::vector< Cut > GetNueQuantCuts2020(const bool isRHC, const caf::Det_t det, const unsigned int nquants, const ExtrapVar var)
const HistAxis kAxisNumuPt("Reco Muon #left|#vec{p}_{t}#right| (GeV)", Binning::Simple(100, 0., 5.), kNumuMuonPt)
OStream cout
Definition: OStream.cxx:6
const HistAxis kAxisNueBinsNueEv("Reco E_{#nu} (GeV)", Binning::Simple(10, 0., 5.), kNueEnergy2020)
const Var kNueEnergy2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;if(kIsRHC(sr)) return kNueEnergyRHC2020(sr);else return kNueEnergyFHC2020(sr);})
Nue energy with 3d prong info. only (old version of vars)
Definition: NueEnergy2020.h:38
const Var kNueElectronPt2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1000.f;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1000.f;TVector3 elecDir=(TVector3) sr->vtx.elastic.fuzzyk.png[0].dir;TVector3 beamDir=ana::NuMIBeamDirection(sr->hdr.det);float ct=elecDir.Dot(beamDir);float p=kNueElectronE2020(sr);if(p< 0||fabs(ct) > 1) return-1000.f;return float(p *sin(acos(ct)));})
Definition: NueEnergy2020.h:43
const Var kRecoThetaMu([](const caf::SRProxy *sr){return acos(kCosNumi(sr))*180.0/3.1416;})
Definition: NumuVars.h:44
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
const HistAxis kAxisNueTheta("Reco #theta_{e} (deg)", Binning::Simple(180, 0., 180.), kFirstProngTheta)
const HistAxis kAxisNumuTheta("Reco #theta_{#mu} (deg)", Binning::Simple(180, 0., 180.), kRecoThetaMu)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
std::vector< Cut > GetNumuEhadFracQuantCuts2020(const bool isRHC, const unsigned int nquants)
enum BeamMode string