PredictionInterp.cxx
Go to the documentation of this file.
2 #include "CAFAna/Core/ISyst.h"
5 #include "CAFAna/Core/Ratio.h"
6 #include "CAFAna/Core/Registry.h"
8
9 #include "TDirectory.h"
10 #include "TH2.h"
11 #include "TMatrixD.h"
12 #include "TObjString.h"
13 #include "TVectorD.h"
14
15 // For debug plots
16 #include "TGraph.h"
17 #include "TCanvas.h"
18
19 #include "OscLib/IOscCalc.h"
20 #include "Utilities/func/MathUtil.h"
21 #include "Utilities/func/Stan.h"
23
25
26 #include <malloc.h>
27 namespace ana
28 {
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);
82  }
83
84  fPreds.emplace(syst, std::move(sp));
85  } // end for syst
86
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  }
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
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  {
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);
686  } // end for shift
687
688  ret->fPreds.emplace(syst, std::move(sp));
689  } // end for systIdx
690  } // end if hSystNames
691
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  malloc_trim(0);
714  }
715
716  //----------------------------------------------------------------------
718  osc::IOscCalc* calc,
719  Flavors::Flavors_t flav,
721  Sign::Sign_t sign) const
722  {
723  InitFits();
724
725  auto it = fPreds.find(syst);
726  if(it == fPreds.end()){
727  std::cout << "PredictionInterp::DebugPlots(): "
729  return;
730  }
731
732  std::unique_ptr<TH1> nom(fPredNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
733  const int nbins = nom->GetNbinsX();
734
735  TGraph* curves[nbins];
736  TGraph* points[nbins];
737
738  for(int i = 0; i <= 80; ++i){
739  const double x = .1*i-4;
740  const SystShifts ss(it->first, x);
741  std::unique_ptr<TH1> h(PredictComponentSyst(calc, ss, flav, curr, sign).ToTH1(18e20));
742
743  for(int bin = 0; bin < nbins; ++bin){
744  if(i == 0){
745  curves[bin] = new TGraph;
746  points[bin] = new TGraph;
747  }
748
749  const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
750
751  if(!std::isnan(ratio)) curves[bin]->SetPoint(curves[bin]->GetN(), x, ratio);
752  else curves[bin]->SetPoint(curves[bin]->GetN(), x, 1);
753  } // end for bin
754  } // end for i (x)
755
756  // As elswhere, to allow BirksC etc that need a different nominal to plot
757  // right.
758  IPrediction* pNom = 0;
759  for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
760  if(it->second.shifts[shiftIdx] == 0) pNom = it->second.preds[shiftIdx];;
761  }
762  if(pNom){ // if not, probably MinimizeMemory() was called
763  std::unique_ptr<TH1> hnom(pNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
764
765  for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
766  if(!it->second.preds[shiftIdx]) continue; // Probably MinimizeMemory()
767  std::unique_ptr<TH1> h(it->second.preds[shiftIdx]->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
768
769  for(int bin = 0; bin < nbins; ++bin){
770  const double ratio = h->GetBinContent(bin+1)/hnom->GetBinContent(bin+1);
771  if(!std::isnan(ratio)) points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], ratio);
772  else points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], 1);
773  }
774  } // end for shiftIdx
775  } // end if pNom
776
777
778  int nx = int(sqrt(nbins));
779  int ny = int(sqrt(nbins));
780  if(nx*ny < nbins) ++nx;
781  if(nx*ny < nbins) ++ny;
782
783  TCanvas* c = new TCanvas;
784  c->Divide(nx, ny);
785
786  for(int bin = 0; bin < nbins; ++bin){
787  c->cd(bin+1);
788  (new TH2F("",
789  TString::Format("%s %g < %s < %g;Shift;Ratio",
790  it->second.systName.c_str(),
791  nom->GetXaxis()->GetBinLowEdge(bin+1),
792  nom->GetXaxis()->GetTitle(),
793  nom->GetXaxis()->GetBinUpEdge(bin+1)),
794  100, -4, +4, 100, .5, 1.5))->Draw();
795  curves[bin]->Draw("l same");
796  points[bin]->SetMarkerStyle(kFullDotMedium);
797  points[bin]->Draw("p same");
798  } // end for bin
799
800  c->cd(0);
801  }
802
803  //----------------------------------------------------------------------
805  const std::string& savePattern,
806  Flavors::Flavors_t flav,
808  Sign::Sign_t sign) const
809  {
810  for(auto& it: fPreds){
811  DebugPlot(it.first, calc, flav, curr, sign);
812
813  if(!savePattern.empty()){
814  assert(savePattern.find("%s") != std::string::npos);
816  }
817  } // end for it
818  }
819
820  //----------------------------------------------------------------------
822  osc::IOscCalc* calc,
823  Flavors::Flavors_t flav,
825  Sign::Sign_t sign) const
826  {
827  InitFits();
828
829  std::unique_ptr<TH1> nom(fPredNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
830  const int nbins = nom->GetNbinsX();
831
832  TH2* h2 = new TH2F("", (syst->LatexName()+";;#sigma").c_str(),
833  nbins, nom->GetXaxis()->GetXmin(), nom->GetXaxis()->GetXmax(),
834  80, -4, +4);
835  h2->GetXaxis()->SetTitle(nom->GetXaxis()->GetTitle());
836
837  for(int i = 1; i <= 80; ++i){
838  const double y = h2->GetYaxis()->GetBinCenter(i);
839  const SystShifts ss(syst, y);
840  std::unique_ptr<TH1> h(PredictComponentSyst(calc, ss, flav, curr, sign).ToTH1(18e20));
841
842  for(int bin = 0; bin < nbins; ++bin){
843  const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
844
845  if(!isnan(ratio) && !isinf(ratio))
846  h2->Fill(h2->GetXaxis()->GetBinCenter(bin), y, ratio);
847  } // end for bin
848  } // end for i (x)
849
850  h2->Draw("colz");
851  h2->SetMinimum(0.5);
852  h2->SetMaximum(1.5);
853  }
854
855  //----------------------------------------------------------------------
857  const std::string& savePattern,
858  Flavors::Flavors_t flav,
860  Sign::Sign_t sign) const
861  {
862  InitFits();
863
864  for(auto& it: fPreds){
865  new TCanvas;
866  DebugPlotColz(it.first, calc, flav, curr, sign);
867
868  if(!savePattern.empty()){
869  assert(savePattern.find("%s") != std::string::npos);
871  }
872  } // end for it
873  }
874
875 } // 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:264
const XML_Char * name
Definition: expat.h:151
tree Draw("slc.nhit")
bool HasAnyStan() const
Definition: SystShifts.h:64
Implements systematic errors by interpolation between shifted templates.
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
double corr
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
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:149
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.
T sqrt(T number)
Definition: d0nt_math.hpp:156
General interface to oscillation calculators.
Definition: StanTypedefs.h:23
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:189
Float_t ss
Definition: plot.C:24
std::vector< std::vector< std::vector< Coeffs > > > fitsNubar
Will be filled if signs are separated, otherwise not.
std::vector< double > Spectrum
Definition: Constants.h:743
void SetOscSeed(osc::IOscCalc *oscSeed)
OStream cerr
Definition: OStream.cxx:7
void Clear()
Definition: Spectrum.cxx:362
Float_t tmp
Definition: plot.C:36
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
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:20
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
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::void_t< T > n
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
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:18
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:18
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:227
Oscillation probability calculators.
Definition: Calcs.h:5
virtual TMD5 * GetParamsHash() const
Use to check two calculators are in the same state.
Definition: IOscCalc.h:39
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
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: StanUtils.h:7
TDirectory * dir
Definition: macro.C:5
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:221
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:263
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:187
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
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:81
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:190
enum BeamMode string