makePrediction.C
Go to the documentation of this file.
1 // This tutorial covers:
2 // - Loaders, including Prod4Loaders
3 // - Generating predictions
4 // - Saving predictions to file
5 
6 #include "CAFAna/Core/Loaders.h" // ECAFType, kRHC, kFHC
7 #include "CAFAna/Analysis/Prod4Loaders.h" // Prod4NomLoaders
8 #include "CAFAna/Cuts/SpillCuts.h" // kStandardSpillCuts
9 #include "CAFAna/Vars/HistAxes.h" // HistAxis
10 #include "CAFAna/Core/Binning.h" // kNumuCCEOptimisedBinning
11 #include "CAFAna/Vars/Vars.h" // kMuE
12 #include "CAFAna/Cuts/NumuCuts2018.h" // kNumuCutFD2018
13 #include "CAFAna/Vars/PPFXWeights.h" // kPPFXFluxCVWgt
14 #include "CAFAna/Vars/GenieWeights.h" // kXSecCVWgt2018
15 #include "CAFAna/Core/SystShifts.h" // kNoShift
16 #include "CAFAna/Prediction/IPrediction.h" // IPrediction
17 #include "CAFAna/Prediction/PredictionGenerator.h" // NoExtrapGenerator
18 
19 #include "TString.h"
20 #include "TFile.h"
21 
22 #include <vector>
23 
24 using namespace ana;
25 
26 // Try running with default options:
27 //
28 // cafe -bq -s 500 -l 5 -ss makePrediction.C
29 //
30 // To run over a single RHC concat file try:
31 //
32 // cafe -bq -l 1 -ss makePrediction.C RHCNumuEnergy_numuconcat
33 //
34 
35 void makePrediction(TString option="FHCNumuEnergy_fullCAFs")
36 {
37  // Step 1: Set up loaders
38 
39  // Get flux (FHC or RHC) from option
40  if (!(option.Contains("FHC") || option.Contains("RHC")))
41  {
42  std:: cerr << "Please mention FHC or RHC in option" << std::endl;
43  abort(); // unable to determine if FHC or RHC
44  }
45  auto flux = option.Contains("RHC") ? Loaders::kRHC : Loaders::kFHC;
46 
47  // Define CAF type --> default is full CAFs
48  auto caftype = ana::ECAFType::kFullCAF;
49  if (option.Contains("numuconcat")) caftype = ana::ECAFType::kNumuConcat;
50 
51  // Always want full dataset
52  const TString period="full";
53 
54  // Create Prod4Loaders object
55  auto loaders = new Prod4NomLoaders(caftype, flux,
56  period.Data(), period.Data());
57  loaders->SetSpillCut(kStandardSpillCuts); // apply standard spill cuts
58 
59  // Step 2: Define a spectrum options
60 
61  // Axis --> numu optimised binning, energy axis
62  const HistAxis axis("Reconstructed Muon Energy (GeV)",
64  // Cut --> numu 2018 FarDet cut
65  auto fdCut = kNumuCutFD2018;
66  // Weight --> standard for 2018 analysis
67  const auto weight = kPPFXFluxCVWgt * kXSecCVWgt2018;
68 
69  // Define shifts
70  // For now we will focus on non-systematically shifted predictions
71  auto shift = kNoShift;
72 
73  // Step 3: Define containers
74  // Always best to define some vectors of structs to put things in
75  // We will do this here too, even though we are only creating on prediction
76 
77  // To hold prediction generators...
78  struct GenDef
79  {
80  const IPredictionGenerator * gen; // generator
81  const TString xpType; // description of extrapolation type
82  const TString description; // analysis description
83  };
84 
85  // To hold predictions
86  struct PredDef
87  {
88  const IPrediction * pred; // prediction
89  const TString xpType; // description of extrapolation type
90  const TString description; // analysis description
91  // For systematics predictions, could also add:
92  // const TString systNname; // systematic name (default 'NoShift')
93  // const TString sigma_name; // sigma name (-2, -1, 0, 1, 2)
94  };
95 
96  std::vector<GenDef> gens;
97  std::vector<PredDef> preds;
98 
99  // Step 4: Create generators
100  // For multiple Vars, Cuts, axes etc. put this in nested loops
101  gens.push_back(
102  {new NoExtrapGenerator(axis, fdCut, weight), "pred_nxp", "NumuEnergy"});
103 
104  // Step 5: Loop over generators and create predictions
105  for (auto gen:gens)
106  {
107  IPrediction * pred = gen.gen->Generate(*loaders, shift).release();
108  preds.push_back({pred, "pred_nxp", "NumuEnergy"});
109  }
110 
111  // Step 6: Once all predictions/spectra have been filled...
112  loaders->Go();
113 
114  // Step 7: Write predictions to file
115  TString filename = option + ".root";
116  auto outFile = new TFile(filename, "recreate");
117 
118  // Loop over vector of predictions
119  std::cout << "writing predictions..." << std::endl;
120  for (auto pred:preds)
121  {
122  const TString name = pred.xpType + "_" + pred.description;
123  std::cout << name << std::endl;
124  pred.pred->SaveTo(outFile->mkdir(name));
125  }
126 
127  std::cout << "Done!" << std::endl;
128 }
129 
130 // Tutorial extensions
131 // 1. How would you extend this macro to produce two predictions, one for Numu
132 // energy as it does currently and one for Nue energy? Hint: create a vector
133 // of Vars and HistAxis obects to loop over.
134 // 2. Think about how you would extend this to produce a prediction from one of
135 // the systematically shifted samples e.g. Absolute Calibration Up.
const XML_Char * name
Definition: expat.h:151
std::map< TString, IPredictionGenerator * > gens
Definition: syst_header.h:387
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
Struct to hold prediction information.
OStream cerr
Definition: OStream.cxx:7
string filename
Definition: shutoffs.py:106
Generates FD-only predictions (no extrapolation)
Loaders::FluxType flux
TFile * outFile
Definition: PlotXSec.C:135
void makePrediction(TString option="FHCNumuEnergy_fullCAFs")
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
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
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Given loaders and an MC shift, Generate() generates an IPrediction.
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
const Var kMuE
Definition: NumuVars.h:22