PredictionSystJointDemo.cxx
Go to the documentation of this file.
2 
5 #include "3FlavorAna/Cuts/NueCuts2018.h" // kNue2018Binning
8 
9 #include "3FlavorAna/Systs/DummySystStorage.h" // norm systs
10 
11 #include "OscLib/IOscCalc.h"
12 
13 #include "TObjString.h"
14 #include "TROOT.h"
15 
16 #include <iostream>
17 #include <unistd.h>
18 
19 namespace ana
20 {
21  REGISTER_LOADFROM("PredictionSystJointDemo", IPrediction, PredictionSystJointDemo);
22 
23  //----------------------------------------------------------------------
24  // Nue constructor
27  const std::vector<const ISyst*>& systs,
29  bool isFakeData,
30  bool mergePeripheral,
31  const std::string fname,
32  bool minimizeMemory)
33  : fBeam(beam)
34  {
35  gROOT->SetMustClean(false);
36 
37  fOscOrigin = osc->Copy();
38 
39  fSplitBySign = IsFHC();
40 
41  std::cout << "mergePeripheral is "
42  << (mergePeripheral ? "true" : "false") << std::endl;
43 
44  std::cout << "use FakeNDData predictions is "
45  << (isFakeData ? "true" : "false") << std::endl;
46  std::string extrapStr = "Nue_";
47  if (extrap == kProportional) extrapStr += "PropExtrap";
48  if (extrap == kCombo) extrapStr += "ComboExtrap";
49  if (extrap == kNoExtrap) extrapStr += "NoExtrap";
50 
51  if (extrap == kFake) extrapStr += "PropExtrap";
52 
53  std::cout << "Loading nominals..." << std::endl;
54 
55  IPrediction *nom;
56  std::string fnameSyst;
57 
58  fnameSyst = fname;
59  TFile fSyst(fnameSyst.c_str());
60  if(fSyst.IsZombie()){
61  std::cout << "Couldn't open " << fnameSyst << std::endl;
62  abort();
63  }
64 
65  // make PredictionExtendOwning
66 
67  IPrediction* nomCore = ana::LoadFrom<IPrediction>(&fSyst, extrapStr + "/Nominal/NoShift").release();
68  PredictionNoExtrap* nomNoExtrap = ana::LoadFrom<PredictionNoExtrap>(&fSyst, "Nue_NoExtrap/Nominal/NoShift").release();
69 
70  fPredNom = std::make_unique<PredictExtendOwning>(nomCore, nomNoExtrap, kNue2018Binning, true);
71  nom = new PredictExtendOwning(nomCore, nomNoExtrap, kNue2018Binning, true);
72 
73 
74  for(const ISyst* syst: systs){
75 
76  if(syst == &kAna2018NormFHC ||
77  syst == &kAna2018NormRHC) continue;
78 
79  std::string systName = syst->ShortName();
80 
81  std::cout << "Loading " << systName << "..." << std::endl;
82 
83  ShiftedPreds sp;
84  sp.systName = systName;
85 
86  for(int sigma: GetDummyShifts(syst))
87  {
88  sp.shifts.push_back(sigma);
89 
90  if (sigma == 0)
91  {
92  if(syst == &kAnaLightlevelSyst ||
93  syst == &kAnaCherenkovSyst)
94  {
95  std::cout << "\tLooking at systName " << systName << ", so pushing back nomLightLevel " << std::endl;
96  // We have exactly 1 syst for the 2018 that is non-standard is some way
97  // The light level syst uses an alternative nominal. Grab that now.
98 
99  IPrediction* nomLightLevel;
100 
102  &fSyst, extrapStr + "/LightNom/NoShift").release();
103  PredictionNoExtrap* nomLLNoExtrap = ana::LoadFrom<PredictionNoExtrap>(
104  &fSyst, "Nue_NoExtrap/LightNom/NoShift").release();
105 
106  nomLightLevel = new PredictExtendOwning(
107  nomLLCore, nomLLNoExtrap, true);
108 
109  sp.preds.push_back(nomLightLevel);
110  }
111  else
112  sp.preds.push_back(nom);
113  continue;
114  }
115 
116  const std::string sigmaStr = "/"+SigmaToString(sigma);
117 
118  IPrediction* shCore = ana::LoadFrom<PredictionExtrap> (&fSyst, extrapStr + "/"+systName + sigmaStr).release();
119  PredictionNoExtrap* shNoExtrap = ana::LoadFrom<PredictionNoExtrap>(&fSyst, "Nue_NoExtrap/" + systName + sigmaStr).release();
120 
121  sp.preds.push_back(new PredictExtendOwning(shCore, shNoExtrap, kNue2018Binning, true));
122 
123  } // end for sigma
124  // std::cout<<" test find "<<systName<<" "<<syst->find(systName)<<std::endl;
125  fPreds.emplace(syst, sp);
126  } // end for syst
127 
128  // Calculate all the ratios etc
129  InitFits();
130 
131  // AddNormSysts(systs);
132 
133  // Clean up. Don't use if you need to make debug plots
134  if(minimizeMemory){
135  MinimizeMemory();
136  }
137  }
138 
139  //----------------------------------------------------------------------
140  // Numu constructor
143  const std::string beam,
144  int WhichQuant,
145  const std::vector<const ISyst*>& systs,
146  const std::string fname,
147  bool minimizeMemory)
148  : fBeam(beam)
149  {
150  gROOT->SetMustClean(false);
151 
152  fOscOrigin = osc->Copy();
153 
154  fSplitBySign = IsFHC();
155 
156  std::string extrapStr = "Numu_";
157  if (extrap == kNuMu) extrapStr += "Extrap";
158  if (extrap == kNuMuNoExtrap) extrapStr += "NoExtrap";
159  extrapStr += "_";
160 
161  std::string VarName = "Quant"+std::to_string(WhichQuant);
162  if (WhichQuant == -1) VarName = "Quant0";
163 
164  std::cout << "\nLooking at VarName " << VarName << " --> Loading nominals now..." << std::endl;
165 
166  TFile fin(fname.c_str());
167 
168  IPrediction *nom;
169  std::string NomVar = extrapStr+VarName+"/Nominal/NoShift";
171  nom = ana::LoadFrom<IPrediction>(&fin, NomVar).release();
172 
173  for(const ISyst* syst: systs){
174 
175  if(syst == &kAna2018NormFHC ||
176  syst == &kAna2018NormRHC ||
177  syst->ShortName() == "RPAShapeRES2018Test" // || // As we aren't using it in Ana18, and so not in the Pred. files.
178  // -- When using Ana18 Predictions, need to comment out the following XSec Systs.
179  //syst->ShortName() == "RPAShapeenh2019" ||
180  //syst->ShortName() == "RPAShapesupp2019" ||
181  //syst->ShortName() == "RPAShapeRES2019" ||
182  //syst->ShortName() == "MECShape2018RPAFixNu" ||
183  //syst->ShortName() == "MECShape2018RPAFixAntiNu"
184  ) continue;
185 
186  std::string systName = syst->ShortName();
187 
188  std::cout << "Loading " << systName << "..." << std::endl;
189 
190  ShiftedPreds sp;
191  sp.systName = systName;
192 
193  for(int sigma: GetDummyShifts(syst)){
194  sp.shifts.push_back(sigma);
195 
196  if (sigma == 0){
197  if (systName == "Lightlevel" || systName == "Cherenkov") {
198  std::cout << "\tLooking at systName " << systName << ", so pushing back nomLightLevel " << std::endl;
199  // --- The light level syst uses an alternative nominal. Grab that now.
200  std::string LNomVar = extrapStr+VarName+"/LightNom/NoShift/"+VarName;
201  IPrediction* nomLightLevel = ana::LoadFrom<IPrediction>(&fin, LNomVar).release();
202 
203  sp.preds.push_back(nomLightLevel);
204  } else {
205  sp.preds.push_back(nom);
206  }
207  continue;
208  }
209 
210  // const std::string sigmaStr = "/"+SigmaToString(sigma)+"/";
211  const std::string sigmaStr = "/"+SigmaToString(sigma);
212 
213  // auto shCore = ana::LoadFrom<PredictionExtrap> (&fin, extrapStr + systName + sigmaStr +VarName).release();
214  auto shCore = ana::LoadFrom<PredictionExtrap> (&fin, extrapStr + VarName+"/"+systName + sigmaStr).release();
215  sp.preds.push_back(shCore);
216 
217  } // end for sigma
218  fPreds.emplace(syst, sp);
219  } // end for syst
220 
221  // Calculate all the ratios etc
222  InitFits();
223 
224  // AddNormSysts(systs);
225 
226  if(minimizeMemory){
227  MinimizeMemory();
228  }
229  }
230 
231  //----------------------------------------------------------------------
233  {
234  }
235 
236  //----------------------------------------------------------------------
238  AddNormSysts(const std::vector<const ISyst*>& systs)
239  {
240  // DocDb-27230-v7 and DocDb-27821 >v1
241  // leave horn currents uncorrelated
242  // e.g. FHC norm applies shift to fhc,
243  // no shift to rhc, and vice versa
244  const double kNormFHC = 0.0144;
245  const double kNormRHC = 0.0064;
246 
247  // Include syst for wrong beam type to keep the fit happy, but set it to
248  // zero magnitude
249  if(IsFHC()){
250  AddNormSyst(&kAna2018NormFHC, systs, kNormFHC);
251  AddNormSyst(&kAna2018NormRHC, systs, 0);
252  }
253  else{
254  AddNormSyst(&kAna2018NormFHC, systs, 0);
255  AddNormSyst(&kAna2018NormRHC, systs, kNormRHC);
256  }
257  }
258 
259  //----------------------------------------------------------------------
261  AddNormSyst(const ISyst* syst,
262  const std::vector<const ISyst*>& systs,
263  double frac)
264  {
265  // Bail out if we're not actually requested to include this systematic
266  if(std::find(systs.begin(), systs.end(), syst) == systs.end()) return;
267 
268  TH1* h = Predict(fOscOrigin).ToTH1(1);
269  const int nBins = h->GetNbinsX();
270  delete h;
271 
272  ShiftedPreds normsp;
273  normsp.systName = syst->ShortName();
274  normsp.shifts = {0};
275  normsp.preds = {fPredNom.get()};
276  normsp.fits.resize(kNCoeffTypes);
277  if(fSplitBySign) normsp.fitsNubar.resize(kNCoeffTypes);
278 
279  for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
280  normsp.fits[coeff].resize(nBins+2);
281  if(fSplitBySign) normsp.fitsNubar[coeff].resize(nBins+2);
282  for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
283  normsp.fits[coeff][binIdx].emplace_back(0, 0, frac, 1);
284  if(fSplitBySign) normsp.fitsNubar[coeff][binIdx].emplace_back(0, 0, frac, 1);
285  }
286  }
287 
288  normsp.FillRemaps();
289 
290  normsp.nCoeffs = 1;
291 
292  fPreds.emplace(syst, normsp);
293  }
294 
295  //----------------------------------------------------------------------
296  std::unique_ptr<PredictionSystJointDemo> PredictionSystJointDemo::LoadFrom(TDirectory* dir, const std::string& name)
297  {
298  dir = dir->GetDirectory(name.c_str()); // switch to subdir
299  assert(dir);
300 
301  TObjString* tag = (TObjString*)dir->Get("type");
302  assert(tag);
303 
304  assert(tag->GetString() == "PredictionSystJointDemo");
305  delete tag;
306 
307  std::unique_ptr<PredictionSystJointDemo> ret(new PredictionSystJointDemo);
308  LoadFromBody(dir, ret.get(), {&kAna2018NormFHC, &kAna2018NormRHC});
309 
310  const TObjString* beam = (TObjString*)dir->Get("beam");
311  assert(beam);
312  ret->fBeam = beam->GetString();
313 
314  ret->fSplitBySign = ret->IsFHC();
315 
316  // Just enable all norm systs unconditionally. Way too hard to do otherwise
317  // NB fBeam was set above. That's necessary for correctness.
318  // ret->AddNormSysts({&kAna2018NormFHC, &kAna2018NormRHC});
319 
320  delete dir;
321 
322  return ret;
323  }
324 
325  //----------------------------------------------------------------------
326  void PredictionSystJointDemo::SaveTo(TDirectory* dir, const std::string& name) const
327  {
328  PredictionInterp::SaveTo(dir, name);
329 
330  TDirectory* tmp = gDirectory;
331 
332  dir = (TDirectory*)dir->GetDirectory(name.c_str());
333  dir->cd();
334  // Overwrite "IPrediction"
335  TObjString("PredictionSystJointDemo").Write("type");
336 
337  TObjString(fBeam.c_str()).Write("beam");
338 
339  dir->Write();
340  delete dir;
341 
342  tmp->cd();
343  }
344 
345  //----------------------------------------------------------------------
347  {
348  if(sigma == -2) return "MinusTwo";
349  if(sigma == -1) return "MinusOne";
350  if(sigma == 0) return "";
351  if(sigma == +1) return "PlusOne";
352  if(sigma == +2) return "PlusTwo";
353 
354  std::cout << "PredictionSystJointDemo: unknown sigma " << sigma << std::endl;
355  abort();
356  }
357 
358  //----------------------------------------------------------------------
360  {
361  if(fBeam == "fhc" ) return true;
362  if(fBeam == "rhc") return false;
363 
364  std::cout << "PredictionSystJointDemo: Bad beam string: "
365  << fBeam << std::endl;
366  abort();
367  }
368 
369 
370 } // namespace
bool IsFHC() const
Interprets fBeam.
TString fin
Definition: Style.C:24
const XML_Char * name
Definition: expat.h:151
int nBins
Definition: plotROC.py:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
bool mergePeripheral
std::string systName
What systematic we&#39;re interpolating.
Loads shifted spectra from files.
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
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
std::vector< std::vector< std::vector< Coeffs > > > fitsNubar
Will be filled if signs are separated, otherwise not.
void SaveTo(TDirectory *dir, const std::string &name) const override
std::vector< int > GetDummyShifts(const ISyst *s)
Float_t tmp
Definition: plot.C:36
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
const DummyAnaSyst kAnaLightlevelSyst("Lightlevel","Lightlevel")
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
std::string SigmaToString(int sigma) const
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
std::vector< double > shifts
Shift values sampled.
std::vector< IPrediction * > preds
void AddNormSyst(const ISyst *syst, const std::vector< const ISyst * > &systs, double frac)
double frac(double x)
Fractional part.
const Binning kNue2018Binning
Definition: NueCuts2018.cxx:92
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:18
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 DummyAnaSyst kAna2018NormRHC("NormRHC2018","RHC. Norm.")
TDirectory * dir
Definition: macro.C:5
static std::unique_ptr< PredictionSystJointDemo > LoadFrom(TDirectory *dir, const std::string &name)
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
assert(nhit_max >=nhit_nbins)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Prediction that just uses FD MC, with no extrapolation.
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 DummyAnaSyst kAna2018NormFHC("NormFHC2018","FHC. Norm.")
void AddNormSysts(const std::vector< const ISyst * > &systs)
Must set fBeam first.
std::vector< std::vector< std::vector< Coeffs > > > fits
Indices: [type][histogram bin][shift bin].
gm Write()
enum BeamMode string