Evaluate_BDTMLP_Algs_PredNoExtrap.C
Go to the documentation of this file.
1 // A macro to test the results that I got from training the Cos Rej BDT.
2 // something like TMVA.
3 
4 #include "CAFAna/Core/Binning.h"
5 #include "CAFAna/Core/Spectrum.h"
6 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Cuts/QuantileCuts.h"
8 
9 #include "CAFAna/Vars/BPFVars.h"
11 #include "CAFAna/Vars/CosmicRejBDTVars.h"
13 #include "CAFAna/Vars/HistAxes.h"
15 #include "CAFAna/Vars/Vars.h"
16 #include "CAFAna/Vars/NumuVars.h"
17 #include "CAFAna/Vars/XsecTunes.h"
18 
20 
22 
23 #include "Utilities/rootlogon.C"
24 
25 #include "TTree.h"
26 #include "TString.h"
27 
28 using namespace ana;
29 
30 // -----------------------------------------------------------
31 // -----------------------------------------------------------
32 // --- Now for my actual module.
33 // -----------------------------------------------------------
34 // -----------------------------------------------------------
35 void Evaluate_BDTMLP_Algs_PredNoExtrap( std::string sIdentifier, bool isFHC, bool isBPF, bool Quant, bool isEval = true ) {
36  // -----------------------------------------------------------
37  // --- First off, lets set some global things and make sure we know what we're running over...
38  // -----------------------------------------------------------
39  std::string sFHC_RHC = (isFHC == true ? "fhc" :"rhc" );
40  std::string sisEval = (isEval == true ? "Test" :"Train" );
41  std::string sisBPF = (isBPF == true ? "BPF" :"Kal" );
42  std::string sisData = "BeamMC";
43 
44  if (!isEval) {
45  std::cout << "\n\nYou're running over the training datasets....Are you sure you want to do this?"
46  << "\nThey should be reserved for training, not evaluating!"
47  << std::endl;
48  }
49 
50  std::string OutName = "BDTSpectra-"+sisData+"-"+sFHC_RHC+"-"+sisBPF+"-"+sIdentifier+"-"+sisEval+".root";
51 
52  // And write them out!
53  std::cout << "\n\nLots of input parameters, so I'm going to write them out now....;"
54  << "\n\t isFHC " << (isFHC == true ? "true --> Run in FHC " :"false --> Run in RHC" ) << " ==> sFHC_RHC is " << sFHC_RHC
55  << "\n\t isBPF " << (isBPF == true ? "true --> BPF training" :"false --> Kal training" ) << " ==> sisBPF is " << sisBPF
56  << "\n\t Making PredNoExtrap, therefore can only run over MC."
57  << "\n\t Identifier is " << sIdentifier
58  << "\n\t isEval " << (isEval == true ? "true --> Evaluation dataset":"false --> Training dataset" ) << " ==> sisEval is " << sisEval
59  << "\n\n OutName is therefore " << OutName
60  << "\n\n" << std::endl;
61 
62  // -----------------------------------------------------------
63  // --- Now lets get our Quartile cuts
64  // -----------------------------------------------------------
65  /*
66  std::string fdspecfile = (std::string)std::getenv("NUMUDATA_DIR")+"ana2018/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
67  //std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/official/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
68  TFile* inFile = TFile::Open( fdspecfile.c_str() );
69  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); //TString("FD_full_") + TString(sFHC_RHC) );
70  // How many quantiles?
71  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used
72  // Make my cuts.
73  //std::vector<Cut> HadEFracQuantCuts = QuantileAndPIDCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles, isFHC);
74  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
75  */
76 
77  // -----------------------------------------------------------
78  // --- Lets figure out which definition I want to run over....
79  // -----------------------------------------------------------
80  std::string NonSwap = "karlwarb-prod4reco-cosmicrej-fd-genie-nonswap-" +sFHC_RHC+"-"+sIdentifier+"_"+sisEval;
81  std::string FluxSwap = "karlwarb-prod4reco-cosmicrej-fd-genie-fluxswap-"+sFHC_RHC+"-"+sIdentifier+"_"+sisEval;
82  std::string TauSwap = "karlwarb-prod4reco-cosmicrej-fd-genie-tau-" +sFHC_RHC+"-"+sIdentifier+"_"+sisEval;
83 
84  // --- And declare my loader.
85  //SpectrumLoader loader( MyDef );
86  auto loader = new Loaders();
87  loader -> SetLoaderPath( NonSwap , caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap );
88  loader -> SetLoaderPath( FluxSwap, caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap );
89  loader -> SetLoaderPath( TauSwap , caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kTauSwap );
90  loader -> SetSpillCut(kStandardSpillCuts);
91 
92  // -----------------------------------------------------------
93  // --- What weights do I want to use?
94  // -----------------------------------------------------------
96 
97  // -----------------------------------------------------------
98  // --- What cuts do I want to use?
99  // -----------------------------------------------------------
101  std::vector<Cut> Cuts;
102  std::vector<std::string> CutNames;
103  std::cout << "Going into my function." << std::endl;
104  MakeCutVec( sIdentifier, isFHC, isBPF, Quant, BaseCut, Cuts, CutNames );
105  unsigned int NumCuts = Cuts.size();
106  std::cout << "Out of my function. CutNames has size " << NumCuts << std::endl;
107 
108  // -----------------------------------------------------------
109  // --- Any custom Hist Axes?
110  // -----------------------------------------------------------
111  const unsigned int NumBins = 100;
112  const Binning NumHitBins = Binning::Simple(NumBins, 0 , 1000 );
113  const Binning LengthBins = Binning::Simple(NumBins, 0 , 20 );
114  const Binning RatioBins = Binning::Simple(NumBins, -1.01, 1.01 );
115  const Binning XRangeBins = Binning::Simple(NumBins, -8 , 8 );
116  const Binning YRangeBins = Binning::Simple(NumBins, -8 , 8 );
117  const Binning ZRangeBins = Binning::Simple(NumBins, -0.1 , 60.1 );
118 
119  // -----------------------------------------------------------
120  // --- Define my spectra.
121  // -----------------------------------------------------------
122  // Energy quanitities
123  //PredictionNoExtrap *sTrueNuE [NumCuts];
124  PredictionNoExtrap *sRecoNuE [NumCuts];
125  //PredictionNoExtrap *sNuMuPtOP [NumCuts];
126  /*
127  // --- My New Variables.
128  PredictionNoExtrap *sCosRejPID;
129  PredictionNoExtrap *sFHCMLPPID;
130  PredictionNoExtrap *sRHCMLPPID;
131  PredictionNoExtrap *sFHCBDTPID;
132  PredictionNoExtrap *sRHCBDTPID;
133  // -- Some 2D spectra vs Reco E and True E
134  // Reco
135  PredictionNoExtrap *sRecoE_CosRejPID;
136  PredictionNoExtrap *sRecoE_FHCMLPPID;
137  PredictionNoExtrap *sRecoE_RHCMLPPID;
138  PredictionNoExtrap *sRecoE_FHCBDTPID;
139  PredictionNoExtrap *sRecoE_RHCBDTPID;
140  // True
141  PredictionNoExtrap *sTrueE_CosRejPID;
142  PredictionNoExtrap *sTrueE_FHCMLPPID;
143  PredictionNoExtrap *sTrueE_RHCMLPPID;
144  PredictionNoExtrap *sTrueE_FHCBDTPID;
145  PredictionNoExtrap *sTrueE_RHCBDTPID;
146  */
147  // -----------------------------------------------------------
148  // --- Create the spectra classes
149  // -----------------------------------------------------------
150  for (unsigned int cc=0; cc<NumCuts; ++cc) {
151  Cut ThisCut = Cuts[cc];
152  // Energy quanitities
153  //sTrueNuE [cc] = new PredictionNoExtrap( *loader, "True Neutrino E", kNumuCCEOptimisedBinning, kTrueE , ThisCut, kNoShift, kEvtWgt );
154  sRecoNuE [cc] = new PredictionNoExtrap( *loader, "Reco Neutrino E", kNumuCCEOptimisedBinning, kCCE , ThisCut, kNoShift, kEvtWgt );
155  //sNuMuPtOP [cc] = new PredictionNoExtrap( *loader, "Muon Pt Over P" , RatioBins , kNumuMuonPtP, ThisCut, kNoShift, kEvtWgt );
156  }
157  /*
158  // --- My New Variables.
159  sCosRejPID = new PredictionNoExtrap( *loader, "Ana18 CosRej" , RatioBins, kNumuContPID , BaseCut, kNoShift, kEvtWgt );
160  sFHCMLPPID = new PredictionNoExtrap( *loader, "FHC MLP CosRej", RatioBins, kFHCMLPCosRej, BaseCut, kNoShift, kEvtWgt );
161  sRHCMLPPID = new PredictionNoExtrap( *loader, "RHC MLP CosRej", RatioBins, kRHCMLPCosRej, BaseCut, kNoShift, kEvtWgt );
162  sFHCBDTPID = new PredictionNoExtrap( *loader, "FHC BDT CosRej", RatioBins, kFHCBDTCosRej, BaseCut, kNoShift, kEvtWgt );
163  sRHCBDTPID = new PredictionNoExtrap( *loader, "RHC BDT CosRej", RatioBins, kRHCBDTCosRej, BaseCut, kNoShift, kEvtWgt );
164  // -- Some 2D spectra vs Reco E and True E
165  // Reco
166  sRecoE_CosRejPID = new PredictionNoExtrap( *loader, "RecoE vs Ana18 CosRej" , kNumuCCEOptimisedBinning, kCCE , RatioBins, kNumuContPID , BaseCut, kNoShift, kEvtWgt );
167  sRecoE_FHCMLPPID = new PredictionNoExtrap( *loader, "RecoE vs FHC MLP CosRej", kNumuCCEOptimisedBinning, kCCE , RatioBins, kFHCMLPCosRej, BaseCut, kNoShift, kEvtWgt );
168  sRecoE_RHCMLPPID = new PredictionNoExtrap( *loader, "RecoE vs RHC MLP CosRej", kNumuCCEOptimisedBinning, kCCE , RatioBins, kRHCMLPCosRej, BaseCut, kNoShift, kEvtWgt );
169  sRecoE_FHCBDTPID = new PredictionNoExtrap( *loader, "RecoE vs FHC BDT CosRej", kNumuCCEOptimisedBinning, kCCE , RatioBins, kFHCBDTCosRej, BaseCut, kNoShift, kEvtWgt );
170  sRecoE_RHCBDTPID = new PredictionNoExtrap( *loader, "RecoE vs RHC BDT CosRej", kNumuCCEOptimisedBinning, kCCE , RatioBins, kRHCBDTCosRej, BaseCut, kNoShift, kEvtWgt );
171  // True
172  sTrueE_CosRejPID = new PredictionNoExtrap( *loader, "TrueE vs Ana18 CosRej" , kNumuCCEOptimisedBinning, kTrueE, RatioBins, kNumuContPID , BaseCut, kNoShift, kEvtWgt );
173  sTrueE_FHCMLPPID = new PredictionNoExtrap( *loader, "TrueE vs FHC MLP CosRej", kNumuCCEOptimisedBinning, kTrueE, RatioBins, kFHCMLPCosRej, BaseCut, kNoShift, kEvtWgt );
174  sTrueE_RHCMLPPID = new PredictionNoExtrap( *loader, "TrueE vs RHC MLP CosRej", kNumuCCEOptimisedBinning, kTrueE, RatioBins, kRHCMLPCosRej, BaseCut, kNoShift, kEvtWgt );
175  sTrueE_FHCBDTPID = new PredictionNoExtrap( *loader, "TrueE vs FHC BDT CosRej", kNumuCCEOptimisedBinning, kTrueE, RatioBins, kFHCBDTCosRej, BaseCut, kNoShift, kEvtWgt );
176  sTrueE_RHCBDTPID = new PredictionNoExtrap( *loader, "TrueE vs RHC BDT CosRej", kNumuCCEOptimisedBinning, kTrueE, RatioBins, kRHCBDTCosRej, BaseCut, kNoShift, kEvtWgt );
177  */
178  // --- Do it!
179  loader -> Go();
180 
181  // -----------------------------------------------------------
182  // --- Finally, write everything to my output file.
183  // -----------------------------------------------------------
184  TFile *outFile = new TFile( OutName.c_str(), "RECREATE" );
185  outFile->cd();
186  for (unsigned int cc=0; cc<NumCuts; ++cc) {
187  // Energy quanitities
188  //sTrueNuE [cc] -> SaveTo( outFile->mkdir( TString("TrueNuE" )+TString(CutNames[cc]) ) );
189  sRecoNuE [cc] -> SaveTo( outFile->mkdir( TString("RecoNuE_" )+TString(CutNames[cc]) ) );
190  //sNuMuPtOP [cc] -> SaveTo( outFile->mkdir( TString("MuPtOverP")+TString(CutNames[cc]) ) );
191  }
192  /*
193  // --- My New Variables.
194  sCosRejPID -> SaveTo( outFile->mkdir( TString("Ana18CosRej" ) ) );
195  sFHCMLPPID -> SaveTo( outFile->mkdir( TString("FHCMLPCosRej") ) );
196  sRHCMLPPID -> SaveTo( outFile->mkdir( TString("RHCMLPCosRej") ) );
197  sFHCBDTPID -> SaveTo( outFile->mkdir( TString("FHCBDTCosRej") ) );
198  sRHCBDTPID -> SaveTo( outFile->mkdir( TString("RHCBDTCosRej") ) );
199  // -- Some 2D spectra vs Reco E and True E
200  // Reco
201  sRecoE_CosRejPID -> SaveTo( outFile->mkdir( TString("RecoE_vs_Ana18CosRej" ) ) );
202  sRecoE_FHCMLPPID -> SaveTo( outFile->mkdir( TString("RecoE_vs_FHCMLPCosRej") ) );
203  sRecoE_RHCMLPPID -> SaveTo( outFile->mkdir( TString("RecoE_vs_RHCMLPCosRej") ) );
204  sRecoE_FHCBDTPID -> SaveTo( outFile->mkdir( TString("RecoE_vs_FHCBDTCosRej") ) );
205  sRecoE_RHCBDTPID -> SaveTo( outFile->mkdir( TString("RecoE_vs_RHCBDTCosRej") ) );
206  // True
207  sTrueE_CosRejPID -> SaveTo( outFile->mkdir( TString("TrueE_vs_Ana18CosRej" ) ) );
208  sTrueE_FHCMLPPID -> SaveTo( outFile->mkdir( TString("TrueE_vs_FHCMLPCosRej") ) );
209  sTrueE_RHCMLPPID -> SaveTo( outFile->mkdir( TString("TrueE_vs_RHCMLPCosRej") ) );
210  sTrueE_FHCBDTPID -> SaveTo( outFile->mkdir( TString("TrueE_vs_FHCBDTCosRej") ) );
211  sTrueE_RHCBDTPID -> SaveTo( outFile->mkdir( TString("TrueE_vs_RHCBDTCosRej") ) );
212  */
213  // --- And now close and we'll be done!
214  outFile->Close();
215 
216 } // --- DONE! ---
const Binning ZRangeBins
Far Detector at Ash River.
Definition: SREnums.h:11
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
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
const Binning YRangeBins
TFile * outFile
Definition: PlotXSec.C:135
const Cut kNumuContainFD2017
Definition: NumuCuts2017.h:21
void MakeCutVec(std::string sIdentifier, bool isFHC, bool isBPF, bool Quants, Cut &kBaseCut, std::vector< Cut > &Cuts, std::vector< std::string > &CutNames)
const Var kCCE
Definition: NumuVars.h:21
loader
Definition: demo0.py:10
void Evaluate_BDTMLP_Algs_PredNoExtrap(std::string sIdentifier, bool isFHC, bool isBPF, bool Quant, bool isEval=true)
static bool isFHC
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const Binning NumHitBins
const Binning LengthBins
void cc()
Definition: test_ana.C:28
const Cut kNumuPID2018([](const caf::SRProxy *sr){std::cout<< "ERROR::kNumuPID2018, cutting on both cvnProd3Train and cvn2017."<< " Neither branch exists anymore. Returning False."<< std::endl;abort();return false;})
Definition: NumuCuts2018.h:22
const Binning kNumuCCEOptimisedBinning
Optimised binning for numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39; in that talk...
Definition: Binnings.cxx:28
bool isBPF
const unsigned int NumBins
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Prediction that just uses FD MC, with no extrapolation.
const Cut kNumuQuality
Definition: NumuCuts.h:18
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kXSecCVWgt2018RPAFix_noDIS
Definition: XsecTunes.h:57
const Binning RatioBins
std::vector< std::string > CutNames
Definition: MakeCutFlow.C:49
enum BeamMode string