Evaluate_header.h
Go to the documentation of this file.
1 #include "CAFAna/Cuts/Cuts.h"
2 #include "CAFAna/Cuts/QuantileCuts.h"
3 #include "CAFAna/Cuts/NumuCuts2018.h"
5 
6 #include "CAFAna/Vars/CosmicRejBDTVars.h"
7 #include "CAFAna/Vars/HistAxes.h"
8 
10 
11 #include "Utilities/rootlogon.C"
12 
13 #include "TString.h"
14 #include "TFile.h"
15 
16 namespace ana {
17 
18  // -----------------------------------------------------------
19  // -----------------------------------------------------------
20  void MakeCutVec (std::string sIdentifier, bool isFHC, bool isBPF, bool Quants, Cut &kBaseCut, std::vector<Cut> &Cuts, std::vector<std::string> &CutNames) {
21  std::cerr << "In my function." << std::endl;
22  // --- First off, lets check that my arguments look correctly.
23  // -----------------------------------------------------------
24  if ( sIdentifier == "HighGain" ||
25  (sIdentifier == "Period1" && isFHC) ||
26  (sIdentifier == "Period2" && isFHC)
27  ){
28 
29  } else {
30  // Passed the wrong argument.
31  std::cout << "\n\n ==========================================================="
32  << "\t I don't recognise the Identifier: " << sIdentifier << "."
33  << "\t Please resolve and try again."
34  << std::endl;
35  return;
36  }
37 
38  // --- Next, lets get our Quartile cuts
39  // -----------------------------------------------------------
40  std::string sFHC_RHC = (isFHC == true ? "fhc":"rhc" );
41  std::string fdspecfile = (std::string)std::getenv("NUMUDATA_DIR")+"ana2018/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
42  //std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/official/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
43  TFile* inFile = TFile::Open( fdspecfile.c_str() );
44  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); //TString("FD_full_") + TString(sFHC_RHC) );
45  // How many quantiles?
46  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used
47  // Make my cuts.
48  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
49 
50 
51  // --- Now I want to define my two base cuts.
52  // -----------------------------------------------------------
53  const Cut kBasePlCos = kBaseCut && kNumuCosmicRej2018;
54  // --- Push the cuts into my vector.
55  Cuts.push_back( kBaseCut ); CutNames.push_back( "MyBaseCut" );
56  Cuts.push_back( kBasePlCos ); CutNames.push_back( "DefCosRej" );
57  if ( Quants ) {
58  for ( size_t cc=0; cc < HadEFracQuantCuts.size(); ++cc ) {
59  const Cut QuantCut = HadEFracQuantCuts[cc];
60  Cuts.push_back( kBaseCut && QuantCut ); CutNames.push_back( "MyBaseCut_Quant"+std::to_string(cc+1) );
61  Cuts.push_back( kBasePlCos && QuantCut ); CutNames.push_back( "DefCosRej_Quant"+std::to_string(cc+1) );
62  }
63  }
64 
65  // --- Make a short vector of doubles for my cut values.
66  // -----------------------------------------------------------
67  std::vector<double> MyCutVals_BDT;
68  std::vector<double> MyCutVals_MLP;
69  Var BDTVar = kUnweighted;
70  Var MLPVar = kUnweighted;
71  if ( isBPF ) {
72  if ( isFHC && sIdentifier == "Period1" ) {
73  MyCutVals_BDT.push_back( 0.4795 ); MyCutVals_BDT.push_back( 0.4965 );
74  MyCutVals_MLP.push_back( 0.9661 ); MyCutVals_MLP.push_back( 0.9530 );
75  BDTVar = kBDTCosRej_BPF_FHCPer1;
76  MLPVar = kMLPCosRej_BPF_FHCPer1;
77  } else if ( isFHC && sIdentifier == "Period2" ) {
78  MyCutVals_BDT.push_back( 0.4715 ); MyCutVals_BDT.push_back( 0.4755 );
79  MyCutVals_MLP.push_back( 0.9473 ); MyCutVals_MLP.push_back( 0.9580 );
80  BDTVar = kBDTCosRej_BPF_FHCPer2;
81  MLPVar = kMLPCosRej_BPF_FHCPer2;
82  } else if ( isFHC && sIdentifier == "HighGain" ) {
83  MyCutVals_BDT.push_back( 0.4671 ); MyCutVals_BDT.push_back( 0.4874 );
84  MyCutVals_MLP.push_back( 0.9763 ); MyCutVals_MLP.push_back( 0.9892 );
85  BDTVar = kBDTCosRej_BPF_FHCHigh;
86  MLPVar = kMLPCosRej_BPF_FHCHigh;
87  } else if (!isFHC && sIdentifier == "HighGain" ) {
88  MyCutVals_BDT.push_back( 0.4629 ); MyCutVals_BDT.push_back( 0.4722 );
89  MyCutVals_MLP.push_back( 0.9820 ); MyCutVals_MLP.push_back( 0.9875 );
90  BDTVar = kBDTCosRej_BPF_RHCHigh;
91  MLPVar = kMLPCosRej_BPF_RHCHigh;
92  }
93  } else {
94  if ( isFHC && sIdentifier == "Period1" ) {
95  MyCutVals_BDT.push_back( 0.4793 ); MyCutVals_BDT.push_back( 0.5015 );
96  MyCutVals_MLP.push_back( 0.9398 ); MyCutVals_MLP.push_back( 0.9750 );
97  BDTVar = kBDTCosRej_Kal_FHCPer1;
98  MLPVar = kMLPCosRej_Kal_FHCPer1;
99  } else if ( isFHC && sIdentifier == "Period2" ) {
100  MyCutVals_BDT.push_back( 0.4748 ); MyCutVals_BDT.push_back( 0.4855 );
101  MyCutVals_MLP.push_back( 0.9721 ); MyCutVals_MLP.push_back( 0.9820 );
102  BDTVar = kBDTCosRej_Kal_FHCPer2;
103  MLPVar = kMLPCosRej_Kal_FHCPer2;
104  } else if ( isFHC && sIdentifier == "HighGain" ) {
105  MyCutVals_BDT.push_back( 0.4724 ); MyCutVals_BDT.push_back( 0.4866 );
106  MyCutVals_MLP.push_back( 0.9708 ); MyCutVals_MLP.push_back( 0.9914 );
107  BDTVar = kBDTCosRej_Kal_FHCHigh;
108  MLPVar = kMLPCosRej_Kal_FHCHigh;
109  } else if (!isFHC && sIdentifier == "HighGain" ) {
110  MyCutVals_BDT.push_back( 0.4638 ); MyCutVals_BDT.push_back( 0.4745 );
111  MyCutVals_MLP.push_back( 0.9908 ); MyCutVals_MLP.push_back( 0.9942 );
112  BDTVar = kBDTCosRej_Kal_RHCHigh;
113  MLPVar = kMLPCosRej_Kal_RHCHigh;
114  }
115  }
116 
117  // --- Now make my cut vector!
118  // -----------------------------------------------------------
119  for (int cut=0; cut<2; ++cut) {
120  // --- Set my Cut Vals and Cut Names.
121  double BDTVal = MyCutVals_BDT[cut];
122  double MLPVal = MyCutVals_MLP[cut];
123 
124  std::string BDTCutName = std::to_string( BDTVal );
125  std::string MLPCutName = std::to_string( MLPVal );
126 
127  // -- BDT Cuts
128  Cuts.push_back( kBaseCut && BDTVar > BDTVal ); CutNames.push_back( "BDTCut_"+BDTCutName );
129  Cuts.push_back( kBasePlCos && BDTVar > BDTVal ); CutNames.push_back( "BDTCut_"+BDTCutName+"_CosRej" );
130  if ( Quants ) {
131  for ( size_t cc=0; cc < HadEFracQuantCuts.size(); ++cc ) {
132  const Cut QuantCut = HadEFracQuantCuts[cc];
133  Cuts.push_back( kBaseCut && QuantCut && BDTVar > BDTVal ); CutNames.push_back( "BDTCut_"+BDTCutName+"_Quant"+std::to_string(cc+1) );
134  Cuts.push_back( kBasePlCos && QuantCut && BDTVar > BDTVal ); CutNames.push_back( "BDTCut_"+BDTCutName+"_Quant"+std::to_string(cc+1)+"_CosRej" );
135  }
136  }
137  // -- MLP Cuts
138  Cuts.push_back( kBaseCut && MLPVar > MLPVal ); CutNames.push_back( "MLPCut_"+MLPCutName );
139  Cuts.push_back( kBasePlCos && MLPVar > MLPVal ); CutNames.push_back( "MLPCut_"+MLPCutName+"_CosRej" );
140  if ( Quants ) {
141  for ( size_t cc=0; cc < HadEFracQuantCuts.size(); ++cc ) {
142  const Cut QuantCut = HadEFracQuantCuts[cc];
143  Cuts.push_back( kBaseCut && QuantCut && MLPVar > MLPVal ); CutNames.push_back( "MLPCut_"+MLPCutName+"_Quant"+std::to_string(cc+1) );
144  Cuts.push_back( kBasePlCos && QuantCut && MLPVar > MLPVal ); CutNames.push_back( "MLPCut_"+MLPCutName+"_Quant"+std::to_string(cc+1)+"_CosRej" );
145  }
146  }
147 
148  } // Loop through cuts.
149 
150  return;
151  }
152 
153 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kBDTCosRej_BPF_FHCHigh
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
OStream cerr
Definition: OStream.cxx:7
const Var kBDTCosRej_Kal_RHCHigh
const Cut kNumuCosmicRej2018([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.53 && sr->slc.nhit< 400 && sr->sel.nuecosrej.pngptp< 0.9 );})
Definition: NumuCuts2018.h:19
const Var kBDTCosRej_BPF_FHCPer1
ifstream inFile
Definition: AnaPlotMaker.h:34
const Var kBDTCosRej_BPF_FHCPer2
void MakeCutVec(std::string sIdentifier, bool isFHC, bool isBPF, bool Quants, Cut &kBaseCut, std::vector< Cut > &Cuts, std::vector< std::string > &CutNames)
const Var kBDTCosRej_BPF_RHCHigh
std::string getenv(std::string const &name)
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
static bool isFHC
const Var kBDTCosRej_Kal_FHCPer1
OStream cout
Definition: OStream.cxx:6
const Var kBDTCosRej_Kal_FHCHigh
const Cut cut
Definition: exporter_fd.C:30
void cc()
Definition: test_ana.C:28
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...
bool isBPF
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const Var kBDTCosRej_Kal_FHCPer2
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
std::vector< std::string > CutNames
Definition: MakeCutFlow.C:49
enum BeamMode string