CalcCutVals.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/Loaders.h"
6 #include "CAFAna/Core/Spectrum.h"
8 #include "CAFAna/Cuts/QuantileCuts.h"
9 
10 #include "NuXAna/Cuts/NusCuts18.h"
11 
12 #include "CAFAna/Vars/BPFVars.h"
14 #include "CAFAna/Vars/CosmicRejBDTVars.h"
16 #include "CAFAna/Vars/HistAxes.h"
18 #include "CAFAna/Vars/Vars.h"
19 #include "CAFAna/Vars/NumuVars.h"
20 #include "CAFAna/Vars/XsecTunes.h"
21 
23 
24 #include "CAFAna/Analysis/Calcs.h"
26 #include "OscLib/OscCalcPMNSOpt.h"
27 
28 #include "Utilities/rootlogon.C"
29 
30 #include "TTree.h"
31 #include "TString.h"
32 
33 using namespace ana;
34 
35 // -----------------------------------------------------------
36 // -----------------------------------------------------------
37 // --- Now for my actual module.
38 // -----------------------------------------------------------
39 // -----------------------------------------------------------
40 void CalcCutVals( bool isFHC, std::string sIdentifier ) {
41  // -----------------------------------------------------------
42  // --- First off, lets set some global things and make sure we know what we're running over...
43  // -----------------------------------------------------------
44  std::string sFHC_RHC = "fhc";
45  if (!isFHC) sFHC_RHC = "rhc";
46 
47  std::string OutName = "CalcCutVals-"+sFHC_RHC+"-"+sIdentifier+".root";
48 
49  // And write them out!
50  std::cout << "\n\nLots of input parameters, so I'm going to write them out now....;"
51  << "\n\t isFHC " << (isFHC == true ? "true --> Run in FHC " :"false --> Run in RHC" ) << " ==> sFHC_RHC is " << sFHC_RHC
52  << "\n\t sIdentifier is " << sIdentifier
53  << "\n\n OutName is therefore " << OutName
54  << "\n\n" << std::endl;
55 
56  // -----------------------------------------------------------
57  // --- Now lets get our Quartile cuts
58  // -----------------------------------------------------------
59  /*
60  std::string fdspecfile = (std::string)std::getenv("NUMUDATA_DIR")+"ana2018/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
61  TFile* inFile = TFile::Open( fdspecfile.c_str() );
62  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); //TString("FD_full_") + TString(sFHC_RHC) );
63  // How many quantiles?
64  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used
65  // Make my cuts.
66  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
67  */
68 
69  // -----------------------------------------------------------
70  // --- Lets figure out which definition I want to run over....
71  // -----------------------------------------------------------
72  // ---- Cosmic Data.
73  std::string MyDataDef = "karlwarb-prod4reco-cosmicrej-fd-cosmics-loose-"+sFHC_RHC+"-"+sIdentifier;
74  SpectrumLoader DataLoader( MyDataDef );
75  DataLoader.SetSpillCut(kStandardSpillCuts);
76 
77  // ---- Monte Carlo.
78  std::string NonSwap = "karlwarb-prod4reco-cosmicrej-fd-genie-nonswap-" +sFHC_RHC+"-"+sIdentifier+"_Test";
79  std::string FluxSwap = "karlwarb-prod4reco-cosmicrej-fd-genie-fluxswap-"+sFHC_RHC+"-"+sIdentifier+"_Test";
80  std::string TauSwap = "karlwarb-prod4reco-cosmicrej-fd-genie-tau-" +sFHC_RHC+"-"+sIdentifier+"_Test";
81 
82  // --- And declare my loader.
83  auto MontLoader = new Loaders();
84  MontLoader -> SetLoaderPath( NonSwap , caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap );
85  MontLoader -> SetLoaderPath( FluxSwap, caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap );
86  MontLoader -> SetLoaderPath( TauSwap , caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kTauSwap );
87  MontLoader -> SetSpillCut(kStandardSpillCuts);
88 
89  // -----------------------------------------------------------
90  // --- What cuts and weights do I want to use?
91  // -----------------------------------------------------------
92  const Cut MyCut = kNumuQuality && kNumuContainFD2017 && kNumuPID2018;
94 
95  // -----------------------------------------------------------
96  // --- Any custom Hist Axes?
97  // -----------------------------------------------------------
98  const Binning RatioBins = Binning::Simple( 1000, 0, 1 );
99 
100  // -----------------------------------------------------------
101  // --- Define my spectra.
102  // -----------------------------------------------------------
103  // ================
104  // = Data Spectra =
105  // ================
106  // --- BPF spectra.
107  //*
108  Spectrum *sData_FHC_BPF_MLP_Per1;
109  Spectrum *sData_FHC_BPF_BDT_Per1;
110  Spectrum *sData_FHC_BPF_MLP_Per2;
111  Spectrum *sData_FHC_BPF_BDT_Per2;
112  Spectrum *sData_FHC_BPF_MLP_High;
113  Spectrum *sData_FHC_BPF_BDT_High;
114  Spectrum *sData_RHC_BPF_MLP_High;
115  Spectrum *sData_RHC_BPF_BDT_High;
116 
117  // --- Kalman spectra.
118  Spectrum *sData_FHC_Kal_MLP_Per1;
119  Spectrum *sData_FHC_Kal_BDT_Per1;
120  Spectrum *sData_FHC_Kal_MLP_Per2;
121  Spectrum *sData_FHC_Kal_BDT_Per2;
122  Spectrum *sData_FHC_Kal_MLP_High;
123  Spectrum *sData_FHC_Kal_BDT_High;
124  Spectrum *sData_RHC_Kal_MLP_High;
125  Spectrum *sData_RHC_Kal_BDT_High;
126  //*/
127  // ================
128  // = Mont Spectra =
129  // ================
130  // --- BPF spectra.
131  PredictionNoExtrap *sMont_FHC_BPF_MLP_Per1;
132  PredictionNoExtrap *sMont_FHC_BPF_BDT_Per1;
133  PredictionNoExtrap *sMont_FHC_BPF_MLP_Per2;
134  PredictionNoExtrap *sMont_FHC_BPF_BDT_Per2;
135  PredictionNoExtrap *sMont_FHC_BPF_MLP_High;
136  PredictionNoExtrap *sMont_FHC_BPF_BDT_High;
137  PredictionNoExtrap *sMont_RHC_BPF_MLP_High;
138  PredictionNoExtrap *sMont_RHC_BPF_BDT_High;
139 
140  // --- Kalman spectra.
141  PredictionNoExtrap *sMont_FHC_Kal_MLP_Per1;
142  PredictionNoExtrap *sMont_FHC_Kal_BDT_Per1;
143  PredictionNoExtrap *sMont_FHC_Kal_MLP_Per2;
144  PredictionNoExtrap *sMont_FHC_Kal_BDT_Per2;
145  PredictionNoExtrap *sMont_FHC_Kal_MLP_High;
146  PredictionNoExtrap *sMont_FHC_Kal_BDT_High;
147  PredictionNoExtrap *sMont_RHC_Kal_MLP_High;
148  PredictionNoExtrap *sMont_RHC_Kal_BDT_High;
149 
150  // -----------------------------------------------------------
151  // --- Create the spectra classes
152  // -----------------------------------------------------------
153  Cut DataCut = MyCut && kInCosmicTimingWindow;
154  Cut MontCut = MyCut;// && kInBeamSpill;
155 
156  // ================
157  // = Data Spectra =
158  // ================
159  // --- BPF spectra.
160  sData_FHC_BPF_MLP_Per1 = new Spectrum( "Data_FHC_BPF_MLP_Per1", RatioBins, DataLoader, kMLPCosRej_BPF_FHCPer1, DataCut );
161  sData_FHC_BPF_BDT_Per1 = new Spectrum( "Data_FHC_BPF_BDT_Per1", RatioBins, DataLoader, kBDTCosRej_BPF_FHCPer1, DataCut );
162  sData_FHC_BPF_MLP_Per2 = new Spectrum( "Data_FHC_BPF_MLP_Per2", RatioBins, DataLoader, kMLPCosRej_BPF_FHCPer2, DataCut );
163  sData_FHC_BPF_BDT_Per2 = new Spectrum( "Data_FHC_BPF_BDT_Per2", RatioBins, DataLoader, kBDTCosRej_BPF_FHCPer2, DataCut );
164  sData_FHC_BPF_MLP_High = new Spectrum( "Data_FHC_BPF_MLP_High", RatioBins, DataLoader, kMLPCosRej_BPF_FHCHigh, DataCut );
165  sData_FHC_BPF_BDT_High = new Spectrum( "Data_FHC_BPF_BDT_High", RatioBins, DataLoader, kBDTCosRej_BPF_FHCHigh, DataCut );
166  sData_RHC_BPF_MLP_High = new Spectrum( "Data_RHC_BPF_MLP_High", RatioBins, DataLoader, kMLPCosRej_BPF_RHCHigh, DataCut );
167  sData_RHC_BPF_BDT_High = new Spectrum( "Data_RHC_BPF_BDT_High", RatioBins, DataLoader, kBDTCosRej_BPF_RHCHigh, DataCut );
168 
169  // --- Kalman spectra.
170  sData_FHC_Kal_MLP_Per1 = new Spectrum( "Data_FHC_Kal_MLP_Per1", RatioBins, DataLoader, kMLPCosRej_Kal_FHCPer1, DataCut );
171  sData_FHC_Kal_BDT_Per1 = new Spectrum( "Data_FHC_Kal_BDT_Per1", RatioBins, DataLoader, kBDTCosRej_Kal_FHCPer1, DataCut );
172  sData_FHC_Kal_MLP_Per2 = new Spectrum( "Data_FHC_Kal_MLP_Per2", RatioBins, DataLoader, kMLPCosRej_Kal_FHCPer2, DataCut );
173  sData_FHC_Kal_BDT_Per2 = new Spectrum( "Data_FHC_Kal_BDT_Per2", RatioBins, DataLoader, kBDTCosRej_Kal_FHCPer2, DataCut );
174  sData_FHC_Kal_MLP_High = new Spectrum( "Data_FHC_Kal_MLP_High", RatioBins, DataLoader, kMLPCosRej_Kal_FHCHigh, DataCut );
175  sData_FHC_Kal_BDT_High = new Spectrum( "Data_FHC_Kal_BDT_High", RatioBins, DataLoader, kBDTCosRej_Kal_FHCHigh, DataCut );
176  sData_RHC_Kal_MLP_High = new Spectrum( "Data_RHC_Kal_MLP_High", RatioBins, DataLoader, kMLPCosRej_Kal_RHCHigh, DataCut );
177  sData_RHC_Kal_BDT_High = new Spectrum( "Data_RHC_Kal_BDT_High", RatioBins, DataLoader, kBDTCosRej_Kal_RHCHigh, DataCut );
178 
179  // ================
180  // = Mont Spectra =
181  // ================
182  // --- BPF spectra.
183  sMont_FHC_BPF_MLP_Per1 = new PredictionNoExtrap( *MontLoader, "Mont_FHC_BPF_MLP_Per1", RatioBins, kMLPCosRej_BPF_FHCPer1, MontCut, kNoShift, kEvtWgt );
184  sMont_FHC_BPF_BDT_Per1 = new PredictionNoExtrap( *MontLoader, "Mont_FHC_BPF_BDT_Per1", RatioBins, kBDTCosRej_BPF_FHCPer1, MontCut, kNoShift, kEvtWgt );
185  sMont_FHC_BPF_MLP_Per2 = new PredictionNoExtrap( *MontLoader, "Mont_FHC_BPF_MLP_Per2", RatioBins, kMLPCosRej_BPF_FHCPer2, MontCut, kNoShift, kEvtWgt );
186  sMont_FHC_BPF_BDT_Per2 = new PredictionNoExtrap( *MontLoader, "Mont_FHC_BPF_BDT_Per2", RatioBins, kBDTCosRej_BPF_FHCPer2, MontCut, kNoShift, kEvtWgt );
187  sMont_FHC_BPF_MLP_High = new PredictionNoExtrap( *MontLoader, "Mont_FHC_BPF_MLP_High", RatioBins, kMLPCosRej_BPF_FHCHigh, MontCut, kNoShift, kEvtWgt );
188  sMont_FHC_BPF_BDT_High = new PredictionNoExtrap( *MontLoader, "Mont_FHC_BPF_BDT_High", RatioBins, kBDTCosRej_BPF_FHCHigh, MontCut, kNoShift, kEvtWgt );
189  sMont_RHC_BPF_MLP_High = new PredictionNoExtrap( *MontLoader, "Mont_RHC_BPF_MLP_High", RatioBins, kMLPCosRej_BPF_RHCHigh, MontCut, kNoShift, kEvtWgt );
190  sMont_RHC_BPF_BDT_High = new PredictionNoExtrap( *MontLoader, "Mont_RHC_BPF_BDT_High", RatioBins, kBDTCosRej_BPF_RHCHigh, MontCut, kNoShift, kEvtWgt );
191 
192  // --- Kalman spectra.
193  sMont_FHC_Kal_MLP_Per1 = new PredictionNoExtrap( *MontLoader, "Mont_FHC_Kal_MLP_Per1", RatioBins, kMLPCosRej_Kal_FHCPer1, MontCut, kNoShift, kEvtWgt );
194  sMont_FHC_Kal_BDT_Per1 = new PredictionNoExtrap( *MontLoader, "Mont_FHC_Kal_BDT_Per1", RatioBins, kBDTCosRej_Kal_FHCPer1, MontCut, kNoShift, kEvtWgt );
195  sMont_FHC_Kal_MLP_Per2 = new PredictionNoExtrap( *MontLoader, "Mont_FHC_Kal_MLP_Per2", RatioBins, kMLPCosRej_Kal_FHCPer2, MontCut, kNoShift, kEvtWgt );
196  sMont_FHC_Kal_BDT_Per2 = new PredictionNoExtrap( *MontLoader, "Mont_FHC_Kal_BDT_Per2", RatioBins, kBDTCosRej_Kal_FHCPer2, MontCut, kNoShift, kEvtWgt );
197  sMont_FHC_Kal_MLP_High = new PredictionNoExtrap( *MontLoader, "Mont_FHC_Kal_MLP_High", RatioBins, kMLPCosRej_Kal_FHCHigh, MontCut, kNoShift, kEvtWgt );
198  sMont_FHC_Kal_BDT_High = new PredictionNoExtrap( *MontLoader, "Mont_FHC_Kal_BDT_High", RatioBins, kBDTCosRej_Kal_FHCHigh, MontCut, kNoShift, kEvtWgt );
199  sMont_RHC_Kal_MLP_High = new PredictionNoExtrap( *MontLoader, "Mont_RHC_Kal_MLP_High", RatioBins, kMLPCosRej_Kal_RHCHigh, MontCut, kNoShift, kEvtWgt );
200  sMont_RHC_Kal_BDT_High = new PredictionNoExtrap( *MontLoader, "Mont_RHC_Kal_BDT_High", RatioBins, kBDTCosRej_Kal_RHCHigh, MontCut, kNoShift, kEvtWgt );
201 
202  // --- Do it!
203  DataLoader.Go();
204  MontLoader->Go();
205 
206  // -----------------------------------------------------------
207  // --- Finally, write everything to my output file.
208  // -----------------------------------------------------------
209  TFile *outFile = new TFile( OutName.c_str(), "RECREATE" );
210  outFile->cd();
211 
212  // ================
213  // = Data Spectra =
214  // ================
215  // --- BPF spectra.
216  sData_FHC_BPF_MLP_Per1 -> SaveTo( outFile->mkdir( TString("Data_FHC_BPF_MLP_Per1") ) );
217  sData_FHC_BPF_BDT_Per1 -> SaveTo( outFile->mkdir( TString("Data_FHC_BPF_BDT_Per1") ) );
218  sData_FHC_BPF_MLP_Per2 -> SaveTo( outFile->mkdir( TString("Data_FHC_BPF_MLP_Per2") ) );
219  sData_FHC_BPF_BDT_Per2 -> SaveTo( outFile->mkdir( TString("Data_FHC_BPF_BDT_Per2") ) );
220  sData_FHC_BPF_MLP_High -> SaveTo( outFile->mkdir( TString("Data_FHC_BPF_MLP_High") ) );
221  sData_FHC_BPF_BDT_High -> SaveTo( outFile->mkdir( TString("Data_FHC_BPF_BDT_High") ) );
222  sData_RHC_BPF_MLP_High -> SaveTo( outFile->mkdir( TString("Data_RHC_BPF_MLP_High") ) );
223  sData_RHC_BPF_BDT_High -> SaveTo( outFile->mkdir( TString("Data_RHC_BPF_BDT_High") ) );
224 
225  // --- Kalman spectra.
226  sData_FHC_Kal_MLP_Per1 -> SaveTo( outFile->mkdir( TString("Data_FHC_Kal_MLP_Per1") ) );
227  sData_FHC_Kal_BDT_Per1 -> SaveTo( outFile->mkdir( TString("Data_FHC_Kal_BDT_Per1") ) );
228  sData_FHC_Kal_MLP_Per2 -> SaveTo( outFile->mkdir( TString("Data_FHC_Kal_MLP_Per2") ) );
229  sData_FHC_Kal_BDT_Per2 -> SaveTo( outFile->mkdir( TString("Data_FHC_Kal_BDT_Per2") ) );
230  sData_FHC_Kal_MLP_High -> SaveTo( outFile->mkdir( TString("Data_FHC_Kal_MLP_High") ) );
231  sData_FHC_Kal_BDT_High -> SaveTo( outFile->mkdir( TString("Data_FHC_Kal_BDT_High") ) );
232  sData_RHC_Kal_MLP_High -> SaveTo( outFile->mkdir( TString("Data_RHC_Kal_MLP_High") ) );
233  sData_RHC_Kal_BDT_High -> SaveTo( outFile->mkdir( TString("Data_RHC_Kal_BDT_High") ) );
234 
235  // ================
236  // = Mont Spectra =
237  // ================
238  // --- BPF spectra.
239  sMont_FHC_BPF_MLP_Per1 -> SaveTo( outFile->mkdir( TString("Mont_FHC_BPF_MLP_Per1") ) );
240  sMont_FHC_BPF_BDT_Per1 -> SaveTo( outFile->mkdir( TString("Mont_FHC_BPF_BDT_Per1") ) );
241  sMont_FHC_BPF_MLP_Per2 -> SaveTo( outFile->mkdir( TString("Mont_FHC_BPF_MLP_Per2") ) );
242  sMont_FHC_BPF_BDT_Per2 -> SaveTo( outFile->mkdir( TString("Mont_FHC_BPF_BDT_Per2") ) );
243  sMont_FHC_BPF_MLP_High -> SaveTo( outFile->mkdir( TString("Mont_FHC_BPF_MLP_High") ) );
244  sMont_FHC_BPF_BDT_High -> SaveTo( outFile->mkdir( TString("Mont_FHC_BPF_BDT_High") ) );
245  sMont_RHC_BPF_MLP_High -> SaveTo( outFile->mkdir( TString("Mont_RHC_BPF_MLP_High") ) );
246  sMont_RHC_BPF_BDT_High -> SaveTo( outFile->mkdir( TString("Mont_RHC_BPF_BDT_High") ) );
247 
248  // --- Kalman spectra.
249  sMont_FHC_Kal_MLP_Per1 -> SaveTo( outFile->mkdir( TString("Mont_FHC_Kal_MLP_Per1") ) );
250  sMont_FHC_Kal_BDT_Per1 -> SaveTo( outFile->mkdir( TString("Mont_FHC_Kal_BDT_Per1") ) );
251  sMont_FHC_Kal_MLP_Per2 -> SaveTo( outFile->mkdir( TString("Mont_FHC_Kal_MLP_Per2") ) );
252  sMont_FHC_Kal_BDT_Per2 -> SaveTo( outFile->mkdir( TString("Mont_FHC_Kal_BDT_Per2") ) );
253  sMont_FHC_Kal_MLP_High -> SaveTo( outFile->mkdir( TString("Mont_FHC_Kal_MLP_High") ) );
254  sMont_FHC_Kal_BDT_High -> SaveTo( outFile->mkdir( TString("Mont_FHC_Kal_BDT_High") ) );
255  sMont_RHC_Kal_MLP_High -> SaveTo( outFile->mkdir( TString("Mont_RHC_Kal_MLP_High") ) );
256  sMont_RHC_Kal_BDT_High -> SaveTo( outFile->mkdir( TString("Mont_RHC_Kal_BDT_High") ) );
257 
258  // --- And now close and we'll be done!
259  outFile->Close();
260 
261 } // --- DONE! ---
262 
263 // LocalWords: RatioBins
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
const Var kBDTCosRej_BPF_FHCHigh
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
void SetSpillCut(const SpillCut &cut)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kBDTCosRej_Kal_RHCHigh
const Cut kInCosmicTimingWindow
Is the event far from the start and ends of the spill ? For FD cosmic selection.
Definition: TimingCuts.cxx:165
const Var kBDTCosRej_BPF_FHCPer1
const Var kBDTCosRej_BPF_FHCPer2
TFile * outFile
Definition: PlotXSec.C:135
const Cut kNumuContainFD2017
Definition: NumuCuts2017.h:21
const Var kBDTCosRej_BPF_RHCHigh
virtual void Go() override
Load all the registered spectra.
std::vector< float > Spectrum
Definition: Constants.h:610
void CalcCutVals(bool isFHC, std::string sIdentifier)
Definition: CalcCutVals.C:40
static bool isFHC
const Var kBDTCosRej_Kal_FHCPer1
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
const Var kBDTCosRej_Kal_FHCHigh
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
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 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
const Var kBDTCosRej_Kal_FHCPer2
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
enum BeamMode string