make_nue_xsec_pred.C
Go to the documentation of this file.
2 //#include "3FlavorAna/Ana2017/nue/SystVariations.h"
3 #include "CAFAna/Core/ISyst.h"
4 #include "CAFAna/Systs/XSecSystLists.h" // getAllXsecSysts_2017
5 #include "CAFAna/Vars/GenieWeights.h" // kXSecCVWgt2017
6 #include "CAFAna/Vars/PPFXWeights.h" // kPPFXFluxCVWgt
7 #include "CAFAna/Analysis/Prod3Loaders.h" // kNueConcat, kNueDecaf
8 #include "CAFAna/Vars/HistAxes.h" // kNumuNonQEAxisFirstAna
9 //#include "CAFAna/Cuts/NueCutsSecondAna.h"
10 #include "3FlavorAna/Vars/NueVarsExtra.h" // kShwCalE, kHadCalE
11 #include "3FlavorAna/Vars/NueVars.h" // kNueEnergy2017
12 #include "CAFAna/Vars/Vars.h" // kCVNe
13 #include "3FlavorAna/Cuts/NueCuts2017.h" // kNue2017Binning, kNue2017AnaBin
19 //#include "CAFAna/Core/Registry.h"
20 //#include "CAFAna/Systs/XSecSysts.h"
21 //#include "CAFAna/Analysis/Calcs.h"
22 //#include "CAFAna/Cuts/Cuts.h"
23 //#include "CAFAna/Cuts/NumuCuts.h"
24 //#include "CAFAna/Cuts/SpillCuts.h"
25 
27 
28 #include "TString.h"
29 #include "TFile.h"
30 
31 #include <vector>
32 #include <iostream>
33 #include <string>
34 
35 using namespace ana;
36 
37 // Structure to hold weights --> only way to combine weights easily
38 struct Weight
39 {
40  const TString name;
42 };
43 
44 
45 void make_nue_xsec_pred(TString suffix="test", bool useConcat=true,
46  int startIdx=0, int stopIdx=-1)
47 {
48  // Filename for output prediction file
49  TString filename = "pred_xsec_fhc_" + suffix + ".root";
50 
51  // Make PredictionWriter instance
52  PredictionWriter* writer = new PredictionWriter(filename);
53 
54  // Set analysis name
55  TString analysisName = "pred_propDecomp";
56 
57  // Get vector of xsec systematics
58  std::vector<const ISyst*> systs = getAllXsecSysts_2017();
59  std::cout << "[make_nue_xsec_pred] Found " << systs.size()
60  << " xsec systematics!" << std::endl;
61  if (stopIdx < 0) // assume all systematics, set as length of vector
62  stopIdx = systs.size();
63  std::cout << "[make_nue_xsec_pred] Processing systematics:" << std::endl;
64  for (std::vector<const ISyst*>::iterator syst = systs.begin()+startIdx;
65  (syst != systs.begin()+stopIdx) && (syst != systs.end()); ++syst)
66  std::cout << std::distance(systs.begin(), syst) << ": "
67  << (*syst)->ShortName() << std::endl;
68 
69  // Set sigmas - only want 0 (no shift) once, so handle separately
70  std::vector<int> sigmas = {-2, -1, 1, 2};
71 
72  // Create vector/map of weights
73  std::vector<Weight> weights;
74  weights.push_back({"ppfx_xsec", kPPFXFluxCVWgt * kXSecCVWgt2017});
75  weights.push_back({"ppfx_only", kPPFXFluxCVWgt});
76 
77  // Setup Loaders
78  ECAFType type;
79  if (useConcat)
80  type = kNueConcat;
81  else
82  type = kDecaf;
84 
85  // Define axes and energy variables for predictions
86  const HistAxis axisNumu = kNumuNonQEAxisFirstAna;
87 
88  struct HistDef
89  {
90  const TString name;
91  HistAxis axis;
92  };
93  std::vector<HistDef> hists;
94 // hists.push_back({"shwCalE",
95 // {"Shower calE (GeV)",
96 // Binning::Simple(30.0, 0.0, 3.0), kShwCalE}});
97 // hists.push_back({"hadCalE",
98 // {"calE_{had} (GeV)",
99 // Binning::Simple(30.0, 0.0, 3.0), kHadCalE}});
100  //hists.push_back({"recoNuE",
101 // {"reco E_{#nu} (GeV)",
102 // Binning::Simple(50.0, 0.0, 5.0), kNueEnergy2017}});
103  //hists.push_back({"id_cvne",
104  // {"CVN e", Binning::Simple(40.0, 0.70, 1.1), kCVNe}});
105  hists.push_back({"EnergyCVN2D",
106  {"Nue Energy / Selection Bin",
107  kNue2017Binning, kNue2017AnaBin}}); // peripheral binning
108 
109  // Create vector/map to store predictions
110  std::vector<std::pair<IPrediction*, TString> > predictions;
111 
112  //TODO: add vector of pointers to Pred Generators
113  //* No extrap
114  //* proportional
115  //* combo
116  // try pred no extrap with numu cut for proof of principle
117 
118  // First do Nominal
119  // Loop over axes and energy Vars
120  for (auto hist:hists)
121  {
122  for (auto weight:weights)
123  {
124  TString path; // In output prediction root file
125  // Extrap first
126  path = writer->ReadPrediction(
127  analysisName, "Nominal", 0, true, hist.name, weight.name);
128 
129  auto extrapGen = NuePropExtrapGenerator(
130  hist.axis, axisNumu, kNue2017FDAllSamples, kNue2017NDCVNSsb,
131  kNumuND, kNoShift, weight.weight);
132  auto predProp = extrapGen.Generate(loaders); // no shift
133  predictions.push_back({predProp.release(), path});
134 
135  // No Extrap
136  path = writer->ReadPrediction(
137  analysisName, "Nominal", 0, false, hist.name, weight.name);
138 
139  auto noExtrapGen = NoExtrapGenerator(
140  hist.axis, kNue2017FDAllSamples, weight.weight);
141  auto predNoExtrap = noExtrapGen.Generate(loaders); // no shift
142  predictions.push_back({predNoExtrap.release(), path});
143 
144  }// Loop over weights
145  }// Loop over axes and energy Vars
146 
147  // Loop over systematics
148  for (std::vector<const ISyst*>::iterator syst = systs.begin()+startIdx;
149  (syst != systs.begin()+stopIdx) && (syst != systs.end()); ++syst)
150  {
151  // Loop over sigmas
152  for (auto sigma:sigmas)
153  {
154  SystShifts shift = SystShifts((*syst), static_cast<double>(sigma));
155  // Loop over axes and energy Vars
156  for (auto hist:hists)
157  {
158  // Loop over weights
159  for (auto weight:weights)
160  {
161  TString path; // In output prediction root file
162  // Extrap first
163  path = writer->ReadPrediction(analysisName, (*syst)->ShortName(),
164  sigma, true, hist.name, weight.name);
165 
166  auto extrapGen = NuePropExtrapGenerator(
167  hist.axis, axisNumu, kNue2017FDAllSamples, kNue2017NDCVNSsb,
168  kNumuND, kNoShift, weight.weight);
169  auto predProp = extrapGen.Generate(loaders, shift);
170  predictions.push_back({predProp.release(), path});
171 
172  // No extrap
173  path = writer->ReadPrediction(analysisName, (*syst)->ShortName(),
174  sigma, false, hist.name, weight.name);
175 
176  auto noExtrapGen = NoExtrapGenerator(
177  hist.axis, kNue2017FDAllSamples, weight.weight);
178  auto predNoExtrap = noExtrapGen.Generate(loaders, shift);
179  predictions.push_back({predNoExtrap.release(), path});
180 
181  }// Loop over weights
182  }// Loop over axes and energy Vars
183  }// Loop over sigmas
184  }// Loop over systNames
185 
186  loaders.Go();
187 
188  // Loop over again to write predictions
189  std::vector<std::pair<IPrediction*, TString> >::iterator iPrediction =
190  predictions.begin();
191 
192  // Do Nominal first
193  // Loop over axes and energy Vars
194  for (auto hist:hists)
195  {
196  // Loop over weights
197  for (auto weight:weights)
198  {
199  // Write extrap prediction
200  TString path;
201  path = writer->ReadPrediction(
202  analysisName, "Nominal", 0, true, hist.name, weight.name);
203 
204  // Check indexing
205  if (!path.EqualTo(iPrediction->second)) // incorrect save location
206  std::cout << "[make_nue_xsec_pred] WARNING: prediction path mismatch!"
207  << std::endl;
208 
209  writer->WritePrediction(iPrediction->first, analysisName, "Nominal",
210  0, true, hist.name, weight.name);
211  ++iPrediction; // increment
212 
213  // Write no extrap prediction
214  path = writer->ReadPrediction(
215  analysisName, "Nominal", 0, false, hist.name, weight.name);
216 
217  // Check indexing
218  if (!path.EqualTo(iPrediction->second)) // incorrect save location
219  std::cout << "[make_nue_xsec_pred] WARNING: prediction path mismatch!"
220  << std::endl;
221 
222  writer->WritePrediction(iPrediction->first, analysisName, "Nominal",
223  0, false, hist.name, weight.name);
224  ++iPrediction; // increment
225  }// Loop over weights
226  }// Loop over axes and energy Vars
227 
228  // Loop over systematics
229  for (std::vector<const ISyst*>::iterator syst = systs.begin()+startIdx;
230  (syst != systs.begin()+stopIdx) && (syst != systs.end()); ++syst)
231  {
232  // Loop over sigmas
233  for (auto sigma:sigmas)
234  {
235  // Loop over axes and energy Vars
236  for (auto hist:hists)
237  {
238  // Loop over weights
239  for (auto weight:weights)
240  {
241  // Write extrap prediction
242  TString path;
243  path = writer->ReadPrediction(analysisName, (*syst)->ShortName(),
244  sigma, true, hist.name, weight.name);
245 
246  // Check indexing
247  if (!path.EqualTo(iPrediction->second)) // incorrect save locatio
248  std::cout << "[make_nue_xsec_pred] WARNING: "
249  << "prediction path mismatch!" << std::endl;
250 
251  writer->WritePrediction(iPrediction->first, analysisName,
252  (*syst)->ShortName(), sigma,
253  true, hist.name, weight.name);
254  ++iPrediction; // increment
255 
256  // Write no extrap prediction
257  path = writer->ReadPrediction(analysisName, (*syst)->ShortName(),
258  sigma, false, hist.name, weight.name);
259 
260  // Check indexing
261  if (!path.EqualTo(iPrediction->second)) // incorrect save locatio
262  std::cout << "[make_nue_xsec_pred] WARNING: "
263  << "prediction path mismatch!" << std::endl;
264 
265  writer->WritePrediction(iPrediction->first, analysisName,
266  (*syst)->ShortName(), sigma,
267  false, hist.name, weight.name);
268  ++iPrediction;
269  }// Loop over weights
270  }// Loop over axes and energy Vars
271  }// Loop over sigmas
272  }// Loop over systNames
273 }// make_nue_xsec_pred
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod3Loaders.h:40
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const TString name
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut kNue2017NDCVNSsb
Definition: NueCuts2017.h:302
string filename
Definition: shutoffs.py:106
const Binning kNue2017Binning
TString hists[nhists]
Definition: bdt_com.C:3
unsigned distance(const T &t1, const T &t2)
Generates FD-only predictions (no extrapolation)
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
Generates extrapolated Nue predictions using ProportionalDecomp.
void WritePrediction(const IPrediction *prediction, const TString &analysisName, const TString &systName, const int &sigma, const bool &extrap=true, const TString &energyVar="EnergyCVN2D", const TString &weight="kXSecCVWgt2017")
void make_nue_xsec_pred(TString suffix="test", bool useConcat=true, int startIdx=0, int stopIdx=-1)
const Cut kNumuND
Definition: NumuCuts.h:55
Var weights
double sigma(TH1F *hist, double percentile)
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
const Var kNue2017AnaBin([](const caf::SRProxy *sr){int selBin=kNue2017SelectionBin(sr);float nuE=kNueEnergy2017(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Definition: NueCuts2017.h:311
std::vector< Loaders * > loaders
Definition: syst_header.h:386
const TString ReadPrediction(const TString &analysisName, const TString &systName, const int &sigma, const bool &extrap=true, const TString &energyVar="EnergyCVN2D", const TString &weight="kXSecCVWgt2017")
const HistAxis axisNumu
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
Template for Var and SpillVar.
std::vector< const ISyst * > getAllXsecSysts_2017()
Get master XSec syst list for 2017 analyses.
const Cut kNue2017FDAllSamples
Our FD selection including all samples, for making predictions, etc.
Definition: NueCuts2017.h:155
ECAFType
Definition: Loaders.h:19