PredictionSystNumu2017.cxx
Go to the documentation of this file.
5 
8 
9 #include "OscLib/IOscCalc.h"
10 
11 #include "TObjString.h"
12 
13 #include <iostream>
14 
15 namespace ana
16 {
17 
18  REGISTER_LOADFROM("PredictionSystNumu2017", IPrediction, PredictionSystNumu2017);
19 
20  //----------------------------------------------------------------------
21 
22  const DummyNumu2017Syst kNumu2017MECq0ShapeSyst ("numu2017MECq0Shape" , "MECq0Shape" , {-2, -1, 0, +1, +2}, false);
23  const DummyNumu2017Syst kNumu2017MaCCQE_reducedSyst ("numu2017MaCCQE_reduced" , "MaCCQE_reduced" , {-2, -1, 0, +1, +2}, false);
24  const DummyNumu2017Syst kNumu2017RPAShapeRESSyst ("numu2017RPAShapeRES" , "RPAShapeRES" , {0, +1, +2}, false);
25  const DummyNumu2017Syst kNumu2017RPAShapeEnhSyst ("numu2017RPAShapeenh" , "RPAShapeEnh" , {-2, -1, 0, +1, +2}, false);
26  const DummyNumu2017Syst kNumu2017MaCCRESSyst ("numu2017MaCCRES" , "MaCCRES" , {-2, -1, 0, +1, +2}, false);
27  const DummyNumu2017Syst kNumu2017MaNCRESSyst ("numu2017MaNCRES" , "MaNCRES" , {-2, -1, 0, +1, +2}, false);
28  const DummyNumu2017Syst kNumu2017MvCCRESSyst ("numu2017MvCCRES" , "MvCCRES" , {-2, -1, 0, +1, +2}, false);
29  const DummyNumu2017Syst kNumu2017DISvnCC1piSyst ("numu2017DISvnCC1pi" , "DISvnCC1pi" , {-2, -1, 0, +1, +2}, false);
30  const DummyNumu2017Syst kNumu2017AllSmallSysts ("numu2017SumSmallXSecJoint2017", "SumSmallXSecJoint2017", {-2, -1, 0, +1, +2}, false);
31  const DummyNumu2017Syst kNumu2017NumuSmallSysts ("numu2017SumSmallXSecNumu2017" , "SumSmallXSecNumu2017" , {-2, -1, 0, +1, +2}, false);
32 
33  const DummyNumu2017Syst kNumu2017AllBeamSysts("numu2017beamTransportComb", "beamTransportComb", {-2, -1, 0, +1, +2}, false);
34 
35  const DummyNumu2017Syst kNumu2017RelHadEScale("numu2017RelHadEScale2017", "RelHadEScale", {-2, -1, 0, +1, +2}, false);
36  const DummyNumu2017Syst kNumu2017AbsHadEScale("numu2017AbsHadEScale2017", "AbsHadEScale", {-2, -1, 0, +1, +2}, false);
37  const DummyNumu2017Syst kNumu2017RelMuEScale ("numu2017RelMuEScale2017" , "RelMuEScale" , {-2, -1, 0, +1, +2}, false);
38  const DummyNumu2017Syst kNumu2017AbsMuEScale ("numu2017AbsMuEScale2017" , "AbsMuEScale" , {-2, -1, 0, +1, +2}, false);
39 
40  const DummyNumu2017Syst kNumu2017AbsCalibSyst ("numu2017Calibration" , "AbsCalib" , {-1, 0, +1}, false);
41  const DummyNumu2017Syst kNumu2017RelCalibSyst ("numu2017RelativeCalib", "RelativeCalib", {-2, 0, +2}, false);
42  const DummyNumu2017Syst kNumu2017CalibShapeSyst("numu2017CalibShape" , "CalibShape" , {0, +1}, false);
43  const DummyNumu2017Syst kNumu2017LightLevelSyst("numu2017Lightlevel" , "LightLevel" , {-1, 0, +1}, false);
44  const DummyNumu2017Syst kNumu2017CherenkovSyst ("numu2017Cherenkov" , "Cherenkov" , {0, +1}, false);
45 
46  const DummyNumu2017Syst kNumu2017PPFXPC00Syst ("numu2017ppfx_pc00", "PPFXPC00Sys", {-2, -1, 0, +1, +2}, false);
47  const DummyNumu2017Syst kNumu2017PPFXPC01Syst ("numu2017ppfx_pc01", "PPFXPC01" , {-2, -1, 0, +1, +2}, false);
48  const DummyNumu2017Syst kNumu2017PPFXPC02Syst ("numu2017ppfx_pc02", "PPFXPC02" , {-2, -1, 0, +1, +2}, false);
49  const DummyNumu2017Syst kNumu2017PPFXPC03Syst ("numu2017ppfx_pc03", "PPFXPC03" , {-2, -1, 0, +1, +2}, false);
50  const DummyNumu2017Syst kNumu2017PPFXPC04Syst ("numu2017ppfx_pc04", "PPFXPC04" , {-2, -1, 0, +1, +2}, false);
51 
52  //const DummyNumu2017Syst kNumu2017AbsCalibSyst2Sig("numu2017Calibration", "AbsCalib2Sig", {-2,-1, 0, +1, +2}, false);
53 
54  const DummyNumu2017Syst kNumu2017NormSystSig ("numu2017normSig" , "Sig. Norm." , {0, +1}, false);
55  const DummyNumu2017Syst kNumu2017NormSystBkg ("numu2017normBkg" , "Bkg. Norm." , {0, +1}, false);
56  const DummyNumu2017Syst kNumu2017NormSystBoth("numu2017NormBoth", "Sig+Bkg Norm.", {0, +1}, false);
57 
58  //----------------------------------------------------------------------
61  int WhichQuant,
62  const std::vector<const DummyNumu2017Syst*>& systs,
63  const std::string fname)
65  {
66  fOscOrigin = osc->Copy();
67 
68  std::string extrapStr = "pred_";
69  if (extrap == ENumu2017ExtrapType::kNuMu) extrapStr += "xp_numu";
70  if (extrap == ENumu2017ExtrapType::kNumuNoExtrap) extrapStr += "nxp";
71  extrapStr += "_";
72 
73  std::string VarName = "numu_pred_Energy_Quant"+std::to_string(WhichQuant);
74 
75  std::cout << "\nLooking at VarName " << VarName << " --> Loading nominals now..." << std::endl;
76 
77  IPrediction *nom;
78  std::string NomVar = extrapStr+"Nominal/noShift/"+VarName;
79  fPredNom = ana::LoadFromFile<IPrediction>(fname, NomVar);
80  nom = ana::LoadFromFile<IPrediction>(fname, NomVar).release();
81  // --- The light level syst uses an alternative nominal. Grab that now.
82  std::string LNomVar = extrapStr+"Lightlevel/noShift/"+VarName;
83  IPrediction* nomLightLevel = ana::LoadFromFile<IPrediction>(fname, LNomVar).release();
84 
85  for(const DummyNumu2017Syst* syst: systs){
86  if(syst == &kNumu2017NormSystBoth ||
87  syst == &kNumu2017NormSystSig ||
88  syst == &kNumu2017NormSystBkg) continue;
89 
90  // Hack. For sake of distinct names, every syst's fist 8 characters are numu2017
91  std::string systName = syst->ShortName();
92  systName.erase(0,8);
93 
94  std::cout << "Loading " << systName << "..." << std::endl;
95 
96  ShiftedPreds sp;
97  sp.systName = systName;
98 
99  for(int sigma: syst->Shifts()){
100  sp.shifts.push_back(sigma);
101 
102  if (sigma == 0){
103  if (systName == "Lightlevel" || systName == "Cherenkov") {
104  std::cout << "\tLooking at systName " << systName << ", so pushing back nomLightLevel " << std::endl;
105  sp.preds.push_back(nomLightLevel);
106  } else {
107  sp.preds.push_back(nom);
108  }
109  continue;
110  }
111 
112  std::string sigmaStr = "/";
113  if(sigma == -2) sigmaStr += "minusTwo/";
114  if(sigma == -1) sigmaStr += "minusOne/";
115  if(sigma == +1) sigmaStr += "plusOne/";
116  if(sigma == +2) sigmaStr += "plusTwo/";
117 
118 
119  auto shCore = ana::LoadFromFile<PredictionExtrap> (fname, extrapStr + systName + sigmaStr +VarName).release();
120  sp.preds.push_back(shCore);
121  } // end for sigma
122 
123  fPreds.emplace(syst, sp);
124  } // end for syst
125 
126  // Calculate all the ratios etc
127  InitFits();
128 
129  AddNormSysts(systs);
130  }
131 
132  //----------------------------------------------------------------------
134  {
135  }
136 
137  //----------------------------------------------------------------------
138  void PredictionSystNumu2017::AddNormSysts(const std::vector<const DummyNumu2017Syst*>& systs)
139  {
140  TH1* h = Predict(fOscOrigin).ToTH1(1);
141  const int nBins = h->GetNbinsX();
142  delete h;
143 
144 
145  // Got this from Alex. TODO add reference
146  // pythag == sum in quadrature
147  // leaving Sig and bkg hooks in place, but they are not used anywhere
148  const double normBothSyst = 3.5/100;
149  const double normBkgSyst = 3.5/100;
150  const double normSigSyst = 3.5/100;
151 
152  // TODO find a way to write these three norm systs as a loop
153  // if(std::count(systs.begin(), systs.end(), &kNumu2017NormSystBoth)){
154  if(std::count_if(systs.begin(), systs.end(),
155  [](const ISyst* s){return s->ShortName() == "numu2017NormBoth";})){
156  ShiftedPreds normsp;
157  normsp.systName = "NormBoth";
158  normsp.shifts = {0};
159  normsp.preds = {fPredNom.get()};
160  for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
161  normsp.fits[coeff].resize(nBins+2);
162  for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
163  normsp.fits[coeff][binIdx].push_back(Coeffs(0, 0, normBothSyst, 1));
164  }
165  }
166  normsp.FillRemaps();
167  fPreds.emplace(&kNumu2017NormSystBoth, normsp);
168  }
169 
170 
171  // if(std::count(systs.begin(), systs.end(), &kNumu2017NormSystBkg)){
172  if(std::count_if(systs.begin(), systs.end(),
173  [](const ISyst* s){return s->ShortName() == "numu2017normBkg";})){
174  ShiftedPreds normbksp;
175  normbksp.systName = "NumuNormBkg";
176  normbksp.shifts = {0};
177  normbksp.preds = {fPredNom.get()};
178  for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
179  normbksp.fits[coeff].resize(nBins+2);
180  for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
181  normbksp.fits[coeff][binIdx].push_back(Coeffs(0, 0, (coeff != kNumuSurv) ? normBkgSyst : 0, 1));
182  }
183  }
184  normbksp.FillRemaps();
185  fPreds.emplace(&kNumu2017NormSystBkg, normbksp);
186  }
187 
188 
189  // if(std::count(systs.begin(), systs.end(), &kNumu2017NormSystSig)){
190  if(std::count_if(systs.begin(), systs.end(),
191  [](const ISyst* s){return s->ShortName() == "numu2017normSig";})){
192  ShiftedPreds normsigsp;
193  normsigsp.systName = "NumuNormSig";
194  normsigsp.shifts ={0};
195  normsigsp.preds = {fPredNom.get()};
196  for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
197  normsigsp.fits[coeff].resize(nBins+2);
198  for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
199  normsigsp.fits[coeff][binIdx].push_back(Coeffs(0, 0, (coeff == kNumuSurv) ? normSigSyst : 0, 1));
200  }
201  }
202  normsigsp.FillRemaps();
203  fPreds.emplace(&kNumu2017NormSystSig, normsigsp);
204  }
205  }
206 
207  //----------------------------------------------------------------------
208  std::unique_ptr<PredictionSystNumu2017> PredictionSystNumu2017::LoadFrom(TDirectory* dir, const std::string& name)
209  {
210  dir = dir->GetDirectory(name.c_str()); // switch to subdir
211  assert(dir);
212 
213  TObjString* tag = (TObjString*)dir->Get("type");
214  assert(tag);
215 
216  assert(tag->GetString() == "PredictionSystNumu2017");
217  delete tag;
218 
219  std::unique_ptr<PredictionSystNumu2017> ret(new PredictionSystNumu2017);
221 
222  // Just enable all norm systs unconditionally. Way too hard to do otherwise
224 
225  delete dir;
226 
227  return ret;
228  }
229 
230  //----------------------------------------------------------------------
231  void PredictionSystNumu2017::SaveTo(TDirectory* dir, const std::string& name) const
232  {
233  PredictionInterp::SaveTo(dir, name);
234 
235  TDirectory* tmp = gDirectory;
236 
237  dir = (TDirectory*)dir->GetDirectory(name.c_str());
238  dir->cd();
239  // Overwrite "IPrediction"
240  TObjString("PredictionSystNumu2017").Write("type");
241 
242  dir->Write();
243  delete dir;
244 
245  tmp->cd();
246  }
247 
248 } // namespace
const XML_Char * name
Definition: expat.h:151
int nBins
Definition: plotROC.py:16
const DummyNumu2017Syst kNumu2017NormSystSig("numu2017normSig","Sig. Norm.",{0,+1}, false)
const DummyNumu2017Syst kNumu2017NormSystBoth("numu2017NormBoth","Sig+Bkg Norm.",{0,+1}, false)
Implements systematic errors by interpolation between shifted templates.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const DummyNumu2017Syst kNumu2017RelMuEScale("numu2017RelMuEScale2017","RelMuEScale",{-2,-1, 0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017MaCCQE_reducedSyst("numu2017MaCCQE_reduced","MaCCQE_reduced",{-2,-1, 0,+1,+2}, false)
std::string systName
What systematic we&#39;re interpolating.
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
const DummyNumu2017Syst kNumu2017LightLevelSyst("numu2017Lightlevel","LightLevel",{-1, 0,+1}, false)
const DummyNumu2017Syst kNumu2017AbsMuEScale("numu2017AbsMuEScale2017","AbsMuEScale",{-2,-1, 0,+1,+2}, false)
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
const DummyNumu2017Syst kNumu2017PPFXPC03Syst("numu2017ppfx_pc03","PPFXPC03",{-2,-1, 0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017RPAShapeRESSyst("numu2017RPAShapeRES","RPAShapeRES",{0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017MaCCRESSyst("numu2017MaCCRES","MaCCRES",{-2,-1, 0,+1,+2}, false)
Float_t tmp
Definition: plot.C:36
const DummyNumu2017Syst kNumu2017MaNCRESSyst("numu2017MaNCRES","MaNCRES",{-2,-1, 0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017MvCCRESSyst("numu2017MvCCRES","MvCCRES",{-2,-1, 0,+1,+2}, false)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
const DummyNumu2017Syst kNumu2017CalibShapeSyst("numu2017CalibShape","CalibShape",{0,+1}, false)
const DummyNumu2017Syst kNumu2017PPFXPC04Syst("numu2017ppfx_pc04","PPFXPC04",{-2,-1, 0,+1,+2}, false)
static std::unique_ptr< PredictionSystNumu2017 > LoadFrom(TDirectory *dir, const std::string &name)
const XML_Char * s
Definition: expat.h:262
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
const DummyNumu2017Syst kNumu2017AbsHadEScale("numu2017AbsHadEScale2017","AbsHadEScale",{-2,-1, 0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017RelCalibSyst("numu2017RelativeCalib","RelativeCalib",{-2, 0,+2}, false)
std::vector< double > shifts
Shift values sampled.
const DummyNumu2017Syst kNumu2017PPFXPC00Syst("numu2017ppfx_pc00","PPFXPC00Sys",{-2,-1, 0,+1,+2}, false)
Loads shifted spectra from files.
std::vector< IPrediction * > preds
const DummyNumu2017Syst kNumu2017CherenkovSyst("numu2017Cherenkov","Cherenkov",{0,+1}, false)
double coeff(double W, int m, int l, bool real)
Oscillation probability calculators.
Definition: Calcs.h:5
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
virtual _IOscCalc * Copy() const =0
Spectrum Predict(osc::IOscCalc *calc) const override
const DummyNumu2017Syst kNumu2017PPFXPC01Syst("numu2017ppfx_pc01","PPFXPC01",{-2,-1, 0,+1,+2}, false)
TDirectory * dir
Definition: macro.C:5
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
const DummyNumu2017Syst kNumu2017NormSystBkg("numu2017normBkg","Bkg. Norm.",{0,+1}, false)
void AddNormSysts(const std::vector< const DummyNumu2017Syst * > &systs)
assert(nhit_max >=nhit_nbins)
const DummyNumu2017Syst kNumu2017AllSmallSysts("numu2017SumSmallXSecJoint2017","SumSmallXSecJoint2017",{-2,-1, 0,+1,+2}, false)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const DummyNumu2017Syst kNumu2017RPAShapeEnhSyst("numu2017RPAShapeenh","RPAShapeEnh",{-2,-1, 0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017AbsCalibSyst("numu2017Calibration","AbsCalib",{-1, 0,+1}, false)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
osc::IOscCalc * fOscOrigin
The oscillation values we assume when evaluating the coefficients.
static void LoadFromBody(TDirectory *dir, PredictionInterp *ret, std::vector< const ISyst * > veto={})
const DummyNumu2017Syst kNumu2017RelHadEScale("numu2017RelHadEScale2017","RelHadEScale",{-2,-1, 0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017AllBeamSysts("numu2017beamTransportComb","beamTransportComb",{-2,-1, 0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017NumuSmallSysts("numu2017SumSmallXSecNumu2017","SumSmallXSecNumu2017",{-2,-1, 0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017DISvnCC1piSyst("numu2017DISvnCC1pi","DISvnCC1pi",{-2,-1, 0,+1,+2}, false)
void SaveTo(TDirectory *dir, const std::string &name) const override
const DummyNumu2017Syst kNumu2017PPFXPC02Syst("numu2017ppfx_pc02","PPFXPC02",{-2,-1, 0,+1,+2}, false)
const DummyNumu2017Syst kNumu2017MECq0ShapeSyst("numu2017MECq0Shape","MECq0Shape",{-2,-1, 0,+1,+2}, false)
std::vector< std::vector< std::vector< Coeffs > > > fits
Indices: [type][histogram bin][shift bin].
enum BeamMode string