makePIDCutTuning.C
Go to the documentation of this file.
1 #include "CAFAna/Core/SpectrumLoader.h" // SpectrumLoader
2 #include "CAFAna/Cuts/SpillCuts.h" // kStandardSpillCuts
3 #include "CAFAna/Prediction/PredictionNoExtrap.h" // PredictionNoExtrap
4 #include "CAFAna/Core/Spectrum.h" // Spectrum
5 #include "CAFAna/Cuts/TimingCuts.h" // kInCosmicTimingWindow,
6  // kTimingSidebandWeight
7 #include "3FlavorAna/Cuts/NueCuts2017.h" // kNue2017FDAllSamples
8 
9 #include "Utilities/func/MathUtil.h"
10 
11 #include "TString.h"
12 #include "TFile.h"
13 
14 #include <iostream>
15 #include <vector>
16 
17 #include "PIDCutTuningDefs.h"
18 
19 using namespace ana;
20 
21 //! Make predictions and cosmic spectra for tuning PID cuts or testing
22 //! efficiency/purity.
23 //!
24 //! \param[in] opt TString describing the options you would like to
25 //! run. We always search for a string contained within `opt`, so you
26 //! can combine multiple options. Some option strings are:
27 //!
28 //! - Analysis e.g. Nue, Numu, Nus
29 //! - Make Cosmics/Preds only
30 //! - PID to tune (in addition to analysis defaults) e.g. CVNe2017,
31 //! AllPIDs
32 //! - PID Cuts to apply e.g. NoPIDCut, AllPIDCuts
33 //!
34 //! Some useful opt combinations:
35 //!
36 //! - `NumuNoPIDCut` - prepares Numu sample for tuning analysis
37 //! default PIDS. We don't apply any PID cuts before tuning.
38 //! - `NueCosmicsCVNe2017NoPIDCut` - prepares Nue cosmic sample for
39 //! tuning analysis default (plus CVNe2017) PIDs.
40 //! - `NusAllPIDsAllPIDCuts` - prepares a Nus sample (both preds
41 //! and cosmic spectra) with all PIDs and all possible
42 //! permutations of PID cuts.
43 void makePIDCutTuning(const TString opt, bool useConcats=false, bool tau=true)
44 {
45  // Create loaders
46  // FD cosmic MC -- use prod4 definitions
47  std::vector<SpectrumLoader> loadersCosmic;
48  if (useConcats)
49  loadersCosmic.push_back(SpectrumLoader(
50  "prod_sumdecaf_R17-11-14-prod4reco.b_"
51  "fd_cosmic_fhcTune_full_v1_goodruns_nue2018"));
52  else
53  loadersCosmic.push_back(SpectrumLoader(
54  "prod_caf_R17-11-14-prod4reco.b_"
55  "fd_cosmic_fhcTune_full_v1_goodruns"));
56  if (useConcats)
57  loadersCosmic.push_back(SpectrumLoader(
58  "prod_sumdecaf_R17-11-14-prod4reco.a_"
59  "fd_cosmic_rhcTune_HighGain_v1_goodruns_nue2018"));
60  else
61  loadersCosmic.push_back(SpectrumLoader(
62  "prod_caf_R17-11-14-prod4reco.a_"
63  "fd_cosmic_rhcTune_HighGain_v1_goodruns_snapshot_170116"));
64 
65  // FD nonswap beam MC
66  std::vector<SpectrumLoader> loaders;
67  if (useConcats)
68  loaders.push_back(SpectrumLoader(
69  "prod_sumdecaf_R17-11-14-prod4reco.d_"
70  "fd_genie_nonswap_fhc_nova_v08_full_v1_nue2018"));
71  else
72  loaders.push_back(SpectrumLoader(
73  "prod_caf_R17-11-14-prod4reco.d_"
74  "fd_genie_nonswap_fhc_nova_v08_full_v1"));
75  if (useConcats)
76  loaders.push_back(SpectrumLoader(
77  "prod_sumdecaf_R17-11-14-prod4reco.e_"
78  "fd_genie_nonswap_rhc_nova_v08_full_v1_nue2018"));
79  else
80  loaders.push_back(SpectrumLoader(
81  "prod_caf_R17-11-14-prod4reco.e_"
82  "fd_genie_nonswap_rhc_nova_v08_full_v1"));
83 
84  // FD fluxswap beam MC
85  std::vector<SpectrumLoader> loadersSwap;
86  if (useConcats)
87  loadersSwap.push_back(SpectrumLoader(
88  "prod_sumdecaf_R17-11-14-prod4reco.d_"
89  "fd_genie_fluxswap_fhc_nova_v08_full_v1_nue2018"));
90  else
91  loadersSwap.push_back(SpectrumLoader(
92  "prod_caf_R17-11-14-prod4reco.d_"
93  "fd_genie_fluxswap_fhc_nova_v08_full_v1"));
94  if (useConcats)
95  loadersSwap.push_back(SpectrumLoader(
96  "prod_sumdecaf_R17-11-14-prod4reco.e_"
97  "fd_genie_fluxswap_rhc_nova_v08_full_v1_nue2018"));
98  else
99  loadersSwap.push_back(SpectrumLoader(
100  "prod_caf_R17-11-14-prod4reco.e_"
101  "fd_genie_fluxswap_rhc_nova_v08_full_v1"));
102 
103  // FD tauswap beam MC
104  std::vector<SpectrumLoader> loadersTau;
105  if (useConcats)
106  // Normal CAFs no concats yet
107  loadersTau.push_back(SpectrumLoader(
108  "prod_caf_R17-11-14-prod4reco.d_fd_genie_tau_fhc_nova_v08_full_v1"));
109  else
110  loadersTau.push_back(SpectrumLoader(
111  "prod_caf_R17-11-14-prod4reco.d_fd_genie_tau_fhc_nova_v08_full_v1"));
112  if (useConcats)
113  // Normal CAFs no concats yet
114  loadersTau.push_back(SpectrumLoader(
115  "prod_caf_R17-11-14-prod4reco.e_fd_genie_tau_rhc_nova_v08_full_v1"));
116  else
117  loadersTau.push_back(SpectrumLoader(
118  "prod_caf_R17-11-14-prod4reco.e_fd_genie_tau_rhc_nova_v08_full_v1"));
119 
120  // Get vectors of PID cuts
121  // noCut should ALWAYS be true here
122  std::vector<std::vector<PIDCutDef> > pidCuts = GetPIDCuts(opt);
123 
124  // Create vector to fill with predictions
125  std::vector<PredDef> preds;
126  // and another to fill with Spectra
127  std::vector<SpecDef> spectra;
128 
129  int rhcMin = 0;
130  int rhcMax = 2;
131  if (opt.Contains("RHC"))
132  rhcMin = 1;
133  if (opt.Contains("FHC"))
134  rhcMax = 1;
135 
136  // Loop over FHC and RHC
137  for(auto rhc = rhcMin; rhc < rhcMax; ++rhc)
138  {
139  loadersCosmic[rhc].SetSpillCut(kStandardSpillCuts);
140  loaders[rhc].SetSpillCut(kStandardSpillCuts);
141  loadersSwap[rhc].SetSpillCut(kStandardSpillCuts);
142  loadersTau[rhc].SetSpillCut(kStandardSpillCuts);
143 
144  TString type;
145  if (opt.Contains("NoPIDCut", TString::kIgnoreCase))
146  type = "pids";
147  else
148  type = "all";
149  std::vector<VarDef> vars = GetVars(opt, type);
150  // You may wish to pass noCut as false
151  // Don't use noCut true if using analysis bin variables
152  // --> these assume you have applied preselection cuts first!
153  std::vector<CutDef> cuts = GetCuts(opt, false, false); // pastCuts false,
154  // noCut false
155  std::vector<WeightDef> weights = GetWeights(opt, false); // noWeight false
156 
157  TString polarity = "fhc";
158  if(rhc)
159  polarity = "rhc";
160 
161  // Loop over vars
162  for(const auto var:vars)
163  {
164  // Loop over cuts
165  for(auto cut:cuts)
166  {
167  // Loop over PID cuts
168  for(auto pidCut:pidCuts[rhc])
169  {
170  // Loop over weights
171  for(auto weight:weights)
172  {
173  // Print permutation we are processing
174  std::cout << "Processing: "
175  << var.name.Data() << " "
176  << pidCut.name.Data() << " "
177  << cut.name.Data() << " "
178  << polarity << " "
179  << weight.name.Data() << std::endl;
180 
181  // Select if we want to create prediction/spectra for permutation
182  // 1. Skip any combinations that have a polarity mismatch
183  // --> only need to check vars, no polarity dependent cuts,
184  // PIDcuts or weights.
185  if (polarity == "fhc" && var.name.Contains("RHC"))
186  {
187  std::cout << "...Polarity mismatch. Skipping!..." << std::endl;
188  continue; // Don't make predictions/spectra
189  }
190  if (polarity == "rhc" && var.name.Contains("FHC"))
191  {
192  std::cout << "...Polarity mismatch. Skipping!..." << std::endl;
193  continue; // Don't make predictions/spectra
194  }
195  // Also skip kNue2017CorePresel in RHC
196  if (polarity == "rhc" && cut.name.Contains("Nue2017CorePresel"))
197  {
198  std::cout << "...Polarity mismatch. Skipping!..." << std::endl;
199  continue; // Don't make predictions/spectra
200  }
201 
202  // 2. Only apply PID cut if relevant to this var.
203  if(var.type == "pid")
204  {
205  // Only apply PID cut if relevant for PID
206  // --> Other Vars --> apply all PID Cuts
207  if(!(pidCut.name.Contains("NoPIDCut") ||
208  pidCut.name.Contains(var.shortName.Data()) ||
209  opt.Contains("AllPIDCuts", TString::kIgnoreCase)))
210  {
211  std::cout << "...Irrelevant PID cut. Skipping!..." << std::endl;
212  continue; // Don't make predictions/spectra
213  }
214  // Also skip making predictions for zoomed PID variables unless
215  // we're making a PID cut
216  if (pidCut.name.Contains("NoPIDCut") &&
217  var.name.Contains("Zoom", TString::kIgnoreCase))
218  {
219  std::cout << "...No Zoom. Skipping!..." << std::endl;
220  continue; // Don't make predictions/spectra
221  }
222  }
223  else if(var.type == "ana18") // 2018 analysis binning tests
224  {
225  // Don't apply any PID cut other than short simple
226  if (!(pidCut.name.Contains("CVNeShortSimple") ||
227  pidCut.name.Contains("NoPIDCut")))
228  {
229  std::cout << "...Irrelevant PID cut. Skipping!..." << std::endl;
230  continue; // Don't make predictions/spectra
231  }
232  }
233 
234  if (!opt.Contains("Cosmics", TString::kIgnoreCase))
235  {
236  // Add prediction
237  if (tau)
238  preds.push_back({
239  new PredictionNoExtrap(
240  loaders[rhc], loadersSwap[rhc], loadersTau[rhc],
241  var.name.Data(), var.binning, var.var,
242  cut.cut && pidCut.cut, kNoShift, weight.weight),
243  cut.name.Data(), pidCut.name.Data(),
244  var.name.Data(), polarity, weight.name.Data()});
245  else
246  preds.push_back({
247  new PredictionNoExtrap(
248  loaders[rhc], loadersSwap[rhc],
249  var.name.Data(), var.binning, var.var,
250  cut.cut && pidCut.cut, kNoShift, weight.weight),
251  cut.name.Data(), pidCut.name.Data(),
252  var.name.Data(), polarity, weight.name.Data()});
253  std::cout << "...Making prediction..." << std::endl;
254  }// if not cosmics
255 
256  if (!opt.Contains("Preds", TString::kIgnoreCase))
257  {
258  // Add cosmic Spectra
259  spectra.push_back({
260  new Spectrum(var.name.Data(), var.binning, loadersCosmic[rhc],
261  var.var, (cut.cut && kInTimingSideband &&
262  pidCut.cut), kNoShift, weight.weight),
263  cut.name.Data(), pidCut.name.Data(),
264  var.name.Data(), polarity, weight.name.Data()});
265  std::cout << "...Making cosmic spectrum..." << std::endl;
266  }// if not Preds
267  }// loop over weights
268  }// loop over PID cuts
269  }// loop over cuts
270  }// loop over vars
271  }// loop over fhc/rhc
272 
273  std::cout << "Loaders Go!" << std::endl;
274  // Loop over FHC and RHC
275  for(auto rhc = rhcMin; rhc < rhcMax; ++rhc)
276  {
277  TString polarity = "fhc";
278  if(rhc)
279  polarity = "rhc";
280 
281  if (!opt.Contains("Preds", TString::kIgnoreCase))
282  {
283  loadersCosmic[rhc].Go();
284  std::cout << " --> loadersCosmic..." << polarity.Data() << std::endl;
285  }
286  if (!opt.Contains("Cosmics", TString::kIgnoreCase))
287  {
288  loaders[rhc].Go();
289  std::cout << " --> loaders..." << polarity.Data() << std::endl;
290  loadersSwap[rhc].Go();
291  std::cout << " --> loadersSwap..." << polarity.Data() << std::endl;
292  if (tau)
293  {
294  loadersTau[rhc].Go();
295  std::cout << " --> loadersTau..." << polarity.Data() << std::endl;
296  }
297  }
298  }
299 
300  // Write predictions and spectra to file
301  TString filename = opt + ".root";
302  auto outFile = new TFile(filename, "recreate");
303 
304  if (!opt.Contains("Preds", TString::kIgnoreCase))
305  {
306  std::cout << "writing cosmic spectra..." << std::endl;
307  for(auto spec:spectra)
308  {
309  std::cout << "cosmic_spec_"
310  << spec.cutName << "_"
311  << spec.pidCutName << "_"
312  << spec.varName << "_"
313  << spec.polarity << "_"
314  << spec.weightName << std::endl;
315  spec.spectrum->SaveTo(outFile,
316  "cosmic_spec_" + spec.cutName + "_" +
317  spec.pidCutName + "_" + spec.varName + "_" +
318  spec.polarity + "_" + spec.weightName);
319  }
320  }
321  if (!opt.Contains("Cosmics", TString::kIgnoreCase))
322  {
323  std::cout << "writing predictions..." << std::endl;
324  for(auto pred:preds)
325  {
326  std::cout << "pred_nxp_"
327  << pred.cutName << "_"
328  << pred.pidCutName << "_"
329  << pred.varName << "_"
330  << pred.polarity << "_"
331  << pred.weightName << std::endl;
332  pred.prediction->SaveTo(outFile,
333  "pred_nxp_" + pred.cutName + "_" +
334  pred.pidCutName + "_" + pred.varName + "_" +
335  pred.polarity + "_" + pred.weightName);
336  }
337  }
338 
339  std::cout << "Ending script..." << std::endl;
340 }// end of makePredictionsPIDTuning
std::vector< WeightDef > GetWeights(TString opt, bool noWeight=false)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void makePIDCutTuning(const TString opt, bool useConcats=false, bool tau=true)
const Var weight
string filename
Definition: shutoffs.py:106
const Cut kInTimingSideband([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return(kInTimingSideband_before(sr)|| kInTimingSideband_after(sr));else return(kInTimingSideband_before(sr)|| kInTimingSideband_afterA(sr)|| kInTimingSideband_afterB(sr));}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_after.Livetime(spill));else return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_afterA.Livetime(spill)+ kInTimingSideband_afterB.Livetime(spill));}, [](const caf::SRSpillProxy *spill){return 0;})
Definition: TimingCuts.h:12
TFile * outFile
Definition: PlotXSec.C:135
std::vector< CutDef > GetCuts(TString opt, bool pastCuts=false, bool noCut=true)
Var weights
const std::map< std::pair< std::string, std::string >, Variable > vars
std::vector< std::vector< PIDCutDef > > GetPIDCuts(TString opt, bool noCut=true)
std::vector< float > Spectrum
Definition: Constants.h:610
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const Cut cut
Definition: exporter_fd.C:30
const std::vector< VarDef > GetVars(TString opt, TString type="all")
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
std::vector< Loaders * > loaders
Definition: syst_header.h:386
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.