PredictionInterp.cxx
Go to the documentation of this file.
2 #include "CAFAna/Core/ISyst.h"
4 #include "CAFAna/Core/Ratio.h"
5 #include "CAFAna/Core/Registry.h"
7 
8 #include "TDirectory.h"
9 #include "TH2.h"
10 #include "TMatrixD.h"
11 #include "TObjString.h"
12 #include "TVectorD.h"
13 
14 // For debug plots
15 #include "TGraph.h"
16 #include "TCanvas.h"
17 
18 #include "OscLib/IOscCalc.h"
19 #include "Utilities/func/MathUtil.h"
20 #include "Utilities/func/Stan.h"
22 
23 #include "CAFAna/Core/Loaders.h"
24 
25 #ifndef DARWINBUILD
26 #include <malloc.h>
27 #endif
28 
29 namespace ana
30 {
31  //----------------------------------------------------------------------
33  {
34  // Copy fits into the remapped indexing order. TODO this is very ugly. Best
35  // would be to generate things in this order natively.
36 
37  for(bool nubar: {false, true}){
38  const auto& from = nubar ? fitsNubar : fits;
39  auto& to = nubar ? fitsNubarRemap : fitsRemap;
40 
41  // Allocate the space
42  to.resize(from.size());
43  for(auto& it: to){
44  it.resize(from[0][0].size());
45  for(auto& it2: it) it2.resize(from[0].size(), Coeffs(0, 0, 0, 0));
46  }
47 
48  // And copy with the transposed indices
49  for(unsigned int i = 0; i < to.size(); ++i)
50  for(unsigned int j = 0; j < to[i].size(); ++j)
51  for(unsigned int k = 0; k < to[i][j].size(); ++k)
52  to[i][j][k] = from[i][k][j];
53 
54  } // end for nubar
55  }
56 
57  //----------------------------------------------------------------------
58  PredictionInterp::PredictionInterp(std::vector<const ISyst*> systs,
60  const IPredictionGenerator& predGen,
62  const SystShifts& shiftMC,
63  EMode_t mode)
64  : fOscOrigin(osc ? osc->Copy() : 0),
65  fBinning(Spectrum::Uninitialized()),
66  fSplitBySign(mode == kSplitBySign)
67  {
68  for(const ISyst* syst: systs){
69  ShiftedPreds sp;
70  sp.systName = syst->ShortName();
71 
72  // Genie reweight shifts aren't able to provide +/-3 sigma
73  if(syst->IsGenieReweight())
74  sp.shifts = {-2, -1, 0, +1, +2};
75  else
76  sp.shifts = {-3, -2, -1, 0, +1, +2, +3};
77 
78  for(int sigma: sp.shifts){
79  SystShifts shiftHere = shiftMC;
80  shiftHere.SetShift(syst, sigma);
81  sp.preds.emplace_back(predGen.Generate(loaders, shiftHere).release());
82  }
83 
84  fPreds.emplace(syst, std::move(sp));
85  } // end for syst
86 
87  fPredNom = predGen.Generate(loaders, shiftMC);
88  }
89 
90  //----------------------------------------------------------------------
92  {
93  }
94 
95  //----------------------------------------------------------------------
96  std::vector<std::vector<PredictionInterp::Coeffs>> PredictionInterp::
97  FitRatios(const std::vector<double>& shifts,
98  const std::vector<Eigen::ArrayXd>& ratios) const
99  {
100  if(ratios.size() < 2){
101  std::cout << "PredictionInterp::FitRatios(): ratios.size() = " << ratios.size() << " - how did that happen?" << std::endl;
102 
103  abort();
104  }
105 
106  assert(shifts.size() == ratios.size());
107 
108  std::vector<std::vector<Coeffs>> ret;
109 
110  const int binMax = ratios[0].size();
111 
112  for(int binIdx = 0; binIdx < binMax; ++binIdx){
113  ret.push_back({});
114 
115  // This is cubic interpolation. For each adjacent set of four points we
116  // determine coefficients for a cubic which will be the curve between the
117  // center two. We constrain the function to match the two center points
118  // and to have the right mean gradient at them. This causes this patch to
119  // match smoothly with the next one along. The resulting function is
120  // continuous and first and second differentiable. At the ends of the
121  // range we fit a quadratic instead with only one constraint on the
122  // slope. The coordinate conventions are that point y1 sits at x=0 and y2
123  // at x=1. The matrices are simply the inverses of writing out the
124  // constraints expressed above.
125 
126  // Special-case for linear interpolation
127  if(ratios.size() == 2){
128  const double y0 = ratios[0][binIdx];
129  const double y1 = ratios[1][binIdx];
130 
131  ret.back().emplace_back(0, 0, y1-y0, y0);
132  continue;
133  }
134 
135  {
136  const double y1 = ratios[0][binIdx];
137  const double y2 = ratios[1][binIdx];
138  const double y3 = ratios[2][binIdx];
139  const double v[3] = {y1, y2, (y3-y1)/2};
140  const double m[9] = { 1, -1, 1,
141  -2, 2, -1,
142  1, 0, 0};
143  const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
144  ret.back().emplace_back(0, res(0), res(1), res(2));
145  }
146 
147  // We're assuming here that the shifts are separated by exactly 1 sigma.
148  for(unsigned int shiftIdx = 1; shiftIdx < ratios.size()-2; ++shiftIdx){
149  const double y0 = ratios[shiftIdx-1][binIdx];
150  const double y1 = ratios[shiftIdx ][binIdx];
151  const double y2 = ratios[shiftIdx+1][binIdx];
152  const double y3 = ratios[shiftIdx+2][binIdx];
153 
154  const double v[4] = {y1, y2, (y2-y0)/2, (y3-y1)/2};
155  const double m[16] = { 2, -2, 1, 1,
156  -3, 3, -2, -1,
157  0, 0, 1, 0,
158  1, 0, 0, 0};
159  const TVectorD res = TMatrixD(4, 4, m) * TVectorD(4, v);
160  ret.back().emplace_back(res(0), res(1), res(2), res(3));
161  } // end for shiftIdx
162 
163  {
164  const int N = ratios.size()-3;
165  const double y0 = ratios[N ][binIdx];
166  const double y1 = ratios[N+1][binIdx];
167  const double y2 = ratios[N+2][binIdx];
168  const double v[3] = {y1, y2, (y2-y0)/2};
169  const double m[9] = {-1, 1, -1,
170  0, 0, 1,
171  1, 0, 0};
172  const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
173  ret.back().emplace_back(0, res(0), res(1), res(2));
174  }
175  } // end for binIdx
176 
177  double stride = -1;
178  for(unsigned int i = 0; i < shifts.size()-1; ++i){
179  const double newStride = shifts[i+1]-shifts[i];
180  assert((stride < 0 || fabs(stride-newStride) < 1e-3) &&
181  "Variably-spaced syst templates are unsupported");
182  stride = newStride;
183  }
184 
185  // If the stride is actually not 1, need to rescale all the coefficients
186  for(std::vector<Coeffs>& cs: ret)
187  for(Coeffs& c: cs){
188  c = Coeffs(c.a/util::cube(stride),
189  c.b/util::sqr(stride),
190  c.c/stride,
191  c.d);}
192  return ret;
193  }
194 
195  //----------------------------------------------------------------------
196  std::vector<std::vector<PredictionInterp::Coeffs>> PredictionInterp::
197  FitComponent(const std::vector<double>& shifts,
198  const std::vector<IPrediction*>& preds,
199  Flavors::Flavors_t flav,
202  const std::string& systName) const
203  {
204  IPrediction* pNom = 0;
205  for(unsigned int i = 0; i < shifts.size(); ++i){
206  if(shifts[i] == 0) {pNom = preds[i]; };
207  }
208  assert(pNom);
209 
210  // Do it this way rather than via fPredNom so that systematics evaluated
211  // relative to some alternate nominal (eg Birks C where the appropriate
212  // nominal is no-rock) can work.
213  const Spectrum nom = pNom->PredictComponent(fOscOrigin,
214  flav, curr, sign);
215 
216  std::vector<Eigen::ArrayXd> ratios;
217  ratios.reserve(preds.size());
218  for(auto& p: preds){
219  ratios.emplace_back(Ratio(p->PredictComponent(fOscOrigin,
220  flav, curr, sign),
221  nom).GetEigen());
222 
223  // Check none of the ratio values is crazy
224  Eigen::ArrayXd& r = ratios.back();
225  for(int i = 0; i < r.size(); ++i){
226  if(r[i] > 500){
227  std::cout << "PredictionInterp: WARNING, ratio in bin "
228  << i << " for " << shifts[&p-&preds.front()]
229  << " sigma shift of " << systName
230  << " is " << r[i] << ". Ignoring." << std::endl;
231  r[i] = 1;
232  }
233  }
234  }
235 
236  return FitRatios(shifts, ratios);
237  }
238 
239  //----------------------------------------------------------------------
240  std::vector<const ISyst*> PredictionInterp::GetAllSysts() const
241  {
242  std::vector<const ISyst*> allsysts;
243  for (const auto &p : fPreds)
244  allsysts.push_back(p.first);
245 
246  return allsysts;
247  }
248 
249  //----------------------------------------------------------------------
251  std::vector<std::vector<std::vector<Coeffs>>>& fits,
252  Sign::Sign_t sign) const
253  {
254  fits.resize(kNCoeffTypes);
255 
259 
261 
263  }
264 
265  //----------------------------------------------------------------------
267  {
268  // No systs
269  if(fPreds.empty()){
270  if(fBinning.POT() > 0 || fBinning.Livetime() > 0) return;
271  }
272  // Already initialized
273  else if(!fPreds.empty() && !fPreds.begin()->second.fits.empty()) return;
274 
275  for(auto& it: fPreds){
276  ShiftedPreds& sp = it.second;
277 
278  if(fSplitBySign){
279  InitFitsHelper(sp, sp.fits, Sign::kNu);
281  }
282  else{
284  }
285  sp.nCoeffs = sp.fits[0][0].size();
286 
287  // Copy the outputs into the remapped indexing order. TODO this is very
288  // ugly. Best would be to generate things in this order natively.
289  sp.FillRemaps();
290  }
291 
292  // Predict something, anything, so that we can know what binning to use
293  fBinning = fPredNom->Predict(fOscOrigin);
294  fBinning.Clear();
295  }
296 
297  //----------------------------------------------------------------------
299  fOscOrigin = oscSeed->Copy();
300  for(auto& it: fPreds) it.second.fits.clear();
301  InitFits();
302  }
303 
304  //----------------------------------------------------------------------
306  {
307  return fPredNom->Predict(calc);
308  }
309 
310  //----------------------------------------------------------------------
312  {
313  return fPredNom->Predict(calc);
314  }
315 
316  //----------------------------------------------------------------------
318  Flavors::Flavors_t flav,
320  Sign::Sign_t sign) const
321  {
322  return fPredNom->PredictComponent(calc, flav, curr, sign);
323  }
324 
325  //----------------------------------------------------------------------
327  Flavors::Flavors_t flav,
329  Sign::Sign_t sign) const
330  {
331  return fPredNom->PredictComponent(calc, flav, curr, sign);
332  }
333 
334  //----------------------------------------------------------------------
336  const SystShifts& shift) const
337  {
338  InitFits();
339 
340  return PredictComponentSyst(calc, shift,
343  Sign::kBoth);
344  }
345 
346  //----------------------------------------------------------------------
348  const SystShifts& shift) const
349  {
350  InitFits();
351 
352  return PredictComponentSyst(calc,
353  shift,
356  Sign::kBoth);
357  }
358 
359  //----------------------------------------------------------------------
362  bool nubar,
363  const SystShifts& shift) const
364  {
365  if(s.HasStan() && (!shift.IsNominal() && !shift.HasAnyStan())){
366  std::cout << "PredictionInterp used with Stan oscillations but non-Stan syst shifts. Is that what you expected?" << std::endl;
367  std::cout << " Systs set: " << std::endl;
368  for (const auto & syst : shift.ActiveSysts())
369  std::cout << " " << syst->ShortName() << " = " << shift.GetShift(syst) << std::endl;
370  }
371 
372  if(s.HasStan() || shift.HasAnyStan()){
373  // Need the second case for the NC component
375 
376  ShiftBins(vec.size(), vec.data(), type, nubar, shift);
377  return Spectrum(std::move(vec), HistAxis(s.GetLabels(), s.GetBinnings()), s.POT(), s.Livetime());
378  }
379  else{
380  Eigen::ArrayXd vec = s.GetEigen(s.POT());
381  ShiftBins(vec.size(), vec.data(), type, nubar, shift);
382  return Spectrum(std::move(vec), HistAxis(s.GetLabels(), s.GetBinnings()), s.POT(), s.Livetime());
383  }
384  }
385 
386  //----------------------------------------------------------------------
387  template <typename T>
388  void PredictionInterp::ShiftBins(unsigned int N,
389  T* arr,
391  bool nubar,
392  const SystShifts& shift) const
393  {
394  static_assert(std::is_same_v<T, double> ||
395  std::is_same_v<T, stan::math::var>,
396  "PredictionInterp::ShiftBins() can only be called using doubles or stan::math::vars");
397  if(nubar) assert(fSplitBySign);
398 
399 
400  if(!std::is_same_v<T, stan::math::var> && shift.HasAnyStan()){
401  std::cout << "PredictionInterp: stan shifts on non-stan spectrum, something is wrong" << std::endl;
402  abort();
403  }
404 
405  T corr[N];
406  for(unsigned int i = 0; i < N; ++i) corr[i] = 1;
407 
408  for(auto& it: fPreds){
409  const ISyst* syst = it.first;
410  const ShiftedPreds& sp = it.second;
411 
412  T x = shift.GetShift<T>(syst);
413 
414  // need to actually do the calculation for the autodiff version
415  // to make sure the gradient is computed correctly
416  if(x == 0 && (!std::is_same_v<T, stan::math::var> || !shift.HasStan(syst))) continue;
417 
418  int shiftBin = (util::GetValAs<double>(x) - sp.shifts[0])/sp.Stride();
419  shiftBin = std::max(0, shiftBin);
420  shiftBin = std::min(shiftBin, sp.nCoeffs - 1);
421 
422  const Coeffs* fits = nubar ? &sp.fitsNubarRemap[type][shiftBin].front() : &sp.fitsRemap[type][shiftBin].front();
423 
424  x -= sp.shifts[shiftBin];
425 
426  const T x_cube = util::cube(x);
427  const T x_sqr = util::sqr(x);
428 
429  for(unsigned int n = 0; n < N; ++n){
430  // Uncomment to debug crashes in this function
431  // assert(type < fits.size());
432  // assert(n < sp.fits[type].size());
433  // assert(shiftBin < int(sp.fits[type][n].size()));
434 
435  const Coeffs& f = fits[n];
436 
437  corr[n] *= f.a*x_cube + f.b*x_sqr + f.c*x + f.d;
438  } // end for n
439  } // end for syst
440 
441  for(unsigned int n = 0; n < N; ++n){
442  // std::max() doesn't work with stan::math::var
443  arr[n] *= (corr[n] > 0.) ? corr[n] : 0.;
444  }
445 
446  }
447 
448  //----------------------------------------------------------------------
450  const TMD5* hash,
451  const SystShifts& shift,
452  Flavors::Flavors_t flav,
455  CoeffsType type) const
456  {
457  return _ShiftedComponent(calc, hash, shift, flav, curr, sign, type);
458  }
459 
460  //----------------------------------------------------------------------
462  const TMD5* hash,
463  const SystShifts& shift,
464  Flavors::Flavors_t flav,
467  CoeffsType type) const
468  {
469  return _ShiftedComponent(calc, hash, shift, flav, curr, sign, type);
470  }
471 
472  //----------------------------------------------------------------------
473  template<typename T>
475  const TMD5* hash,
476  const SystShifts& shift,
477  Flavors::Flavors_t flav,
480  CoeffsType type) const
481  {
483  "PredictionInterp::ShiftedComponent() can only be called using doubles or stan::math::vars");
484 
485  if(fSplitBySign && sign == Sign::kBoth){
486  return (ShiftedComponent(calc, hash, shift, flav, curr, Sign::kAntiNu, type)+
487  ShiftedComponent(calc, hash, shift, flav, curr, Sign::kNu, type));
488  }
489 
490  // Must be the base case of the recursion to use the cache. Otherwise we
491  // can cache systematically shifted versions of our children, which is
492  // wrong. Also, some calculators won't hash themselves.
493  // Moreover, caching is not going to work with stan::math::vars
494  // since they get reset every time Stan's log_prob() is called.
495  const bool canCache = (hash != 0) && !std::is_same<T, stan::math::var>::value;
496 
497  const Key_t key = {flav, curr, sign};
498  auto it = fNomCache->find(key);
499 
500  // Should the interpolation use the nubar fits?
501  const bool nubar = (fSplitBySign && sign == Sign::kAntiNu);
502 
503  // We have the nominal for this exact combination of flav, curr, sign, calc
504  // stored. Shift it and return.
505  if(canCache && it != fNomCache->end() && it->second.hash == *hash){
506  return ShiftSpectrum(it->second.nom, type, nubar, shift);
507  }
508 
509  // We need to compute the nominal again for whatever reason
510  const Spectrum nom = fPredNom->PredictComponent(calc, flav, curr, sign);
511 
512  if(canCache){
513  // Insert into the cache if not already there, or update if there but
514  // with old oscillation parameters.
515  if(it == fNomCache->end())
516  fNomCache->emplace(key, Val_t({*hash, nom}));
517  else
518  it->second = {*hash, nom};
519  }
520 
521  return ShiftSpectrum(nom, type, nubar, shift);
522  }
523 
524  //----------------------------------------------------------------------
525  template<typename T>
527  const SystShifts& shift,
528  Flavors::Flavors_t flav,
530  Sign::Sign_t sign) const
531  {
532  InitFits();
533 
535  ret.Clear();
536 
537  // Check that we're able to handle all the systs we were passed
538  for(const ISyst* syst: shift.ActiveSysts()){
539  if(fPreds.find(syst) == fPreds.end()){
540  std::cerr << "This PredictionInterp is not set up to handle the requested systematic: " << syst->ShortName() << std::endl;
541  abort();
542  }
543  } // end for syst
544 
545 
546  const TMD5* hash = calc ? calc->GetParamsHash() : 0;
547 
548  if(curr & Current::kCC){
549  if(flav & Flavors::kNuEToNuE) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuE, Current::kCC, sign, kNueSurv);
550  if(flav & Flavors::kNuEToNuMu) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuMu, Current::kCC, sign, kOther );
551  if(flav & Flavors::kNuEToNuTau) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuTau, Current::kCC, sign, kOther );
552 
553  if(flav & Flavors::kNuMuToNuE) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuE, Current::kCC, sign, kNueApp );
554  if(flav & Flavors::kNuMuToNuMu) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuMu, Current::kCC, sign, kNumuSurv);
555  if(flav & Flavors::kNuMuToNuTau) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuTau, Current::kCC, sign, kOther );
556  }
557  if(curr & Current::kNC){
558  assert(flav == Flavors::kAll); // Don't know how to calculate anything else
559 
560  ret += ShiftedComponent(calc, hash, shift, Flavors::kAll, Current::kNC, sign, kNC);
561  }
562 
563  delete hash;
564 
565  return ret;
566  }
567 
568  //----------------------------------------------------------------------
569  // can't template these directly since the interface isn't templated
571  const SystShifts& shift,
572  Flavors::Flavors_t flav,
574  Sign::Sign_t sign) const
575  {
576  return _PredictComponentSyst(calc, shift, flav, curr, sign);
577  }
578 
579  //----------------------------------------------------------------------
581  const SystShifts& shift,
582  Flavors::Flavors_t flav,
584  Sign::Sign_t sign) const
585  {
586  return _PredictComponentSyst(calc, shift, flav, curr, sign);
587  }
588 
589  void PredictionInterp::SaveTo(TDirectory* dir, const std::string& name) const
590  {
591  InitFits();
592 
593  TDirectory* tmp = gDirectory;
594 
595  dir = dir->mkdir(name.c_str()); // switch to subdir
596  dir->cd();
597 
598  TObjString("PredictionInterp").Write("type");
599 
600  fPredNom->SaveTo(dir, "pred_nom");
601 
602  for(auto& it: fPreds){
603  const ShiftedPreds& sp = it.second;
604 
605  for(unsigned int i = 0; i < sp.shifts.size(); ++i){
606  if(!sp.preds[i]){
607  std::cout << "Can't save a PredictionInterp after MinimizeMemory()" << std::endl;
608  abort();
609  }
610 
611  sp.preds[i]->SaveTo(dir, TString::Format("pred_%s_%+d",
612  sp.systName.c_str(),
613  int(sp.shifts[i])).Data());
614  } // end for i
615  } // end for it
616 
617  ana::SaveTo(*fOscOrigin, dir, "osc_origin");
618 
619  if(!fPreds.empty()){
620  TH1F hSystNames("syst_names", ";Syst names", fPreds.size(), 0, fPreds.size());
621  int binIdx = 1;
622  for(auto& it: fPreds){
623  hSystNames.GetXaxis()->SetBinLabel(binIdx++, it.second.systName.c_str());
624  }
625  dir->cd();
626  hSystNames.Write("syst_names");
627  }
628 
629  TObjString split_sign = fSplitBySign ? "yes" : "no";
630  dir->cd();
631  split_sign.Write("split_sign");
632 
633  dir->Write();
634  delete dir;
635 
636  tmp->cd();
637  }
638 
639  //----------------------------------------------------------------------
640  std::unique_ptr<PredictionInterp> PredictionInterp::LoadFrom(TDirectory* dir, const std::string& name)
641  {
642  dir = dir->GetDirectory(name.c_str()); // switch to subdir
643  assert(dir);
644 
645  TObjString* tag = (TObjString*)dir->Get("type");
646  assert(tag);
647  assert(tag->GetString() == "PredictionInterp");
648 
649  std::unique_ptr<PredictionInterp> ret(new PredictionInterp);
650 
651  LoadFromBody(dir, ret.get());
652 
653  TObjString* split_sign = (TObjString*)dir->Get("split_sign");
654  // Can be missing from old files
655  ret->fSplitBySign = (split_sign && split_sign->String() == "yes");
656 
657  delete dir;
658 
659  return ret;
660  }
661 
662  //----------------------------------------------------------------------
664  std::vector<const ISyst*> veto)
665  {
666  ret->fPredNom = ana::LoadFrom<IPrediction>(dir, "pred_nom");
667 
668  TH1* hSystNames = (TH1*)dir->Get("syst_names");
669  if(hSystNames){
670  for(int systIdx = 0; systIdx < hSystNames->GetNbinsX(); ++systIdx){
671  ShiftedPreds sp;
672  sp.systName = hSystNames->GetXaxis()->GetBinLabel(systIdx+1);
673 
675 
676  if(std::find(veto.begin(), veto.end(), syst) != veto.end()) continue;
677 
678  for(int shift = -3; shift <= +3; ++shift){
679  const std::string subname = TString::Format("pred_%s_%+d", sp.systName.c_str(), shift).Data();
680  TDirectory *preddir = dir->GetDirectory(subname.c_str());
681  if(!preddir) continue; // Can happen for genie systs
682  delete preddir;
683 
684  sp.shifts.push_back(shift);
685  sp.preds.emplace_back(ana::LoadFrom<IPrediction>(dir, subname).release());
686  } // end for shift
687 
688  ret->fPreds.emplace(syst, std::move(sp));
689  } // end for systIdx
690  } // end if hSystNames
691 
692  ret->fOscOrigin = ana::LoadFrom<osc::IOscCalc>(dir, "osc_origin").release();
693  }
694 
695  //----------------------------------------------------------------------
697  {
698  std::set<IPrediction*> todel;
699  for(auto& it: fPreds){
700  std::vector<IPrediction*>& preds = it.second.preds;
701  for(unsigned int i = 0; i < preds.size(); ++i){
702  if(preds[i] != fPredNom.get()){
703  todel.insert(preds[i]);
704  preds[i] = 0;
705  }
706  }
707  }
708 
709  for(IPrediction* p: todel) delete p;
710 
711  // We probably just freed up a lot of memory, but malloc by default hangs
712  // on to all of it as cache.
713  #ifndef DARWINBUILD
714  malloc_trim(0);
715  #endif
716  }
717 
718  //----------------------------------------------------------------------
720  osc::IOscCalc* calc,
721  Flavors::Flavors_t flav,
723  Sign::Sign_t sign) const
724  {
725  InitFits();
726 
727  auto it = fPreds.find(syst);
728  if(it == fPreds.end()){
729  std::cout << "PredictionInterp::DebugPlots(): "
730  << syst->ShortName() << " not found" << std::endl;
731  return;
732  }
733 
734  std::unique_ptr<TH1> nom(fPredNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
735  const int nbins = nom->GetNbinsX();
736 
737  TGraph* curves[nbins];
738  TGraph* points[nbins];
739 
740  for(int i = 0; i <= 80; ++i){
741  const double x = .1*i-4;
742  const SystShifts ss(it->first, x);
743  std::unique_ptr<TH1> h(PredictComponentSyst(calc, ss, flav, curr, sign).ToTH1(18e20));
744 
745  for(int bin = 0; bin < nbins; ++bin){
746  if(i == 0){
747  curves[bin] = new TGraph;
748  points[bin] = new TGraph;
749  }
750 
751  const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
752 
753  if(!std::isnan(ratio)) curves[bin]->SetPoint(curves[bin]->GetN(), x, ratio);
754  else curves[bin]->SetPoint(curves[bin]->GetN(), x, 1);
755  } // end for bin
756  } // end for i (x)
757 
758  // As elswhere, to allow BirksC etc that need a different nominal to plot
759  // right.
760  IPrediction* pNom = 0;
761  for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
762  if(it->second.shifts[shiftIdx] == 0) pNom = it->second.preds[shiftIdx];;
763  }
764  if(pNom){ // if not, probably MinimizeMemory() was called
765  std::unique_ptr<TH1> hnom(pNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
766 
767  for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
768  if(!it->second.preds[shiftIdx]) continue; // Probably MinimizeMemory()
769  std::unique_ptr<TH1> h;
770  h = std::move(std::unique_ptr<TH1>(it->second.preds[shiftIdx]->PredictComponent(calc, flav, curr, sign).ToTH1(18e20)));
771 
772  for(int bin = 0; bin < nbins; ++bin){
773  const double ratio = h->GetBinContent(bin+1)/hnom->GetBinContent(bin+1);
774  if(!std::isnan(ratio)) points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], ratio);
775  else points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], 1);
776  }
777  } // end for shiftIdx
778  } // end if pNom
779 
780 
781  int nx = int(sqrt(nbins));
782  int ny = int(sqrt(nbins));
783  if(nx*ny < nbins) ++nx;
784  if(nx*ny < nbins) ++ny;
785 
786  TCanvas* c = new TCanvas;
787  c->Divide(nx, ny);
788 
789  for(int bin = 0; bin < nbins; ++bin){
790  c->cd(bin+1);
791  (new TH2F("",
792  TString::Format("%s %g < %s < %g;Shift;Ratio",
793  it->second.systName.c_str(),
794  nom->GetXaxis()->GetBinLowEdge(bin+1),
795  nom->GetXaxis()->GetTitle(),
796  nom->GetXaxis()->GetBinUpEdge(bin+1)),
797  100, -4, +4, 100, .5, 1.5))->Draw();
798  curves[bin]->Draw("l same");
799  points[bin]->SetMarkerStyle(kFullDotMedium);
800  points[bin]->Draw("p same");
801  } // end for bin
802 
803  c->cd(0);
804  }
805 
806  //----------------------------------------------------------------------
808  const std::string& savePattern,
809  Flavors::Flavors_t flav,
811  Sign::Sign_t sign) const
812  {
813  for(auto& it: fPreds){
814  DebugPlot(it.first, calc, flav, curr, sign);
815 
816  if(!savePattern.empty()){
817  assert(savePattern.find("%s") != std::string::npos);
818  gPad->Print(TString::Format(savePattern.c_str(), it.second.systName.c_str()).Data());
819  }
820  } // end for it
821  }
822 
823  //----------------------------------------------------------------------
825  osc::IOscCalc* calc,
826  Flavors::Flavors_t flav,
828  Sign::Sign_t sign) const
829  {
830  InitFits();
831 
832  std::unique_ptr<TH1> nom(fPredNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
833  const int nbins = nom->GetNbinsX();
834 
835  TH2* h2 = new TH2F("", (syst->LatexName()+";;#sigma").c_str(),
836  nbins, nom->GetXaxis()->GetXmin(), nom->GetXaxis()->GetXmax(),
837  80, -4, +4);
838  h2->GetXaxis()->SetTitle(nom->GetXaxis()->GetTitle());
839 
840  for(int i = 1; i <= 80; ++i){
841  const double y = h2->GetYaxis()->GetBinCenter(i);
842  const SystShifts ss(syst, y);
843  std::unique_ptr<TH1> h(PredictComponentSyst(calc, ss, flav, curr, sign).ToTH1(18e20));
844 
845  for(int bin = 0; bin < nbins; ++bin){
846  const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
847 
848  if(!isnan(ratio) && !isinf(ratio))
849  h2->Fill(h2->GetXaxis()->GetBinCenter(bin), y, ratio);
850  } // end for bin
851  } // end for i (x)
852 
853  h2->Draw("colz");
854  h2->SetMinimum(0.5);
855  h2->SetMaximum(1.5);
856  }
857 
858  //----------------------------------------------------------------------
860  const std::string& savePattern,
861  Flavors::Flavors_t flav,
863  Sign::Sign_t sign) const
864  {
865  InitFits();
866 
867  for(auto& it: fPreds){
868  new TCanvas;
869  DebugPlotColz(it.first, calc, flav, curr, sign);
870 
871  if(!savePattern.empty()){
872  assert(savePattern.find("%s") != std::string::npos);
873  gPad->Print(TString::Format(savePattern.c_str(), it.second.systName.c_str()).Data());
874  }
875  } // end for it
876  }
877 
878 } // namespace
T max(const caf::Proxy< T > &a, T b)
void ShiftBins(unsigned int N, T *arr, CoeffsType type, bool nubar, const SystShifts &shift) const
Helper for ShiftSpectrum.
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:268
const XML_Char * name
Definition: expat.h:151
tree Draw("slc.nhit")
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
bool HasAnyStan() const
Definition: SystShifts.h:64
Implements systematic errors by interpolation between shifted templates.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
set< int >::iterator it
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
std::string systName
What systematic we&#39;re interpolating.
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
bool IsNominal() const
Definition: SystShifts.h:43
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:209
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
ThreadLocal< std::map< Key_t, Val_t > > fNomCache
std::vector< std::vector< std::vector< Coeffs > > > fitsNubarRemap
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
(&#39;beam &#39;)
Definition: IPrediction.h:15
const char * p
Definition: xmltok.h:285
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ratio(TH1 *h1, TH1 *h2)
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:200
Float_t ss
Definition: plot.C:24
std::vector< std::vector< std::vector< Coeffs > > > fitsNubar
Will be filled if signs are separated, otherwise not.
void SetOscSeed(osc::IOscCalc *oscSeed)
OStream cerr
Definition: OStream.cxx:7
void Clear()
Definition: Spectrum.cxx:433
double corr(TGraph *g, double thresh)
Float_t tmp
Definition: plot.C:36
T cube(T x)
More efficient cube function than pow(x,3)
Definition: MathUtil.h:26
void DebugPlotsColz(osc::IOscCalc *calc, const std::string &savePattern="", Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
osc::OscCalcDumb calc
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
TVectorT< double > TVectorD
Definition: Utilities.h:18
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
::xsd::cxx::tree::type type
Definition: Database.h:110
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
const XML_Char * s
Definition: expat.h:262
const int nbins
Definition: cellShifts.C:15
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Charged-current interactions.
Definition: IPrediction.h:39
Interactions of both types.
Definition: IPrediction.h:42
const XML_Char int const XML_Char * value
Definition: expat.h:331
std::vector< double > shifts
Shift values sampled.
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
std::vector< std::vector< Coeffs > > FitComponent(const std::vector< double > &shifts, const std::vector< IPrediction * > &preds, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign, const std::string &systName) const
Find coefficients describing the ratios from this component.
T GetShift(const ISyst *syst) const
std::vector< IPrediction * > preds
const double j
Definition: BetheBloch.cxx:29
Eigen::VectorXd vec
Spectrum fBinning
Dummy spectrum to provide binning.
Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
TMatrixT< double > TMatrixD
Definition: Utilities.h:16
bool HasStan(const ISyst *s) const
Definition: SystShifts.h:63
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:28
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:231
Oscillation probability calculators.
Definition: Calcs.h:5
virtual TMD5 * GetParamsHash() const
Use to check two calculators are in the same state.
Definition: IOscCalc.h:34
double sigma(TH1F *hist, double percentile)
void DebugPlots(osc::IOscCalc *calc, const std::string &savePattern="", Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
Represent the ratio between two spectra.
Definition: Ratio.h:8
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
Taus, numu appearance.
TH1F * h2
Definition: plot.C:45
virtual _IOscCalc * Copy() const =0
Spectrum Predict(osc::IOscCalc *calc) const override
Spectrum _ShiftedComponent(osc::_IOscCalc< T > *calc, const TMD5 *hash, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign, CoeffsType type) const
Templated helper for ShiftedComponent.
virtual const std::string & LatexName() const final
The name used on plots (ROOT&#39;s TLatex syntax)
Definition: ISyst.h:30
void DebugPlotColz(const ISyst *syst, osc::IOscCalc *calc, Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
std::vector< std::vector< std::vector< Coeffs > > > fitsRemap
Eigen::Array< stan::math::var, Eigen::Dynamic, 1 > ArrayXstan
Definition: Hist.h:20
TDirectory * dir
Definition: macro.C:5
virtual std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const =0
std::vector< Loaders * > loaders
Definition: syst_header.h:386
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
std::vector< std::string > fits
Neutral-current interactions.
Definition: IPrediction.h:40
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
std::vector< const ISyst * > GetAllSysts() const
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::vector< const ISyst * > ActiveSysts() const
Definition: SystShifts.cxx:207
void DebugPlot(const ISyst *syst, osc::IOscCalc *calc, Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
void InitFitsHelper(ShiftedPreds &sp, std::vector< std::vector< std::vector< Coeffs >>> &fits, Sign::Sign_t sign) const
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
double T
Definition: Xdiff_gwt.C:5
Given loaders and an MC shift, Generate() generates an IPrediction.
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={})
T min(const caf::Proxy< T > &a, T b)
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:267
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
Float_t e
Definition: plot.C:35
std::vector< std::vector< Coeffs > > FitRatios(const std::vector< double > &shifts, const std::vector< Eigen::ArrayXd > &ratios) const
Find coefficients describing this set of shifts.
bool HasStan() const
Definition: Spectrum.h:198
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:234
Spectrum ShiftedComponent(osc::IOscCalc *calc, const TMD5 *hash, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign, CoeffsType type) const
Helper for PredictComponentSyst.
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
Spectrum _PredictComponentSyst(osc::_IOscCalc< T > *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Templated helper for PredictComponentSyst.
Spectrum ShiftSpectrum(const Spectrum &s, CoeffsType type, bool nubar, const SystShifts &shift) const
def sign(x)
Definition: canMan.py:197
static const T * ShortNameToPtr(const std::string &s, bool allowFail=false)
Definition: Registry.cxx:60
std::vector< std::vector< std::vector< Coeffs > > > fits
Indices: [type][histogram bin][shift bin].
const Eigen::ArrayXstan & GetEigenStan() const
Definition: Spectrum.h:201