PredictionInterp.cxx
Go to the documentation of this file.
5 #include "CAFAna/Core/Ratio.h"
6 #include "CAFAna/Core/Registry.h"
8
9 #include "TDirectory.h"
10
11 #include "TH1.h"
12 #include "TH1D.h"
13 #include "TH2.h"
14 #include "TMatrixD.h"
15 #include "TObjString.h"
16 #include "TVectorD.h"
17
18 // For debug plots
19 #include "TGraph.h"
20 #include "TCanvas.h"
21
23 #include "Utilities/func/MathUtil.h"
24 #include "Utilities/func/Stan.h"
26
28
29 #ifndef DARWINBUILD
30 #include <malloc.h>
31 #endif
32
33 namespace ana
34 {
35  //----------------------------------------------------------------------
37  {
38  // Copy fits into the remapped indexing order. TODO this is very ugly. Best
39  // would be to generate things in this order natively.
40
41  for(bool nubar: {false, true}){
42  const auto& from = nubar ? fitsNubar : fits;
43  auto& to = nubar ? fitsNubarRemap : fitsRemap;
44
45  // Allocate the space
46  to.resize(from.size());
47  for(auto& it: to){
48  it.resize(from[0][0].size());
49  for(auto& it2: it) it2.resize(from[0].size(), Coeffs(0, 0, 0, 0));
50  }
51
52  // And copy with the transposed indices
53  for(unsigned int i = 0; i < to.size(); ++i)
54  for(unsigned int j = 0; j < to[i].size(); ++j)
55  for(unsigned int k = 0; k < to[i][j].size(); ++k)
56  to[i][j][k] = from[i][k][j];
57
58  } // end for nubar
59  }
60
61  //----------------------------------------------------------------------
62  PredictionInterp::PredictionInterp(std::vector<const ISyst*> systs,
64  const IPredictionGenerator& predGen,
66  const SystShifts& shiftMC,
67  EMode_t mode)
68  : fOscOrigin(osc ? osc->Copy() : 0),
69  fBinning(0, {}, {}, 0, 0),
71  {
72  for(const ISyst* syst: systs){
73  ShiftedPreds sp;
74  sp.systName = syst->ShortName();
75
76  // Genie reweight shifts aren't able to provide +/-3 sigma
77  if(syst->IsGenieReweight())
78  sp.shifts = {-2, -1, 0, +1, +2};
79  else
80  sp.shifts = {-3, -2, -1, 0, +1, +2, +3};
81
82  for(int sigma: sp.shifts){
83  SystShifts shiftHere = shiftMC;
84  shiftHere.SetShift(syst, sigma);
86  }
87
88  fPreds.emplace(syst, sp);
89  } // end for syst
90
92  }
93
94  //----------------------------------------------------------------------
96  {
97  // for(auto it: fPreds) for(IPrediction* p: it.second.preds) delete p;
98  // delete fOscOrigin;
99
100  // It isn't really a unique ptr when we use PredictionInterpTemplates
101  fPredNom.release();
102  }
103
104  //----------------------------------------------------------------------
105  std::vector<std::vector<PredictionInterp::Coeffs>>
106  PredictionInterp::FitRatios(const std::vector<double>& shifts,
107  const std::vector<std::unique_ptr<TH1>>& ratios) const
108  {
109  assert(shifts.size() == ratios.size());
110  std::vector<std::vector<Coeffs>> ret;
111
112  const int binMax = ratios[0]->GetNbinsX();
113
114  for(int binIdx = 0; binIdx < binMax+2; ++binIdx){
115  ret.push_back({});
116  // This is cubic interpolation. For each adjacent set of four points we
117  // determine coefficients for a cubic which will be the curve between the
118  // center two. We constrain the function to match the two center points
119  // and to have the right mean gradient at them. This causes this patch to
120  // match smoothly with the next one along. The resulting function is
121  // continuous and first and second differentiable. At the ends of the
122  // range we fit a quadratic instead with only one constraint on the
123  // slope. The coordinate conventions are that point y1 sits at x=0 and y2
124  // at x=1. The matrices are simply the inverses of writing out the
125  // constraints expressed above.
126
127  // Special-case for linear interpolation
128  if(ratios.size() == 2){
129  const double y0 = ratios[0]->GetBinContent(binIdx);
130  const double y1 = ratios[1]->GetBinContent(binIdx);
131
132  ret.back().emplace_back(0, 0, y1-y0, y0);
133  continue;
134  }
135
136  {
137  const double y1 = ratios[0]->GetBinContent(binIdx);
138  const double y2 = ratios[1]->GetBinContent(binIdx);
139  const double y3 = ratios[2]->GetBinContent(binIdx);
140  const double v[3] = {y1, y2, (y3-y1)/2};
141  const double m[9] = { 1, -1, 1,
142  -2, 2, -1,
143  1, 0, 0};
144  const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
145  ret.back().emplace_back(0, res(0), res(1), res(2));
146  }
147
148  // We're assuming here that the shifts are separated by exactly 1 sigma.
149  for(unsigned int shiftIdx = 1; shiftIdx < ratios.size()-2; ++shiftIdx){
150  const double y0 = ratios[shiftIdx-1]->GetBinContent(binIdx);
151  const double y1 = ratios[shiftIdx ]->GetBinContent(binIdx);
152  const double y2 = ratios[shiftIdx+1]->GetBinContent(binIdx);
153  const double y3 = ratios[shiftIdx+2]->GetBinContent(binIdx);
154
155  const double v[4] = {y1, y2, (y2-y0)/2, (y3-y1)/2};
156  const double m[16] = { 2, -2, 1, 1,
157  -3, 3, -2, -1,
158  0, 0, 1, 0,
159  1, 0, 0, 0};
160  const TVectorD res = TMatrixD(4, 4, m) * TVectorD(4, v);
161  ret.back().emplace_back(res(0), res(1), res(2), res(3));
162  } // end for shiftIdx
163
164  {
165  const int N = ratios.size()-3;
166  const double y0 = ratios[N ]->GetBinContent(binIdx);
167  const double y1 = ratios[N+1]->GetBinContent(binIdx);
168  const double y2 = ratios[N+2]->GetBinContent(binIdx);
169  const double v[3] = {y1, y2, (y2-y0)/2};
170  const double m[9] = {-1, 1, -1,
171  0, 0, 1,
172  1, 0, 0};
173  const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
174  ret.back().emplace_back(0, res(0), res(1), res(2));
175  }
176  } // end for binIdx
177
178  double stride = -1;
179  for(unsigned int i = 0; i < shifts.size()-1; ++i){
180  const double newStride = shifts[i+1]-shifts[i];
181  assert((stride < 0 || fabs(stride-newStride) < 1e-3) &&
182  "Variably-spaced syst templates are unsupported");
183  stride = newStride;
184  }
185
186  // If the stride is actually not 1, need to rescale all the coefficients
187  for(std::vector<Coeffs>& cs: ret)
188  for(Coeffs& c: cs){
189  c = Coeffs(c.a/util::cube(stride),
190  c.b/util::sqr(stride),
191  c.c/stride,
192  c.d);}
193  return ret;
194  }
195
196  //----------------------------------------------------------------------
197  std::vector<std::vector<PredictionInterp::Coeffs>>
198  PredictionInterp::FitComponent(const std::vector<double>& shifts,
199  const std::vector<IPrediction*>& preds,
200  Flavors::Flavors_t flav,
202  Sign::Sign_t sign) 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<std::unique_ptr<TH1>> ratios;
217  for(auto& p: preds){
218  ratios.emplace_back(Ratio(p->PredictComponent(fOscOrigin,
219  flav, curr, sign),
220  nom).ToTH1());
221
222  // Check none of the ratio values is utterly crazy
223  std::unique_ptr<TH1>& r = ratios.back();
224  for(int i = 0; i < r->GetNbinsX()+2; ++i){
225  const double y = r->GetBinContent(i);
226  if(y > 500){
227  std::cout << "PredictionInterp: WARNING, ratio in bin "
228  << i << " is " << y << ". Ignoring." << std::endl;
229  r->SetBinContent(i, 1);
230  }
231  }
232  }
233
234  return FitRatios(shifts, ratios);
235  }
236
237  //----------------------------------------------------------------------
238  std::vector<std::vector<PredictionInterp::Coeffs>>
239  PredictionInterp::FitSumComponent(const std::vector<ShiftedPreds>& sps,
240  Flavors::Flavors_t flav,
242  Sign::Sign_t sign) const
243  {
244  IPrediction* pNom = 0;
245  //take 1st one as template
246  std::vector<double> shifts = sps[0].shifts;
247  std::vector<IPrediction*> preds = sps[0].preds;
248
249  //for now, genie only take +-2
250  assert(shifts[2] == 0 && "PredictionSumInterp only implemented for +-2 shifts");
251
252  for(unsigned int i = 0; i < shifts.size(); ++i){
253  if(shifts[i] == 0) pNom = preds[i];
254  }
255  assert(pNom);
256
257  const Spectrum nom = pNom->PredictComponent(fOscOrigin,
258  flav, curr, sign);
259
260  std::vector<std::unique_ptr<TH1>> sumratios;
261
262  double pot = nom.POT();
263  const Binning binning = nom.GetBinnings()[0];
265  TH1D* hnom = nom.ToTH1(pot);
266  TH1D* hsummed = HistCache::Copy(hnom);
267
269  int nbin = hnom->GetNbinsX();
270  double* from = hnom->GetArray();
271
272  for(unsigned int i = 0; i < shifts.size(); i++){
273  for(auto& sp: sps){
274
275  IPrediction* p = sp.preds[i];
276  // force-convert to _Spectrum<double> since we don't need autodiff here
277  TH1D* hpred = Spectrum(p->PredictComponent(fOscOrigin, flav, curr, sign)).ToTH1(pot);
278
279  double* to = hpred->GetArray();
280  for(int bin = 0; bin < nbin+2; ++bin){
281  to[bin] -= from[bin];
282  to[bin] *= to[bin];
284  }//end all bins
285  HistCache::Delete(hpred);
286  }//end all systpreds
287
288  //sqrt of sum of diff
289  for(int bin = 0; bin < nbin+2; ++bin){
294  }
295
296  std::unique_ptr<TH1> hret(HistCache::New("", binning));
297
298  for( int i = 0; i< nbin+2; ++i){ hret->SetBinContent( i, quads[i]);}
299
300  sumratios.emplace_back(std::move(hret));
301  } //end shifts
302  HistCache::Delete(hsummed);
303  HistCache::Delete(hnom);
304  return FitRatios(shifts, sumratios);
305  }
306
307  //----------------------------------------------------------------------
309  std::vector<std::vector<std::vector<Coeffs>>>& fits,
310  Sign::Sign_t sign) const
311  {
312  fits.resize(kNCoeffTypes);
313
317
319
321  }
322
323  //----------------------------------------------------------------------
324  void PredictionInterp::InitFitsSumHelper(std::vector<ShiftedPreds>& sp,
325  std::vector<std::vector<std::vector<Coeffs>>>& fits,
326  Sign::Sign_t sign) const
327  {
328  fits.resize(kNCoeffTypes);
332
334
336  }
337  //----------------------------------------------------------------------
339  {
340
341  // No systs
342  if(fPreds.empty() && fSumPreds.empty()){
343  if(fBinning.POT() > 0 || fBinning.Livetime() > 0) return;
344  }
345  // Already initialized, Summed merged to fPreds
346  else if(!fPreds.empty() && !fPreds.begin()->second.fits.empty()) return;
347
348  for(auto& it: fPreds){
349  ShiftedPreds& sp = it.second;
350
351  if(fSplitBySign){
352  InitFitsHelper(sp, sp.fits, Sign::kNu);
354  }
355  else{
357  }
358  sp.nCoeffs = sp.fits[0][0].size();
359
360  // Copy the outputs into the remapped indexing order. TODO this is very
361  // ugly. Best would be to generate things in this order natively.
362  sp.FillRemaps();
363  }
364
365  if(!fSumPreds.empty()){
366  std::cout<<"PredictionInterp: Summing "<<fSumPreds.size()<<
367  " systematics..."<<std::endl;
368  std::vector<ShiftedPreds> sp;
369  for(auto& it: fSumPreds) sp.push_back( it.second);
370
371  fSummedsp.systName = "Summed Syst";
372  fSummedsp.shifts = {-2, -1, 0, 1, 2};
373
374  if(fSplitBySign){
377  }
378  else{
380  }
381  fSummedsp.nCoeffs = fSummedsp.fits[0][0].size();
382
383  fPreds.emplace(&kSummedSyst, fSummedsp);
384  }
385
386  // Predict something, anything, so that we can know what binning to use
387  fBinning = fPredNom->Predict(fOscOrigin);
388  fBinning.Clear();
389  }
390
391  //----------------------------------------------------------------------
393  fOscOrigin = oscSeed->Copy();
394  for(auto& it: fPreds) it.second.fits.clear();
395  InitFits();
396  }
397
398  //----------------------------------------------------------------------
400  {
401  return fPredNom->Predict(calc);
402  }
403
404  //----------------------------------------------------------------------
406  {
407  return fPredNom->Predict(calc);
408  }
409
410  //----------------------------------------------------------------------
412  Flavors::Flavors_t flav,
414  Sign::Sign_t sign) const
415  {
416  return fPredNom->PredictComponent(calc, flav, curr, sign);
417  }
418
419  //----------------------------------------------------------------------
421  Flavors::Flavors_t flav,
423  Sign::Sign_t sign) const
424  {
425  return fPredNom->PredictComponent(calc, flav, curr, sign);
426  }
427
428  //----------------------------------------------------------------------
430  const SystShifts& shift) const
431  {
432  InitFits();
433
434  return PredictComponentSyst(calc,
435  shift,
438  Sign::kBoth);
439  }
440
441  //----------------------------------------------------------------------
443  const SystShifts& shift) const
444  {
445  InitFits();
446
447  return PredictComponentSyst(calc,
448  shift,
451  Sign::kBoth);
452  }
453
454  //----------------------------------------------------------------------
457  bool nubar,
458  const SystShifts& shift) const
459  {
460  // TODO histogram operations could be too slow
461  TH1D *h = s.ToTH1(s.POT());
462
463  ShiftBins(h->GetNbinsX()+2, h->GetArray(), type, nubar, shift);
464
465  return Spectrum(std::unique_ptr<TH1D>(h), s.GetLabels(), s.GetBinnings(), s.POT(), s.Livetime());
466  }
467
468  //----------------------------------------------------------------------
471  bool nubar,
472  const SystShifts & shift) const
473  {
474
475  auto binContents = s.ToBins(s.POT());
476
477  ShiftBins(binContents.size(), &binContents.front(), type, nubar, shift);
478
479  return SpectrumStan(std::move(binContents), s.GetLabels(), s.GetBinnings(), s.POT(), s.Livetime());
480  }
481
482  //----------------------------------------------------------------------
483  template <typename T>
484  void PredictionInterp::ShiftBins(unsigned int N,
485  T* arr,
487  bool nubar,
488  const SystShifts& shift) const
489  {
491  "PredictionInterp::ShiftBins() can only be called using doubles or stan::math::vars");
492  if(nubar) assert(fSplitBySign);
493
494  T corr[N];
495  for(unsigned int i = 0; i < N; ++i) corr[i] = 1;
496
497  for(auto& it: fPreds){
498  const ISyst* syst = it.first;
499  const ShiftedPreds& sp = it.second;
500
501  T x = shift.GetShift<T>(syst);
502
503  // need to actually do the calculation for the autodiff version
504  // to make sure the gradient is computed correctly
505  if(x == 0 && !std::is_same<T, stan::math::var>::value) continue;
506
507  int shiftBin = (util::GetValAs<double>(x) - sp.shifts[0])/sp.Stride();
508  shiftBin = std::max(0, shiftBin);
509  shiftBin = std::min(shiftBin, sp.nCoeffs-1);
510
511  const Coeffs* fits = nubar ? &sp.fitsNubarRemap[type][shiftBin].front() : &sp.fitsRemap[type][shiftBin].front();
512
513  x -= sp.shifts[shiftBin];
514
515  const T x_cube = util::cube(x);
516  const T x_sqr = util::sqr(x);
517
518  for(unsigned int n = 0; n < N; ++n){
519  // Uncomment to debug crashes in this function
520  // assert(type < fits.size());
521  // assert(n < sp.fits[type].size());
522  // assert(shiftBin < int(sp.fits[type][n].size()));
523
524  const Coeffs& f = fits[n];
525
526  corr[n] *= f.a*x_cube + f.b*x_sqr + f.c*x + f.d;
527  } // end for n
528  } // end for syst
529
530  for(unsigned int n = 0; n < N; ++n){
531  // std::max() doesn't work with stan::math::var
532  arr[n] *= (corr[n] > 0.) ? corr[n] : 0.;
533  }
534
535  }
536
537  //----------------------------------------------------------------------
539  const TMD5* hash,
540  const SystShifts& shift,
541  Flavors::Flavors_t flav,
544  CoeffsType type) const
545  {
546  return _ShiftedComponent<Spectrum>(calc, hash, shift, flav, curr, sign, type);
547  }
548
549  //----------------------------------------------------------------------
551  const TMD5* hash,
552  const SystShifts& shift,
553  Flavors::Flavors_t flav,
556  CoeffsType type) const
557  {
558  return _ShiftedComponent<SpectrumStan>(calc, hash, shift, flav, curr, sign, type);
559  }
560
561  //----------------------------------------------------------------------
562  template <typename U, typename T>
564  const TMD5* hash,
565  const SystShifts& shift,
566  Flavors::Flavors_t flav,
569  CoeffsType type) const
570  {
572  "PredictionInterp::ShiftedComponent() can only be called using doubles or stan::math::vars");
573
574  if(fSplitBySign && sign == Sign::kBoth){
575  return (ShiftedComponent(calc, hash, shift, flav, curr, Sign::kAntiNu, type)+
576  ShiftedComponent(calc, hash, shift, flav, curr, Sign::kNu, type));
577  }
578
579  // Must be the base case of the recursion to use the cache. Otherwise we
580  // can cache systematically shifted versions of our children, which is
581  // wrong. Also, some calculators won't hash themselves.
582  // Moreover, caching is not going to work with stan::math::vars
583  // since they get reset every time Stan's log_prob() is called.
584  const bool canCache = (hash != 0) && !std::is_same<T, stan::math::var>::value;
585
586  const Key_t key = {flav, curr, sign};
587  auto it = fNomCache.find(key);
588
589  // Should the interpolation use the nubar fits?
590  const bool nubar = (fSplitBySign && sign == Sign::kAntiNu);
591
592  // We have the nominal for this exact combination of flav, curr, sign, calc
593  // stored. Shift it and return.
594  if(canCache && it != fNomCache.end() && it->second.hash == *hash){
595  return ShiftSpectrum(it->second.nom, type, nubar, shift);
596  }
597
598  // We need to compute the nominal again for whatever reason
599  const auto nom = fPredNom->PredictComponent(calc, flav, curr, sign);
600
601  if(canCache){
602  // Insert into the cache if not already there, or update if there but
603  // with old oscillation parameters.
604  if(it == fNomCache.end())
605  fNomCache.emplace(key, Val_t({*hash, nom}));
606  else
607  it->second = {*hash, nom};
608  }
609
610  return ShiftSpectrum(nom, type, nubar, shift);
611  }
612
613  //----------------------------------------------------------------------
614  template <typename U, typename T>
616  const SystShifts& shift,
617  Flavors::Flavors_t flav,
619  Sign::Sign_t sign) const
620  {
621  InitFits();
622
623  U ret = fBinning;
624  ret.Clear();
625
626  // Check that we're able to handle all the systs we were passed
627  for(const ISyst* syst: shift.ActiveSysts()){
628  if(fPreds.find(syst) == fPreds.end()){
629  std::cerr << "This PredictionInterp is not set up to handle the requested systematic: " << syst->ShortName() << std::endl;
630  abort();
631  }
632  } // end for syst
633
634
635  const TMD5* hash = calc ? calc->GetParamsHash() : 0;
636
637  if(curr & Current::kCC){
638  if(flav & Flavors::kNuEToNuE) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuE, Current::kCC, sign, kNueSurv);
639  if(flav & Flavors::kNuEToNuMu) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuMu, Current::kCC, sign, kOther );
640  if(flav & Flavors::kNuEToNuTau) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuTau, Current::kCC, sign, kOther );
641
642  if(flav & Flavors::kNuMuToNuE) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuE, Current::kCC, sign, kNueApp );
643  if(flav & Flavors::kNuMuToNuMu) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuMu, Current::kCC, sign, kNumuSurv);
644  if(flav & Flavors::kNuMuToNuTau) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuTau, Current::kCC, sign, kOther );
645  }
646  if(curr & Current::kNC){
647  assert(flav == Flavors::kAll); // Don't know how to calculate anything else
648
649  ret += ShiftedComponent(calc, hash, shift, Flavors::kAll, Current::kNC, sign, kNC);
650  }
651
652  delete hash;
653
654  return ret;
655  }
656
657  //----------------------------------------------------------------------
658  // can't template these directly since the interface isn't templated
660  const SystShifts& shift,
661  Flavors::Flavors_t flav,
663  Sign::Sign_t sign) const
664  {
665  return _PredictComponentSyst<Spectrum>(calc, shift, flav, curr, sign);
666  }
667
668  //----------------------------------------------------------------------
670  const SystShifts& shift,
671  Flavors::Flavors_t flav,
673  Sign::Sign_t sign) const
674  {
675  return _PredictComponentSyst<SpectrumStan>(calc, shift, flav, curr, sign);
676  }
677
678  //----------------------------------------------------------------------
680  Flavors::Flavors_t flav,
684  const SystShifts& shift,
685  double pot,
686  std::unordered_map<const ISyst*, std::vector<double>>& dp) const
687  {
688  if(fSplitBySign && sign == Sign::kBoth){
689  ComponentDerivative(calc, flav, curr, Sign::kNu, type, shift, pot, dp);
690  ComponentDerivative(calc, flav, curr, Sign::kAntiNu, type, shift, pot, dp);
691  return;
692  }
693
694  // this method is only ever going to be used by the "Frequentist" version of the class.
695  // don't bother templating.
696  const Spectrum base = PredictComponentSyst(calc, shift, flav, curr, sign);
697  TH1D* h = base.ToTH1(pot);
698  double* arr = h->GetArray();
699  const unsigned int N = h->GetNbinsX()+2;
700
701
702  // Should the interpolation use the nubar fits?
703  const bool nubar = (fSplitBySign && sign == Sign::kAntiNu);
704
705  for(auto& it: dp){
706  const ISyst* syst = it.first;
707  std::vector<double>& diff = it.second;
708  if(diff.empty()) diff.resize(N);
709  assert(diff.size() == N);
710  assert(fPreds.find(syst) != fPreds.end());
711  const ShiftedPreds& sp = fPreds[syst];
712
713  // not instrumenting this method for stan::math::vars
714  // since we're using autodiff anyway
715  double x = util::GetValAs<double>(shift.GetShift(syst));
716
717  int shiftBin = (x - sp.shifts[0])/sp.Stride();
718  shiftBin = std::max(0, shiftBin);
719  shiftBin = std::min(shiftBin, sp.nCoeffs-1);
720
721  const Coeffs* fits = nubar ? &sp.fitsNubarRemap[type][shiftBin].front() : &sp.fitsRemap[type][shiftBin].front();
722
723  x -= sp.shifts[shiftBin];
724
725  const double x_cube = util::cube(x);
726  const double x_sqr = util::sqr(x);
727
728  for(unsigned int n = 0; n < N; ++n){
729  // Uncomment to debug crashes in this function
730  // assert(type < fits.size());
731  // assert(n < sp.fits[type].size());
732  // assert(shiftBin < int(sp.fits[type][n].size()));
733  const Coeffs& f = fits[n];
734
735  const double corr = f.a*x_cube + f.b*x_sqr + f.c*x + f.d;
736  if(corr > 0) diff[n] += (3*f.a*x_sqr + 2*f.b*x + f.c)/corr*arr[n];
737  } // end for n
738  } // end for syst
739
741  }
742
743  //----------------------------------------------------------------------
745  const SystShifts& shift,
746  double pot,
747  std::unordered_map<const ISyst*, std::vector<double>>& dp) const
748  {
749  InitFits();
750
751  // Check that we're able to handle all the systs we were passed
752  for(auto& it: dp){
753  if(fPreds.find(it.first) == fPreds.end()){
754  std::cerr << "This PredictionInterp is not set up to handle the requested systematic: " << it.first->ShortName() << std::endl;
755  abort();
756  }
757  it.second.clear();
758  } // end for syst
759
763
767
769
770  // Simpler (much slower) implementation in terms of finite differences for
771  // test purposes
772  /*
773  const Spectrum p0 = PredictSyst(calc, shift);
774  TH1D* h0 = p0.ToTH1(pot);
775
776  const double dx = 1e-9;
777  for(auto& it: dp){
778  const ISyst* s = it.first;
779  std::vector<double>& v = it.second;
780  SystShifts s2 = shift;
781  s2.SetShift(s, s2.GetShift(s)+dx);
782
783  const Spectrum p1 = PredictSyst(calc, s2);
784
785  TH1D* h1 = p1.ToTH1(pot);
786
787  v.resize(h1->GetNbinsX()+2);
788  for(int i = 0; i < h1->GetNbinsX()+2; ++i){
789  v[i] = (h1->GetBinContent(i) - h0->GetBinContent(i))/dx;
790  }
791
792  HistCache::Delete(h1);
793  }
794
795  HistCache::Delete(h0);
796  */
797  }
798
799  //----------------------------------------------------------------------
800  void PredictionInterp::SaveTo(TDirectory* dir) const
801  {
802  InitFits();
803
804  TDirectory* tmp = gDirectory;
805
806  dir->cd();
807  TObjString("PredictionInterp").Write("type");
808
809
810  fPredNom->SaveTo(dir->mkdir("pred_nom"));
811
812
813  for(auto it: fPreds){
814  const ShiftedPreds& sp = it.second;
815
816  for(unsigned int i = 0; i < sp.shifts.size(); ++i){
817  if(!sp.preds[i]){
818  std::cout << "Can't save a PredictionInterp after MinimizeMemory()" << std::endl;
819  abort();
820  }
821  sp.preds[i]->SaveTo(dir->mkdir(TString::Format("pred_%s_%+d",
822  sp.systName.c_str(),
823  int(sp.shifts[i])).Data()));
824  } // end for i
825  } // end for it
826
827  ana::SaveTo(*fOscOrigin, dir->mkdir("osc_origin"));
828
829  if(!fPreds.empty()){
830  TH1F hSystNames("syst_names", ";Syst names", fPreds.size(), 0, fPreds.size());
831  int binIdx = 1;
832  for(auto it: fPreds){
833  hSystNames.GetXaxis()->SetBinLabel(binIdx++, it.second.systName.c_str());
834  }
835  hSystNames.Write("syst_names");
836  }
837
838  TObjString split_sign = fSplitBySign ? "yes" : "no";
839  split_sign.Write("split_sign");
840
841  tmp->cd();
842  }
843
844  //----------------------------------------------------------------------
846  {
847  TObjString* tag = (TObjString*)dir->Get("type");
848  assert(tag);
849  assert(tag->GetString() == "PredictionInterp" ||
850  tag->GetString() == "PredictionInterp2"); // backwards compatibility
851
852  auto ret = std::make_unique<PredictionInterp>();
853
855
856  TObjString* split_sign = (TObjString*)dir->Get("split_sign");
857  // Can be missing from old files
858  ret->fSplitBySign = (split_sign && split_sign->String() == "yes");
859
860  return ret;
861  }
862
863  //----------------------------------------------------------------------
865  std::vector<const ISyst*> veto)
866  {
868
869  TH1* hSystNames = (TH1*)dir->Get("syst_names");
870  if(hSystNames){
871  for(int systIdx = 0; systIdx < hSystNames->GetNbinsX(); ++systIdx){
872  ShiftedPreds sp;
873  sp.systName = hSystNames->GetXaxis()->GetBinLabel(systIdx+1);
874
876
877  if(std::find(veto.begin(), veto.end(), syst) != veto.end()) continue;
878
879  for(int shift = -3; shift <= +3; ++shift){
880  TDirectory* preddir = dir->GetDirectory(TString::Format("pred_%s_%+d", sp.systName.c_str(), shift).Data());
881  if(!preddir) continue; // Can happen for genie systs
882
884
885  sp.shifts.push_back(shift);
886  sp.preds.push_back(pred);
887  } // end for shift
888
889  ret->fPreds.emplace(syst, sp);
890  } // end for systIdx
891  } // end if hSystNames
892
894  }
895
896  //----------------------------------------------------------------------
898  {
899  std::set<IPrediction*> todel;
900  for(auto& it: fPreds){
901  std::vector<IPrediction*>& preds = it.second.preds;
902  for(unsigned int i = 0; i < preds.size(); ++i){
903  if(preds[i] != fPredNom.get()){
904  todel.insert(preds[i]);
905  preds[i] = 0;
906  }
907  }
908  }
909
910  for(IPrediction* p: todel) delete p;
911
912  // We probably just freed up a lot of memory, but malloc by default hangs
913  // on to all of it as cache.
914  #ifndef DARWINBUILD
915  malloc_trim(0);
916  #endif
917  }
918
919  //----------------------------------------------------------------------
922  Flavors::Flavors_t flav,
924  Sign::Sign_t sign) const
925  {
926  InitFits();
927
928  auto it = fPreds.find(syst);
929  if(it == fPreds.end()){
930  std::cout << "PredictionInterp::DebugPlots(): "
932  return;
933  }
934
935
936  // force-convert to _Spectrum<double> for ease here since in DebugPlots()
937  // we're not concerned with preserving autodiff or speed.
938  // (more instances further below.)
939  std::unique_ptr<TH1> nom(Spectrum(fPredNom->PredictComponent(calc, flav, curr, sign)).ToTH1(18e20));
940  const int nbins = nom->GetNbinsX();
941
942  TGraph* curves[nbins];
943  TGraph* points[nbins];
944
945  for(int i = 0; i <= 80; ++i){
946  const double x = .1*i-4;
947  const SystShifts ss(it->first, x);
948  std::unique_ptr<TH1> h(Spectrum(PredictComponentSyst(calc, ss, flav, curr, sign)).ToTH1(18e20));
949
950  for(int bin = 0; bin < nbins; ++bin){
951  if(i == 0){
952  curves[bin] = new TGraph;
953  points[bin] = new TGraph;
954  }
955
956  const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
957
958  if(!std::isnan(ratio)) curves[bin]->SetPoint(curves[bin]->GetN(), x, ratio);
959  else curves[bin]->SetPoint(curves[bin]->GetN(), x, 1);
960  } // end for bin
961  } // end for i (x)
962
963  // As elswhere, to allow BirksC etc that need a different nominal to plot
964  // right.
965  IPrediction* pNom = 0;
966  for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
967  if(it->second.shifts[shiftIdx] == 0) pNom = it->second.preds[shiftIdx];
968  }
969  if(pNom){ // if not, probably MinimizeMemory() was called
970  std::unique_ptr<TH1> hnom(Spectrum(pNom->PredictComponent(calc, flav, curr, sign)).ToTH1(18e20));
971
972  for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
973  if(!it->second.preds[shiftIdx]) continue; // Probably MinimizeMemory()
974  std::unique_ptr<TH1> h;
975  h = std::move(std::unique_ptr<TH1>(Spectrum(it->second.preds[shiftIdx]->PredictComponent(calc, flav, curr, sign)).ToTH1(18e20)));
976
977  for(int bin = 0; bin < nbins; ++bin){
978  const double ratio = h->GetBinContent(bin+1)/hnom->GetBinContent(bin+1);
979  if(!std::isnan(ratio)) points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], ratio);
980  else points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], 1);
981  }
982  } // end for shiftIdx
983  } // end if pNom
984
985
986  int nx = int(sqrt(nbins));
987  int ny = int(sqrt(nbins));
988  if(nx*ny < nbins) ++nx;
989  if(nx*ny < nbins) ++ny;
990
991  TCanvas* c = new TCanvas;
992  c->Divide(nx, ny);
993
994  for(int bin = 0; bin < nbins; ++bin){
995  c->cd(bin+1);
996  (new TH2F("",
997  TString::Format("%s %g < %s < %g;Shift;Ratio",
998  it->second.systName.c_str(),
999  nom->GetXaxis()->GetBinLowEdge(bin+1),
1000  nom->GetXaxis()->GetTitle(),
1001  nom->GetXaxis()->GetBinUpEdge(bin+1)),
1002  100, -4, +4, 100, .5, 1.5))->Draw();
1003  curves[bin]->Draw("l same");
1004  points[bin]->SetMarkerStyle(kFullDotMedium);
1005  points[bin]->Draw("p same");
1006  } // end for bin
1007
1008  c->cd(0);
1009  }
1010
1011  //----------------------------------------------------------------------
1013  const std::string& savePattern,
1014  Flavors::Flavors_t flav,
1016  Sign::Sign_t sign) const
1017  {
1018  for(auto& it: fPreds){
1019  DebugPlot(it.first, calc, flav, curr, sign);
1020
1021  if(!savePattern.empty()){
1022  assert(savePattern.find("%s") != std::string::npos);
1024  }
1025  } // end for it
1026  }
1027
1028  //----------------------------------------------------------------------
1031  Flavors::Flavors_t flav,
1033  Sign::Sign_t sign) const
1034  {
1035  InitFits();
1036
1037  std::unique_ptr<TH1> nom(fPredNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
1038  const int nbins = nom->GetNbinsX();
1039
1040  TH2* h2 = new TH2F("", (syst->LatexName()+";;#sigma").c_str(),
1041  nbins, nom->GetXaxis()->GetXmin(), nom->GetXaxis()->GetXmax(),
1042  80, -4, +4);
1043  h2->GetXaxis()->SetTitle(nom->GetXaxis()->GetTitle());
1044
1045  for(int i = 1; i <= 80; ++i){
1046  const double y = h2->GetYaxis()->GetBinCenter(i);
1047  const SystShifts ss(syst, y);
1048  std::unique_ptr<TH1> h(PredictComponentSyst(calc, ss, flav, curr, sign).ToTH1(18e20));
1049
1050  for(int bin = 0; bin < nbins; ++bin){
1051  const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
1052
1053  if(!isnan(ratio) && !isinf(ratio))
1054  h2->Fill(h2->GetXaxis()->GetBinCenter(bin), y, ratio);
1055  } // end for bin
1056  } // end for i (x)
1057
1058  h2->Draw("colz");
1059  h2->SetMinimum(0.5);
1060  h2->SetMaximum(1.5);
1061  }
1062
1063  //----------------------------------------------------------------------
1065  const std::string& savePattern,
1066  Flavors::Flavors_t flav,
1068  Sign::Sign_t sign) const
1069  {
1070  InitFits();
1071
1072  for(auto it: fPreds){
1073  new TCanvas;
1074  DebugPlotColz(it.first, calc, flav, curr, sign);
1075
1076  if(!savePattern.empty()){
1077  assert(savePattern.find("%s") != std::string::npos);
1079  }
1080  } // end for it
1081  }
1082
1083 } // 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.
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
Implements systematic errors by interpolation between shifted templates.
osc::OscCalculatorDumb calc
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Oscillation analysis framework, runs over CAF files outside of ART.
void DebugPlot(const ISyst *syst, osc::IOscCalculator *calc, Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
set< int >::iterator it
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::string systName
What systematic we&#39;re interpolating.
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
Float_t y1[n_points_granero]
Definition: compare.C:5
static TH1D * Copy(const TH1D *h)
Definition: HistCache.cxx:175
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
std::vector< std::vector< std::vector< Coeffs > > > fitsNubarRemap
std::vector< stan::math::var > ToBins(double pot) const
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SetOscSeed(osc::IOscCalculator *oscSeed)
(&#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
std::vector< std::vector< Coeffs > > FitSumComponent(const std::vector< ShiftedPreds > &spf, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Float_t ss
Definition: plot.C:24
std::vector< std::vector< std::vector< Coeffs > > > fitsNubar
Will be filled if signs are separated, otherwise not.
virtual _IOscCalculator * Copy() const =0
OStream cerr
Definition: OStream.cxx:7
std::vector< std::vector< Coeffs > > FitRatios(const std::vector< double > &shifts, const std::vector< std::unique_ptr< TH1 >> &ratios) const
Find coefficients describing this set of shifts.
void Clear()
Definition: Spectrum.cxx:878
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 ComponentDerivative(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign, CoeffsType type, const SystShifts &shift, double pot, std::unordered_map< const ISyst *, std::vector< double >> &dp) const
Helper for Derivative.
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
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
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
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.
void SaveTo(const osc::IOscCalculator &x, TDirectory *dir)
std::map< Key_t, Val_t > fNomCache
U _ShiftedComponent(osc::_IOscCalculator< 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.
const std::vector< std::string > & GetLabels() const
Definition: SpectrumStan.h:44
const XML_Char * s
Definition: expat.h:262
const int nbins
Definition: cellShifts.C:15
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.
static void Delete(TH1D *&h)
Definition: HistCache.cxx:215
virtual void SaveTo(TDirectory *dir) const override
virtual TMD5 * GetParamsHash() const
Use to check two calculators are in the same state.
#define pot
T GetShift(const ISyst *syst) const
std::vector< IPrediction * > preds
osc::IOscCalculator * fOscOrigin
The oscillation values we assume when evaluating the coefficients.
const double j
Definition: BetheBloch.cxx:29
void DebugPlotColz(const ISyst *syst, osc::IOscCalculator *calc, Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
const DummyAnaSyst kSummedSyst("SummedSyst","Summed Syst.")
Spectrum fBinning
Dummy spectrum to provide binning.
double Livetime() const
Definition: SpectrumStan.h:46
std::unordered_map< const ISyst *, ShiftedPreds > fSumPreds
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:263
Spectrum ShiftedComponent(osc::IOscCalculator *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.
Oscillation probability calculators.
Definition: Calcs.h:5
Represent the ratio between two spectra.
Definition: Ratio.h:8
Neutrinos-only.
Definition: IPrediction.h:49
def sign(x)
Definition: canMan.py:204
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir)
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
Taus, numu appearance.
TH1F * h2
Definition: plot.C:45
void InitFitsSumHelper(std::vector< ShiftedPreds > &sps, std::vector< std::vector< std::vector< Coeffs >>> &fits, Sign::Sign_t sign) const
virtual const std::string & LatexName() const final
The name used on plots (ROOT&#39;s TLatex syntax)
Definition: ISyst.h:30
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
Find coefficients describing the ratios from this component.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
TVectorT< double > TVectorD
Definition: Utilities.h:18
void DebugPlots(osc::IOscCalculator *calc, const std::string &savePattern="", Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
std::vector< std::vector< std::vector< Coeffs > > > fitsRemap
Spectrum Predict(osc::IOscCalculator *calc) const override
TDirectory * dir
Definition: macro.C:5
string syst
Definition: plotSysts.py:176
Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
std::vector< std::string > fits
recTree Draw("energy.numu.trkccE>>precosmics","fdcuts&&preshutdown")
Neutral-current interactions.
Definition: IPrediction.h:40
TRandom3 r(0)
void Derivative(osc::IOscCalculator *calc, const SystShifts &shift, double pot, std::unordered_map< const ISyst *, std::vector< double >> &dp) const override
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
static TH1D * New(const std::string &title, const Binning &bins)
Definition: HistCache.cxx:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::vector< const ISyst * > ActiveSysts() const
Definition: SystShifts.cxx:199
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
void DebugPlotsColz(osc::IOscCalculator *calc, const std::string &savePattern="", Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
double T
Definition: Xdiff_gwt.C:5
Given loaders and an MC shift, Generate() generates an IPrediction.
std::vector< std::string > GetLabels() const
Definition: Spectrum.h:298
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
static void LoadFromBody(TDirectory *dir, PredictionInterp *ret, std::vector< const ISyst * > veto={})
T min(const caf::Proxy< T > &a, T b)
TMatrixT< double > TMatrixD
Definition: Utilities.h:16
Spectrum PredictSyst(osc::IOscCalculator *calc, const SystShifts &shift) const override
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
double POT() const
Definition: SpectrumStan.h:47
Prevent histograms being added to the current directory.
Definition: Utilities.h:66
Float_t e
Definition: plot.C:35
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
Spectrum PredictComponentSyst(osc::IOscCalculator *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:70
Spectrum ShiftSpectrum(const Spectrum &s, CoeffsType type, bool nubar, const SystShifts &shift) const
U _PredictComponentSyst(osc::_IOscCalculator< T > *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Templated helper for PredictComponentSyst.
std::vector< Binning > GetBinnings() const
Definition: Spectrum.h:299
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].
h
Definition: demo3.py:41
const std::vector< Binning > & GetBinnings() const
Definition: SpectrumStan.h:43