Classes | Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
ana::PredictionInterp Class Reference

Implements systematic errors by interpolation between shifted templates. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-10-23/CAFAna/Prediction/PredictionInterp.h"

Inheritance diagram for ana::PredictionInterp:
ana::IPrediction ana::PredictionSyst3Flavor2020 ana::PredictionSystJoint2018 ana::PredictionSystJointDemo ana::PredictionSystNue2017 ana::PredictionSystNumu2017

Classes

struct  Coeffs
 
struct  Key_t
 
struct  ShiftedPreds
 
struct  Val_t
 

Public Types

enum  EMode_t { kCombineSigns, kSplitBySign }
 
enum  CoeffsType {
  kNueApp, kNueSurv, kNumuSurv, kNC,
  kOther, kNCoeffTypes
}
 

Public Member Functions

 PredictionInterp (std::vector< const ISyst * > systs, osc::IOscCalc *osc, const IPredictionGenerator &predGen, Loaders &loaders, const SystShifts &shiftMC=kNoShift, EMode_t mode=kCombineSigns)
 
virtual ~PredictionInterp ()
 
Spectrum Predict (osc::IOscCalc *calc) const override
 
Spectrum Predict (osc::IOscCalcStan *calc) const override
 
Spectrum PredictSyst (osc::IOscCalc *calc, const SystShifts &shift) const override
 
Spectrum PredictSyst (osc::IOscCalcStan *calc, const SystShifts &shift) const override
 
Spectrum PredictComponent (osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
 
Spectrum PredictComponent (osc::IOscCalcStan *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
 
Spectrum PredictComponentSyst (osc::IOscCalc *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
 
Spectrum PredictComponentSyst (osc::IOscCalcStan *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
 
virtual void SaveTo (TDirectory *dir, const std::string &name) const override
 
void MinimizeMemory ()
 
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 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
 
void SetOscSeed (osc::IOscCalc *oscSeed)
 
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
 
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
 
bool SplitBySign () const
 
 PredictionInterp ()
 
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. More...
 
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. More...
 
std::vector< const ISyst * > GetAllSysts () const
 
Spectrum ShiftSpectrum (const Spectrum &s, CoeffsType type, bool nubar, const SystShifts &shift) const
 
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. More...
 
Spectrum ShiftedComponent (osc::IOscCalcStan *calc, const TMD5 *hash, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign, CoeffsType type) const
 
virtual Spectrum PredictUnoscillated () const
 
virtual OscillatableSpectrum ComponentCC (int from, int to) const
 
virtual Spectrum ComponentNCTotal () const
 
virtual Spectrum ComponentNC () const
 
virtual Spectrum ComponentNCAnti () const
 

Static Public Member Functions

static std::unique_ptr< PredictionInterpLoadFrom (TDirectory *dir, const std::string &name)
 
static void LoadFromBody (TDirectory *dir, PredictionInterp *ret, std::vector< const ISyst * > veto={})
 

Protected Member Functions

void InitFits () const
 
void InitFitsHelper (ShiftedPreds &sp, std::vector< std::vector< std::vector< Coeffs >>> &fits, Sign::Sign_t sign) const
 
template<typename T >
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. More...
 
template<typename T >
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. More...
 
template<typename T >
void ShiftBins (unsigned int N, T *arr, CoeffsType type, bool nubar, const SystShifts &shift) const
 Helper for ShiftSpectrum. More...
 

Protected Attributes

std::unique_ptr< IPredictionfPredNom
 The nominal prediction. More...
 
std::unordered_map< const ISyst *, ShiftedPredsfPreds
 
std::unordered_map< const ISyst *, ShiftedPredsfSumPreds
 
osc::IOscCalcfOscOrigin
 The oscillation values we assume when evaluating the coefficients. More...
 
Spectrum fBinning
 Dummy spectrum to provide binning. More...
 
ThreadLocal< std::map< Key_t, Val_t > > fNomCache
 
bool fSplitBySign
 

Detailed Description

Implements systematic errors by interpolation between shifted templates.

Definition at line 23 of file PredictionInterp.h.

Member Enumeration Documentation

Enumerator
kNueApp 
kNueSurv 
kNumuSurv 
kNC 
kOther 

Taus, numu appearance.

kNCoeffTypes 

Definition at line 117 of file PredictionInterp.h.

Enumerator
kCombineSigns 
kSplitBySign 

Definition at line 26 of file PredictionInterp.h.

Constructor & Destructor Documentation

ana::PredictionInterp::PredictionInterp ( std::vector< const ISyst * >  systs,
osc::IOscCalc osc,
const IPredictionGenerator predGen,
Loaders loaders,
const SystShifts shiftMC = kNoShift,
EMode_t  mode = kCombineSigns 
)
Parameters
systsWhat systematics we should be capable of interpolating
oscThe oscillation point to expand around
predGenConstruct an IPrediction from the following information.
loadersThe loaders to be passed on to the underlying prediction
shiftMCUnderlying shift. Use with care. Mostly for PredictionNumuFAHadE. Should not contain any of of systs
modeIn FHC the wrong-sign has bad stats and probably the fits can't be split out reasonably. For RHC it's important not to conflate them.

Definition at line 61 of file PredictionInterp.cxx.

References fPredNom, fPreds, ana::IPredictionGenerator::Generate(), loaders, ana::SystShifts::SetShift(), ana::PredictionInterp::ShiftedPreds::shifts, sigma(), and ana::PredictionInterp::ShiftedPreds::systName.

67  : fOscOrigin(osc ? osc->Copy() : 0),
70  {
71  for(const ISyst* syst: systs){
72  ShiftedPreds sp;
73  sp.systName = syst->ShortName();
74 
75  // Genie reweight shifts aren't able to provide +/-3 sigma
76  if(syst->IsGenieReweight())
77  sp.shifts = {-2, -1, 0, +1, +2};
78  else
79  sp.shifts = {-3, -2, -1, 0, +1, +2, +3};
80 
81  for(int sigma: sp.shifts){
82  SystShifts shiftHere = shiftMC;
83  shiftHere.SetShift(syst, sigma);
84  sp.preds.emplace_back(predGen.Generate(loaders, shiftHere).release());
85  }
86 
87  fPreds.emplace(syst, std::move(sp));
88  } // end for syst
89 
90  fPredNom = predGen.Generate(loaders, shiftMC);
91  }
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
static Spectrum Uninitialized()
Definition: Spectrum.h:145
Spectrum fBinning
Dummy spectrum to provide binning.
double sigma(TH1F *hist, double percentile)
virtual _IOscCalc * Copy() const =0
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
osc::IOscCalc * fOscOrigin
The oscillation values we assume when evaluating the coefficients.
ana::PredictionInterp::~PredictionInterp ( )
virtual

Definition at line 94 of file PredictionInterp.cxx.

References FitRatios().

95  {
96  }
ana::PredictionInterp::PredictionInterp ( )
inline

Definition at line 123 of file PredictionInterp.h.

References dir, LoadFromBody(), and runNovaSAM::ret.

static Spectrum Uninitialized()
Definition: Spectrum.h:145
Spectrum fBinning
Dummy spectrum to provide binning.
osc::IOscCalc * fOscOrigin
The oscillation values we assume when evaluating the coefficients.

Member Function Documentation

template<typename T >
Spectrum ana::PredictionInterp::_PredictComponentSyst ( osc::_IOscCalc< T > *  calc,
const SystShifts shift,
Flavors::Flavors_t  flav,
Current::Current_t  curr,
Sign::Sign_t  sign 
) const
protected

Templated helper for PredictComponentSyst.

Definition at line 529 of file PredictionInterp.cxx.

References ana::SystShifts::ActiveSysts(), ana::assert(), om::cerr, ana::Spectrum::Clear(), allTimeWatchdog::endl, fBinning, fPreds, osc::_IOscCalc< T >::GetParamsHash(), samweb_client.utility::hash, InitFits(), ana::Flavors::kAll, ana::Current::kCC, ana::Current::kNC, kNC, kNueApp, kNueSurv, ana::Flavors::kNuEToNuE, ana::Flavors::kNuEToNuMu, ana::Flavors::kNuEToNuTau, kNumuSurv, ana::Flavors::kNuMuToNuE, ana::Flavors::kNuMuToNuMu, ana::Flavors::kNuMuToNuTau, kOther, runNovaSAM::ret, and ShiftedComponent().

Referenced by PredictComponentSyst().

534  {
535  InitFits();
536 
538  ret.Clear();
539 
540  // Check that we're able to handle all the systs we were passed
541  for(const ISyst* syst: shift.ActiveSysts()){
542  if(fPreds.find(syst) == fPreds.end()){
543  std::cerr << "This PredictionInterp is not set up to handle the requested systematic: " << syst->ShortName() << std::endl;
544  abort();
545  }
546  } // end for syst
547 
548 
549  const TMD5* hash = calc ? calc->GetParamsHash() : 0;
550 
551  if(curr & Current::kCC){
552  if(flav & Flavors::kNuEToNuE) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuE, Current::kCC, sign, kNueSurv);
553  if(flav & Flavors::kNuEToNuMu) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuMu, Current::kCC, sign, kOther );
554  if(flav & Flavors::kNuEToNuTau) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuTau, Current::kCC, sign, kOther );
555 
556  if(flav & Flavors::kNuMuToNuE) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuE, Current::kCC, sign, kNueApp );
557  if(flav & Flavors::kNuMuToNuMu) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuMu, Current::kCC, sign, kNumuSurv);
558  if(flav & Flavors::kNuMuToNuTau) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuTau, Current::kCC, sign, kOther );
559  }
560  if(curr & Current::kNC){
561  assert(flav == Flavors::kAll); // Don't know how to calculate anything else
562 
563  ret += ShiftedComponent(calc, hash, shift, Flavors::kAll, Current::kNC, sign, kNC);
564  }
565 
566  delete hash;
567 
568  return ret;
569  }
(&#39; appearance&#39;)
Definition: IPrediction.h:18
(&#39;beam &#39;)
Definition: IPrediction.h:15
OStream cerr
Definition: OStream.cxx:7
void Clear()
Definition: Spectrum.cxx:372
Charged-current interactions.
Definition: IPrediction.h:39
Spectrum fBinning
Dummy spectrum to provide binning.
std::vector< float > Spectrum
Definition: Constants.h:570
virtual TMD5 * GetParamsHash() const
Use to check two calculators are in the same state.
Definition: IOscCalc.h:34
(&#39; survival&#39;)
Definition: IPrediction.h:19
Taus, numu appearance.
Neutral-current interactions.
Definition: IPrediction.h:40
assert(nhit_max >=nhit_nbins)
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
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.
def sign(x)
Definition: canMan.py:197
template<typename T >
Spectrum ana::PredictionInterp::_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
protected

Templated helper for ShiftedComponent.

Definition at line 477 of file PredictionInterp.cxx.

References febshutoff_auto::curr, fNomCache, fPredNom, fSplitBySign, samweb_client.utility::hash, it, ana::Sign::kAntiNu, ana::Sign::kBoth, findDuplicateFiles::key, ana::Sign::kNu, ShiftedComponent(), and ShiftSpectrum().

Referenced by ShiftedComponent().

484  {
486  "PredictionInterp::ShiftedComponent() can only be called using doubles or stan::math::vars");
487 
488  if(fSplitBySign && sign == Sign::kBoth){
489  return (ShiftedComponent(calc, hash, shift, flav, curr, Sign::kAntiNu, type)+
490  ShiftedComponent(calc, hash, shift, flav, curr, Sign::kNu, type));
491  }
492 
493  // Must be the base case of the recursion to use the cache. Otherwise we
494  // can cache systematically shifted versions of our children, which is
495  // wrong. Also, some calculators won't hash themselves.
496  // Moreover, caching is not going to work with stan::math::vars
497  // since they get reset every time Stan's log_prob() is called.
498  const bool canCache = (hash != 0) && !std::is_same<T, stan::math::var>::value;
499 
500  const Key_t key = {flav, curr, sign};
501  auto it = fNomCache->find(key);
502 
503  // Should the interpolation use the nubar fits?
504  const bool nubar = (fSplitBySign && sign == Sign::kAntiNu);
505 
506  // We have the nominal for this exact combination of flav, curr, sign, calc
507  // stored. Shift it and return.
508  if(canCache && it != fNomCache->end() && it->second.hash == *hash){
509  return ShiftSpectrum(it->second.nom, type, nubar, shift);
510  }
511 
512  // We need to compute the nominal again for whatever reason
513  const Spectrum nom = fPredNom->PredictComponent(calc, flav, curr, sign);
514 
515  if(canCache){
516  // Insert into the cache if not already there, or update if there but
517  // with old oscillation parameters.
518  if(it == fNomCache->end())
519  fNomCache->emplace(key, Val_t({*hash, nom}));
520  else
521  it->second = {*hash, nom};
522  }
523 
524  return ShiftSpectrum(nom, type, nubar, shift);
525  }
set< int >::iterator it
Antineutrinos-only.
Definition: IPrediction.h:50
ThreadLocal< std::map< Key_t, Val_t > > fNomCache
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
const XML_Char int const XML_Char * value
Definition: expat.h:331
std::vector< float > Spectrum
Definition: Constants.h:570
Neutrinos-only.
Definition: IPrediction.h:49
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
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.
Spectrum ShiftSpectrum(const Spectrum &s, CoeffsType type, bool nubar, const SystShifts &shift) const
def sign(x)
Definition: canMan.py:197
OscillatableSpectrum ana::IPrediction::ComponentCC ( int  from,
int  to 
) const
virtualinherited

Reimplemented in ana::PredictionExtendToPeripheral, ana::PredictionAddRock, ana::PredictionCombinePeriods, ana::PredictionExtrapSum, and ana::PredictionExtrap.

Definition at line 105 of file IPrediction.cxx.

References om::cout, and allTimeWatchdog::endl.

Referenced by ana::PredictionCombinePeriods::ComponentCC(), ana::PredictionAddRock::ComponentCC(), and ND_DataMC().

106  {
107  std::cout << "WARNING! ComponentCC is unimplemented in IPrediction" << std::endl; abort();
108  }
OStream cout
Definition: OStream.cxx:6
Spectrum ana::IPrediction::ComponentNC ( ) const
virtualinherited

Reimplemented in ana::PredictionExtendToPeripheral, ana::PredictionAddRock, ana::PredictionCombinePeriods, ana::PredictionExtrapSum, and ana::PredictionExtrap.

Definition at line 115 of file IPrediction.cxx.

References om::cout, and allTimeWatchdog::endl.

Referenced by ana::PredictionAddRock::ComponentNC().

116  {
117  std::cout << "WARNING! ComponentNC is unimplemented in IPrediction" << std::endl; abort();
118  }
OStream cout
Definition: OStream.cxx:6
Spectrum ana::IPrediction::ComponentNCAnti ( ) const
virtualinherited

Reimplemented in ana::PredictionExtendToPeripheral, ana::PredictionAddRock, ana::PredictionCombinePeriods, ana::PredictionExtrapSum, and ana::PredictionExtrap.

Definition at line 120 of file IPrediction.cxx.

References om::cout, and allTimeWatchdog::endl.

Referenced by ana::PredictionAddRock::ComponentNCAnti().

121  {
122  std::cout << "WARNING! ComponentNCAnti is unimplemented in IPrediction" << std::endl; abort();
123  }
OStream cout
Definition: OStream.cxx:6
Spectrum ana::IPrediction::ComponentNCTotal ( ) const
virtualinherited

Reimplemented in ana::PredictionExtendToPeripheral, ana::PredictionAddRock, ana::PredictionCombinePeriods, ana::PredictionExtrapSum, and ana::PredictionExtrap.

Definition at line 110 of file IPrediction.cxx.

References om::cout, and allTimeWatchdog::endl.

Referenced by ana::PredictionAddRock::ComponentNCTotal(), and ND_DataMC().

111  {
112  std::cout << "WARNING! ComponentNCTotal is unimplemented in IPrediction" << std::endl; abort();
113  }
OStream cout
Definition: OStream.cxx:6
void ana::PredictionInterp::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

Definition at line 722 of file PredictionInterp.cxx.

References bin, plot_validation_datamc::c, om::cout, Draw(), allTimeWatchdog::endl, genie::utils::style::Format(), fPredNom, fPreds, make_syst_table_plots::h, MECModelEnuComparisons::i, InitFits(), makeTrainCVSamples::int, std::isnan(), it, nbins, ana::IPrediction::PredictComponent(), PredictComponentSyst(), ratio(), ana::ISyst::ShortName(), std::sqrt(), ss, ana::Spectrum::ToTH1(), and submit_syst::x.

Referenced by DebugPlots().

727  {
728  InitFits();
729 
730  auto it = fPreds.find(syst);
731  if(it == fPreds.end()){
732  std::cout << "PredictionInterp::DebugPlots(): "
733  << syst->ShortName() << " not found" << std::endl;
734  return;
735  }
736 
737  std::unique_ptr<TH1> nom(fPredNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
738  const int nbins = nom->GetNbinsX();
739 
740  TGraph* curves[nbins];
741  TGraph* points[nbins];
742 
743  for(int i = 0; i <= 80; ++i){
744  const double x = .1*i-4;
745  const SystShifts ss(it->first, x);
746  std::unique_ptr<TH1> h(PredictComponentSyst(calc, ss, flav, curr, sign).ToTH1(18e20));
747 
748  for(int bin = 0; bin < nbins; ++bin){
749  if(i == 0){
750  curves[bin] = new TGraph;
751  points[bin] = new TGraph;
752  }
753 
754  const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
755 
756  if(!std::isnan(ratio)) curves[bin]->SetPoint(curves[bin]->GetN(), x, ratio);
757  else curves[bin]->SetPoint(curves[bin]->GetN(), x, 1);
758  } // end for bin
759  } // end for i (x)
760 
761  // As elswhere, to allow BirksC etc that need a different nominal to plot
762  // right.
763  IPrediction* pNom = 0;
764  for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
765  if(it->second.shifts[shiftIdx] == 0) pNom = it->second.preds[shiftIdx];;
766  }
767  if(pNom){ // if not, probably MinimizeMemory() was called
768  std::unique_ptr<TH1> hnom(pNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
769 
770  for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
771  if(!it->second.preds[shiftIdx]) continue; // Probably MinimizeMemory()
772  std::unique_ptr<TH1> h;
773  h = std::move(std::unique_ptr<TH1>(it->second.preds[shiftIdx]->PredictComponent(calc, flav, curr, sign).ToTH1(18e20)));
774 
775  for(int bin = 0; bin < nbins; ++bin){
776  const double ratio = h->GetBinContent(bin+1)/hnom->GetBinContent(bin+1);
777  if(!std::isnan(ratio)) points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], ratio);
778  else points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], 1);
779  }
780  } // end for shiftIdx
781  } // end if pNom
782 
783 
784  int nx = int(sqrt(nbins));
785  int ny = int(sqrt(nbins));
786  if(nx*ny < nbins) ++nx;
787  if(nx*ny < nbins) ++ny;
788 
789  TCanvas* c = new TCanvas;
790  c->Divide(nx, ny);
791 
792  for(int bin = 0; bin < nbins; ++bin){
793  c->cd(bin+1);
794  (new TH2F("",
795  TString::Format("%s %g < %s < %g;Shift;Ratio",
796  it->second.systName.c_str(),
797  nom->GetXaxis()->GetBinLowEdge(bin+1),
798  nom->GetXaxis()->GetTitle(),
799  nom->GetXaxis()->GetBinUpEdge(bin+1)),
800  100, -4, +4, 100, .5, 1.5))->Draw();
801  curves[bin]->Draw("l same");
802  points[bin]->SetMarkerStyle(kFullDotMedium);
803  points[bin]->Draw("p same");
804  } // end for bin
805 
806  c->cd(0);
807  }
tree Draw("slc.nhit")
set< int >::iterator it
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ratio(TH1 *h1, TH1 *h2)
Float_t ss
Definition: plot.C:24
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
const int nbins
Definition: cellShifts.C:15
Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
def sign(x)
Definition: canMan.py:197
void ana::PredictionInterp::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

Definition at line 827 of file PredictionInterp.cxx.

References bin, fPredNom, make_syst_table_plots::h, h2, MECModelEnuComparisons::i, InitFits(), std::isinf(), std::isnan(), ana::ISyst::LatexName(), nbins, PredictComponentSyst(), ratio(), ss, and submit_syst::y.

Referenced by DebugPlotsColz().

832  {
833  InitFits();
834 
835  std::unique_ptr<TH1> nom(fPredNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
836  const int nbins = nom->GetNbinsX();
837 
838  TH2* h2 = new TH2F("", (syst->LatexName()+";;#sigma").c_str(),
839  nbins, nom->GetXaxis()->GetXmin(), nom->GetXaxis()->GetXmax(),
840  80, -4, +4);
841  h2->GetXaxis()->SetTitle(nom->GetXaxis()->GetTitle());
842 
843  for(int i = 1; i <= 80; ++i){
844  const double y = h2->GetYaxis()->GetBinCenter(i);
845  const SystShifts ss(syst, y);
846  std::unique_ptr<TH1> h(PredictComponentSyst(calc, ss, flav, curr, sign).ToTH1(18e20));
847 
848  for(int bin = 0; bin < nbins; ++bin){
849  const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
850 
851  if(!isnan(ratio) && !isinf(ratio))
852  h2->Fill(h2->GetXaxis()->GetBinCenter(bin), y, ratio);
853  } // end for bin
854  } // end for i (x)
855 
856  h2->Draw("colz");
857  h2->SetMinimum(0.5);
858  h2->SetMaximum(1.5);
859  }
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
TH1 * ratio(TH1 *h1, TH1 *h2)
Float_t ss
Definition: plot.C:24
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
const int nbins
Definition: cellShifts.C:15
Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
float bin[41]
Definition: plottest35.C:14
TH1F * h2
Definition: plot.C:45
def sign(x)
Definition: canMan.py:197
void ana::PredictionInterp::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

Definition at line 810 of file PredictionInterp.cxx.

References ana::assert(), DebugPlot(), genie::utils::style::Format(), fPreds, and it.

Referenced by check_predinterp(), check_predinterp_numu(), and test_prediction_interp().

815  {
816  for(auto& it: fPreds){
817  DebugPlot(it.first, calc, flav, curr, sign);
818 
819  if(!savePattern.empty()){
820  assert(savePattern.find("%s") != std::string::npos);
821  gPad->Print(TString::Format(savePattern.c_str(), it.second.systName.c_str()).Data());
822  }
823  } // end for it
824  }
set< int >::iterator it
assert(nhit_max >=nhit_nbins)
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 Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
def sign(x)
Definition: canMan.py:197
void ana::PredictionInterp::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

Definition at line 862 of file PredictionInterp.cxx.

References ana::assert(), DebugPlotColz(), genie::utils::style::Format(), fPreds, InitFits(), and it.

867  {
868  InitFits();
869 
870  for(auto& it: fPreds){
871  new TCanvas;
872  DebugPlotColz(it.first, calc, flav, curr, sign);
873 
874  if(!savePattern.empty()){
875  assert(savePattern.find("%s") != std::string::npos);
876  gPad->Print(TString::Format(savePattern.c_str(), it.second.systName.c_str()).Data());
877  }
878  } // end for it
879  }
set< int >::iterator it
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
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
def sign(x)
Definition: canMan.py:197
std::vector< std::vector< PredictionInterp::Coeffs > > ana::PredictionInterp::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.

Definition at line 200 of file PredictionInterp.cxx.

References ana::assert(), om::cout, allTimeWatchdog::endl, FitRatios(), fOscOrigin, MECModelEnuComparisons::i, ana::IPrediction::PredictComponent(), and r().

Referenced by FitRatios(), and InitFitsHelper().

206  {
207  IPrediction* pNom = 0;
208  for(unsigned int i = 0; i < shifts.size(); ++i){
209  if(shifts[i] == 0) {pNom = preds[i]; };
210  }
211  assert(pNom);
212 
213  // Do it this way rather than via fPredNom so that systematics evaluated
214  // relative to some alternate nominal (eg Birks C where the appropriate
215  // nominal is no-rock) can work.
216  const Spectrum nom = pNom->PredictComponent(fOscOrigin,
217  flav, curr, sign);
218 
219  std::vector<Eigen::ArrayXd> ratios;
220  ratios.reserve(preds.size());
221  for(auto& p: preds){
222  ratios.emplace_back(Ratio(p->PredictComponent(fOscOrigin,
223  flav, curr, sign),
224  nom).GetEigen());
225 
226  // Check none of the ratio values is crazy
227  Eigen::ArrayXd& r = ratios.back();
228  for(int i = 0; i < r.size(); ++i){
229  if(r[i] > 500){
230  std::cout << "PredictionInterp: WARNING, ratio in bin "
231  << i << " for " << shifts[&p-&preds.front()]
232  << " sigma shift of " << systName
233  << " is " << r[i] << ". Ignoring." << std::endl;
234  r[i] = 1;
235  }
236  }
237  }
238 
239  return FitRatios(shifts, ratios);
240  }
const char * p
Definition: xmltok.h:285
std::vector< float > Spectrum
Definition: Constants.h:570
OStream cout
Definition: OStream.cxx:6
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
osc::IOscCalc * fOscOrigin
The oscillation values we assume when evaluating the coefficients.
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.
def sign(x)
Definition: canMan.py:197
std::vector< std::vector< PredictionInterp::Coeffs > > ana::PredictionInterp::FitRatios ( const std::vector< double > &  shifts,
const std::vector< Eigen::ArrayXd > &  ratios 
) const

Find coefficients describing this set of shifts.

Definition at line 100 of file PredictionInterp.cxx.

References ana::assert(), plot_validation_datamc::c, om::cout, lem_server::cs, util::cube(), e, allTimeWatchdog::endl, stan::math::fabs(), FitComponent(), MECModelEnuComparisons::i, m, runNovaSAM::ret, util::sqr(), update_sam_good_runs_metadata::stride, registry_explorer::v, y1, and submit_syst::y2.

Referenced by FitComponent(), and ~PredictionInterp().

102  {
103  if(ratios.size() < 2){
104  std::cout << "PredictionInterp::FitRatios(): ratios.size() = " << ratios.size() << " - how did that happen?" << std::endl;
105 
106  abort();
107  }
108 
109  assert(shifts.size() == ratios.size());
110 
111  std::vector<std::vector<Coeffs>> ret;
112 
113  const int binMax = ratios[0].size();
114 
115  for(int binIdx = 0; binIdx < binMax; ++binIdx){
116  ret.push_back({});
117 
118  // This is cubic interpolation. For each adjacent set of four points we
119  // determine coefficients for a cubic which will be the curve between the
120  // center two. We constrain the function to match the two center points
121  // and to have the right mean gradient at them. This causes this patch to
122  // match smoothly with the next one along. The resulting function is
123  // continuous and first and second differentiable. At the ends of the
124  // range we fit a quadratic instead with only one constraint on the
125  // slope. The coordinate conventions are that point y1 sits at x=0 and y2
126  // at x=1. The matrices are simply the inverses of writing out the
127  // constraints expressed above.
128 
129  // Special-case for linear interpolation
130  if(ratios.size() == 2){
131  const double y0 = ratios[0][binIdx];
132  const double y1 = ratios[1][binIdx];
133 
134  ret.back().emplace_back(0, 0, y1-y0, y0);
135  continue;
136  }
137 
138  {
139  const double y1 = ratios[0][binIdx];
140  const double y2 = ratios[1][binIdx];
141  const double y3 = ratios[2][binIdx];
142  const double v[3] = {y1, y2, (y3-y1)/2};
143  const double m[9] = { 1, -1, 1,
144  -2, 2, -1,
145  1, 0, 0};
146  const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
147  ret.back().emplace_back(0, res(0), res(1), res(2));
148  }
149 
150  // We're assuming here that the shifts are separated by exactly 1 sigma.
151  for(unsigned int shiftIdx = 1; shiftIdx < ratios.size()-2; ++shiftIdx){
152  const double y0 = ratios[shiftIdx-1][binIdx];
153  const double y1 = ratios[shiftIdx ][binIdx];
154  const double y2 = ratios[shiftIdx+1][binIdx];
155  const double y3 = ratios[shiftIdx+2][binIdx];
156 
157  const double v[4] = {y1, y2, (y2-y0)/2, (y3-y1)/2};
158  const double m[16] = { 2, -2, 1, 1,
159  -3, 3, -2, -1,
160  0, 0, 1, 0,
161  1, 0, 0, 0};
162  const TVectorD res = TMatrixD(4, 4, m) * TVectorD(4, v);
163  ret.back().emplace_back(res(0), res(1), res(2), res(3));
164  } // end for shiftIdx
165 
166  {
167  const int N = ratios.size()-3;
168  const double y0 = ratios[N ][binIdx];
169  const double y1 = ratios[N+1][binIdx];
170  const double y2 = ratios[N+2][binIdx];
171  const double v[3] = {y1, y2, (y2-y0)/2};
172  const double m[9] = {-1, 1, -1,
173  0, 0, 1,
174  1, 0, 0};
175  const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
176  ret.back().emplace_back(0, res(0), res(1), res(2));
177  }
178  } // end for binIdx
179 
180  double stride = -1;
181  for(unsigned int i = 0; i < shifts.size()-1; ++i){
182  const double newStride = shifts[i+1]-shifts[i];
183  assert((stride < 0 || fabs(stride-newStride) < 1e-3) &&
184  "Variably-spaced syst templates are unsupported");
185  stride = newStride;
186  }
187 
188  // If the stride is actually not 1, need to rescale all the coefficients
189  for(std::vector<Coeffs>& cs: ret)
190  for(Coeffs& c: cs){
191  c = Coeffs(c.a/util::cube(stride),
192  c.b/util::sqr(stride),
193  c.c/stride,
194  c.d);}
195  return ret;
196  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
T cube(T x)
More efficient cube function than pow(x,3)
Definition: MathUtil.h:26
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
TVectorT< double > TVectorD
Definition: Utilities.h:20
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
OStream cout
Definition: OStream.cxx:6
assert(nhit_max >=nhit_nbins)
Float_t e
Definition: plot.C:35
std::vector< const ISyst * > ana::PredictionInterp::GetAllSysts ( ) const

Definition at line 243 of file PredictionInterp.cxx.

References fPreds.

244  {
245  std::vector<const ISyst*> allsysts;
246  for (const auto &p : fPreds)
247  allsysts.push_back(p.first);
248 
249  return allsysts;
250  }
const char * p
Definition: xmltok.h:285
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
void ana::PredictionInterp::InitFits ( ) const
protected

Definition at line 269 of file PredictionInterp.cxx.

References ana::Spectrum::Clear(), fBinning, ana::PredictionInterp::ShiftedPreds::FillRemaps(), ana::PredictionInterp::ShiftedPreds::fits, ana::PredictionInterp::ShiftedPreds::fitsNubar, fOscOrigin, fPredNom, fPreds, fSplitBySign, InitFitsHelper(), it, ana::Sign::kAntiNu, ana::Sign::kBoth, ana::Sign::kNu, ana::Spectrum::Livetime(), ana::PredictionInterp::ShiftedPreds::nCoeffs, and ana::Spectrum::POT().

Referenced by _PredictComponentSyst(), DebugPlot(), DebugPlotColz(), DebugPlotsColz(), ana::PredictionSyst3Flavor2020::PredictionSyst3Flavor2020(), ana::PredictionSystJoint2018::PredictionSystJoint2018(), ana::PredictionSystJointDemo::PredictionSystJointDemo(), ana::PredictionSystNue2017::PredictionSystNue2017(), ana::PredictionSystNumu2017::PredictionSystNumu2017(), PredictSyst(), SaveTo(), and SetOscSeed().

270  {
271  // No systs
272  if(fPreds.empty()){
273  if(fBinning.POT() > 0 || fBinning.Livetime() > 0) return;
274  }
275  // Already initialized
276  else if(!fPreds.empty() && !fPreds.begin()->second.fits.empty()) return;
277 
278  for(auto& it: fPreds){
279  ShiftedPreds& sp = it.second;
280 
281  if(fSplitBySign){
282  InitFitsHelper(sp, sp.fits, Sign::kNu);
283  InitFitsHelper(sp, sp.fitsNubar, Sign::kAntiNu);
284  }
285  else{
286  InitFitsHelper(sp, sp.fits, Sign::kBoth);
287  }
288  sp.nCoeffs = sp.fits[0][0].size();
289 
290  // Copy the outputs into the remapped indexing order. TODO this is very
291  // ugly. Best would be to generate things in this order natively.
292  sp.FillRemaps();
293  }
294 
295  // Predict something, anything, so that we can know what binning to use
296  fBinning = fPredNom->Predict(fOscOrigin);
297  fBinning.Clear();
298  }
set< int >::iterator it
Antineutrinos-only.
Definition: IPrediction.h:50
void Clear()
Definition: Spectrum.cxx:372
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
Spectrum fBinning
Dummy spectrum to provide binning.
double POT() const
Definition: Spectrum.h:219
Neutrinos-only.
Definition: IPrediction.h:49
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
void InitFitsHelper(ShiftedPreds &sp, std::vector< std::vector< std::vector< Coeffs >>> &fits, Sign::Sign_t sign) const
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
osc::IOscCalc * fOscOrigin
The oscillation values we assume when evaluating the coefficients.
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:222
void ana::PredictionInterp::InitFitsHelper ( ShiftedPreds sp,
std::vector< std::vector< std::vector< Coeffs >>> &  fits,
Sign::Sign_t  sign 
) const
protected

Definition at line 253 of file PredictionInterp.cxx.

References FitComponent(), fits, ana::Flavors::kAll, ana::Flavors::kAllNuTau, ana::Current::kCC, ana::Current::kNC, kNC, kNCoeffTypes, kNueApp, kNueSurv, ana::Flavors::kNuEToNuE, ana::Flavors::kNuEToNuMu, kNumuSurv, ana::Flavors::kNuMuToNuE, ana::Flavors::kNuMuToNuMu, kOther, ana::PredictionInterp::ShiftedPreds::preds, ana::PredictionInterp::ShiftedPreds::shifts, and ana::PredictionInterp::ShiftedPreds::systName.

Referenced by InitFits().

256  {
257  fits.resize(kNCoeffTypes);
258 
259  fits[kNueApp] = FitComponent(sp.shifts, sp.preds, Flavors::kNuMuToNuE, Current::kCC, sign, sp.systName);
260  fits[kNueSurv] = FitComponent(sp.shifts, sp.preds, Flavors::kNuEToNuE, Current::kCC, sign, sp.systName);
261  fits[kNumuSurv] = FitComponent(sp.shifts, sp.preds, Flavors::kNuMuToNuMu, Current::kCC, sign, sp.systName);
262 
263  fits[kNC] = FitComponent(sp.shifts, sp.preds, Flavors::kAll, Current::kNC, sign, sp.systName);
264 
265  fits[kOther] = FitComponent(sp.shifts, sp.preds, Flavors::kNuEToNuMu | Flavors::kAllNuTau, Current::kCC, sign, sp.systName);
266  }
(&#39; appearance&#39;)
Definition: IPrediction.h:18
(&#39;beam &#39;)
Definition: IPrediction.h:15
Charged-current interactions.
Definition: IPrediction.h:39
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.
(&#39; survival&#39;)
Definition: IPrediction.h:19
Taus, numu appearance.
std::vector< std::string > fits
Neutral-current interactions.
Definition: IPrediction.h:40
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
def sign(x)
Definition: canMan.py:197
std::unique_ptr< PredictionInterp > ana::PredictionInterp::LoadFrom ( TDirectory *  dir,
const std::string &  name 
)
static

Definition at line 643 of file PredictionInterp.cxx.

References ana::assert(), dir, LoadFromBody(), runNovaSAM::ret, and getGoodRuns4SAM::tag.

Referenced by CalcRWithSystsNus17(), FCContour(), futureSig_reach_singlePOTcombo_syst(), genie_contours(), get_univ_chisq(), LoadPrediction(), make_nus17_fc_surfs(), make_nus_fc_surfs(), makeMatrixElementSurface(), PlotNus17PredSystsData(), and PlotNusSensAna01().

644  {
645  dir = dir->GetDirectory(name.c_str()); // switch to subdir
646  assert(dir);
647 
648  TObjString* tag = (TObjString*)dir->Get("type");
649  assert(tag);
650  assert(tag->GetString() == "PredictionInterp");
651 
652  std::unique_ptr<PredictionInterp> ret(new PredictionInterp);
653 
654  LoadFromBody(dir, ret.get());
655 
656  TObjString* split_sign = (TObjString*)dir->Get("split_sign");
657  // Can be missing from old files
658  ret->fSplitBySign = (split_sign && split_sign->String() == "yes");
659 
660  delete dir;
661 
662  return ret;
663  }
const XML_Char * name
Definition: expat.h:151
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
static void LoadFromBody(TDirectory *dir, PredictionInterp *ret, std::vector< const ISyst * > veto={})
void ana::PredictionInterp::LoadFromBody ( TDirectory *  dir,
PredictionInterp ret,
std::vector< const ISyst * >  veto = {} 
)
static

Definition at line 666 of file PredictionInterp.cxx.

References dir, genie::utils::style::Format(), fOscOrigin, fPredNom, fPreds, ana::LoadFrom< IPrediction >(), ana::PredictionInterp::ShiftedPreds::preds, runNovaSAM::release, ana::PredictionInterp::ShiftedPreds::shifts, ana::Registry< T >::ShortNameToPtr(), and ana::PredictionInterp::ShiftedPreds::systName.

Referenced by ana::PredictionSystJoint2018::LoadFrom(), ana::PredictionSystJointDemo::LoadFrom(), LoadFrom(), ana::PredictionSystNue2017::LoadFrom(), ana::PredictionSystNumu2017::LoadFrom(), and PredictionInterp().

668  {
669  ret->fPredNom = ana::LoadFrom<IPrediction>(dir, "pred_nom");
670 
671  TH1* hSystNames = (TH1*)dir->Get("syst_names");
672  if(hSystNames){
673  for(int systIdx = 0; systIdx < hSystNames->GetNbinsX(); ++systIdx){
674  ShiftedPreds sp;
675  sp.systName = hSystNames->GetXaxis()->GetBinLabel(systIdx+1);
676 
677  const ISyst* syst = Registry<ISyst>::ShortNameToPtr(sp.systName);
678 
679  if(std::find(veto.begin(), veto.end(), syst) != veto.end()) continue;
680 
681  for(int shift = -3; shift <= +3; ++shift){
682  const std::string subname = TString::Format("pred_%s_%+d", sp.systName.c_str(), shift).Data();
683  TDirectory *preddir = dir->GetDirectory(subname.c_str());
684  if(!preddir) continue; // Can happen for genie systs
685  delete preddir;
686 
687  sp.shifts.push_back(shift);
688  sp.preds.emplace_back(ana::LoadFrom<IPrediction>(dir, subname).release());
689  } // end for shift
690 
691  ret->fPreds.emplace(syst, std::move(sp));
692  } // end for systIdx
693  } // end if hSystNames
694 
695  ret->fOscOrigin = ana::LoadFrom<osc::IOscCalc>(dir, "osc_origin").release();
696  }
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:18
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
static const T * ShortNameToPtr(const std::string &s, bool allowFail=false)
Definition: Registry.cxx:60
void ana::PredictionInterp::MinimizeMemory ( )

After calling this DebugPlots won't work fully and SaveTo won't work at all.

Definition at line 699 of file PredictionInterp.cxx.

References fPredNom, fPreds, MECModelEnuComparisons::i, and it.

Referenced by ana::PredictionSyst3Flavor2020::PredictionSyst3Flavor2020(), ana::PredictionSystJoint2018::PredictionSystJoint2018(), and ana::PredictionSystJointDemo::PredictionSystJointDemo().

700  {
701  std::set<IPrediction*> todel;
702  for(auto& it: fPreds){
703  std::vector<IPrediction*>& preds = it.second.preds;
704  for(unsigned int i = 0; i < preds.size(); ++i){
705  if(preds[i] != fPredNom.get()){
706  todel.insert(preds[i]);
707  preds[i] = 0;
708  }
709  }
710  }
711 
712  for(IPrediction* p: todel) delete p;
713 
714  // We probably just freed up a lot of memory, but malloc by default hangs
715  // on to all of it as cache.
716  #ifndef DARWINBUILD
717  malloc_trim(0);
718  #endif
719  }
set< int >::iterator it
const char * p
Definition: xmltok.h:285
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
Spectrum ana::PredictionInterp::Predict ( osc::IOscCalc calc) const
overridevirtual
Spectrum ana::PredictionInterp::Predict ( osc::IOscCalcStan calc) const
overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 314 of file PredictionInterp.cxx.

References fPredNom.

315  {
316  return fPredNom->Predict(calc);
317  }
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
Spectrum ana::PredictionInterp::PredictComponent ( osc::IOscCalc calc,
Flavors::Flavors_t  flav,
Current::Current_t  curr,
Sign::Sign_t  sign 
) const
overridevirtual

Implements ana::IPrediction.

Definition at line 320 of file PredictionInterp.cxx.

References fPredNom.

Referenced by syst_plot_test().

324  {
325  return fPredNom->PredictComponent(calc, flav, curr, sign);
326  }
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
def sign(x)
Definition: canMan.py:197
Spectrum ana::PredictionInterp::PredictComponent ( osc::IOscCalcStan calc,
Flavors::Flavors_t  flav,
Current::Current_t  curr,
Sign::Sign_t  sign 
) const
overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 329 of file PredictionInterp.cxx.

References fPredNom.

333  {
334  return fPredNom->PredictComponent(calc, flav, curr, sign);
335  }
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
def sign(x)
Definition: canMan.py:197
Spectrum ana::PredictionInterp::PredictComponentSyst ( osc::IOscCalc calc,
const SystShifts shift,
Flavors::Flavors_t  flav,
Current::Current_t  curr,
Sign::Sign_t  sign 
) const
overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 573 of file PredictionInterp.cxx.

References _PredictComponentSyst().

Referenced by check_predinterp(), check_predinterp_numu(), DebugPlot(), DebugPlotColz(), and PredictSyst().

578  {
579  return _PredictComponentSyst(calc, shift, flav, curr, sign);
580  }
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.
def sign(x)
Definition: canMan.py:197
Spectrum ana::PredictionInterp::PredictComponentSyst ( osc::IOscCalcStan calc,
const SystShifts shift,
Flavors::Flavors_t  flav,
Current::Current_t  curr,
Sign::Sign_t  sign 
) const
overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 583 of file PredictionInterp.cxx.

References _PredictComponentSyst().

588  {
589  return _PredictComponentSyst(calc, shift, flav, curr, sign);
590  }
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.
def sign(x)
Definition: canMan.py:197
Spectrum ana::PredictionInterp::PredictSyst ( osc::IOscCalc calc,
const SystShifts shift 
) const
overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 338 of file PredictionInterp.cxx.

References InitFits(), ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth, and PredictComponentSyst().

Referenced by check_predinterp(), check_predinterp_numu(), demo::DrawUpDownRatioCanvas(), mec_tuning(), mec_tuning_fitter_2020(), starPlot(), test_stanfit_systpulls(), and TestPred().

340  {
341  InitFits();
342 
343  return PredictComponentSyst(calc, shift,
346  Sign::kBoth);
347  }
Interactions of both types.
Definition: IPrediction.h:42
Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
Spectrum ana::PredictionInterp::PredictSyst ( osc::IOscCalcStan calc,
const SystShifts shift 
) const
overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 350 of file PredictionInterp.cxx.

References InitFits(), ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth, and PredictComponentSyst().

352  {
353  InitFits();
354 
355  return PredictComponentSyst(calc,
356  shift,
359  Sign::kBoth);
360  }
Interactions of both types.
Definition: IPrediction.h:42
Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
Spectrum ana::IPrediction::PredictUnoscillated ( ) const
virtualinherited

Reimplemented in ana::PredictionSterile.

Definition at line 33 of file IPrediction.cxx.

References noosc, and ana::IPrediction::Predict().

Referenced by cc(), demo5(), demo::DrawUpDownRatioCanvas(), efficiency(), efficiencySA(), plot_nd_data_mc(), template_basic(), and test_ana().

34  {
35  // Default implementation
37  return Predict(&noosc);
38  }
Pass neutrinos through unchanged.
Definition: IOscCalc.h:40
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
void ana::PredictionInterp::SaveTo ( TDirectory *  dir,
const std::string &  name 
) const
overridevirtual

Reimplemented from ana::IPrediction.

Reimplemented in ana::PredictionSystNumu2017, ana::PredictionSystNue2017, ana::PredictionSyst3Flavor2020, ana::PredictionSystJointDemo, and ana::PredictionSystJoint2018.

Definition at line 592 of file PredictionInterp.cxx.

References om::cout, dir, allTimeWatchdog::endl, genie::utils::style::Format(), fOscOrigin, fPredNom, fPreds, fSplitBySign, MECModelEnuComparisons::i, InitFits(), makeTrainCVSamples::int, it, ana::PredictionInterp::ShiftedPreds::preds, ana::SaveTo(), ana::PredictionInterp::ShiftedPreds::shifts, ana::PredictionInterp::ShiftedPreds::systName, and tmp.

Referenced by FDDataMCSystBandLoad(), getPredictions(), getPredictions_updatedAna(), make_prediction_rhc(), MakeCovarSim(), MakeNus17PredictionSysts(), MakeNusPredictionSystsAna01(), makePred(), MakePrediction(), MakePredictionNoOsc_FHC_FD(), MakePredictionNoOsc_FHC_ND(), MakePredictionNoOsc_RHC_FD(), MakePredictionNoOsc_RHC_ND(), mec_tuning_preds_2020(), SavePrediction(), savePrediction_complete(), savePrediction_extrap(), savePrediction_systs(), ana::PredictionSystJoint2018::SaveTo(), ana::PredictionSystJointDemo::SaveTo(), ana::PredictionSystNue2017::SaveTo(), and ana::PredictionSystNumu2017::SaveTo().

593  {
594  InitFits();
595 
596  TDirectory* tmp = gDirectory;
597 
598  dir = dir->mkdir(name.c_str()); // switch to subdir
599  dir->cd();
600 
601  TObjString("PredictionInterp").Write("type");
602 
603  fPredNom->SaveTo(dir, "pred_nom");
604 
605  for(auto& it: fPreds){
606  const ShiftedPreds& sp = it.second;
607 
608  for(unsigned int i = 0; i < sp.shifts.size(); ++i){
609  if(!sp.preds[i]){
610  std::cout << "Can't save a PredictionInterp after MinimizeMemory()" << std::endl;
611  abort();
612  }
613 
614  sp.preds[i]->SaveTo(dir, TString::Format("pred_%s_%+d",
615  sp.systName.c_str(),
616  int(sp.shifts[i])).Data());
617  } // end for i
618  } // end for it
619 
620  ana::SaveTo(*fOscOrigin, dir, "osc_origin");
621 
622  if(!fPreds.empty()){
623  TH1F hSystNames("syst_names", ";Syst names", fPreds.size(), 0, fPreds.size());
624  int binIdx = 1;
625  for(auto& it: fPreds){
626  hSystNames.GetXaxis()->SetBinLabel(binIdx++, it.second.systName.c_str());
627  }
628  dir->cd();
629  hSystNames.Write("syst_names");
630  }
631 
632  TObjString split_sign = fSplitBySign ? "yes" : "no";
633  dir->cd();
634  split_sign.Write("split_sign");
635 
636  dir->Write();
637  delete dir;
638 
639  tmp->cd();
640  }
const XML_Char * name
Definition: expat.h:151
set< int >::iterator it
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
Float_t tmp
Definition: plot.C:36
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
OStream cout
Definition: OStream.cxx:6
TDirectory * dir
Definition: macro.C:5
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
osc::IOscCalc * fOscOrigin
The oscillation values we assume when evaluating the coefficients.
void ana::PredictionInterp::SetOscSeed ( osc::IOscCalc oscSeed)

Definition at line 301 of file PredictionInterp.cxx.

References osc::_IOscCalc< T >::Copy(), fOscOrigin, fPreds, InitFits(), and it.

301  {
302  fOscOrigin = oscSeed->Copy();
303  for(auto& it: fPreds) it.second.fits.clear();
304  InitFits();
305  }
set< int >::iterator it
virtual _IOscCalc * Copy() const =0
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
osc::IOscCalc * fOscOrigin
The oscillation values we assume when evaluating the coefficients.
template<typename T >
void ana::PredictionInterp::ShiftBins ( unsigned int  N,
T arr,
CoeffsType  type,
bool  nubar,
const SystShifts shift 
) const
protected

Helper for ShiftSpectrum.

Definition at line 391 of file PredictionInterp.cxx.

References ana::PredictionInterp::Coeffs::a, ana::assert(), ana::PredictionInterp::Coeffs::b, ana::PredictionInterp::Coeffs::c, calib::corr(), om::cout, util::cube(), ana::PredictionInterp::Coeffs::d, allTimeWatchdog::endl, MakeMiniprodValidationCuts::f, fits, ana::PredictionInterp::ShiftedPreds::fitsNubarRemap, ana::PredictionInterp::ShiftedPreds::fitsRemap, fPreds, fSplitBySign, ana::SystShifts::GetShift(), ana::SystShifts::HasAnyStan(), ana::SystShifts::HasStan(), MECModelEnuComparisons::i, it, std::max(), std::min(), getGoodRuns4SAM::n, ana::PredictionInterp::ShiftedPreds::nCoeffs, ana::PredictionInterp::ShiftedPreds::shifts, util::sqr(), ana::PredictionInterp::ShiftedPreds::Stride(), T, and submit_syst::x.

Referenced by ShiftSpectrum().

396  {
397  static_assert(std::is_same_v<T, double> ||
398  std::is_same_v<T, stan::math::var>,
399  "PredictionInterp::ShiftBins() can only be called using doubles or stan::math::vars");
400  if(nubar) assert(fSplitBySign);
401 
402 
403  if(!std::is_same_v<T, stan::math::var> && shift.HasAnyStan()){
404  std::cout << "PredictionInterp: stan shifts on non-stan spectrum, something is wrong" << std::endl;
405  abort();
406  }
407 
408  T corr[N];
409  for(unsigned int i = 0; i < N; ++i) corr[i] = 1;
410 
411  for(auto& it: fPreds){
412  const ISyst* syst = it.first;
413  const ShiftedPreds& sp = it.second;
414 
415  T x = shift.GetShift<T>(syst);
416 
417  // need to actually do the calculation for the autodiff version
418  // to make sure the gradient is computed correctly
419  if(x == 0 && (!std::is_same_v<T, stan::math::var> || !shift.HasStan(syst))) continue;
420 
421  int shiftBin = (util::GetValAs<double>(x) - sp.shifts[0])/sp.Stride();
422  shiftBin = std::max(0, shiftBin);
423  shiftBin = std::min(shiftBin, sp.nCoeffs - 1);
424 
425  const Coeffs* fits = nubar ? &sp.fitsNubarRemap[type][shiftBin].front() : &sp.fitsRemap[type][shiftBin].front();
426 
427  x -= sp.shifts[shiftBin];
428 
429  const T x_cube = util::cube(x);
430  const T x_sqr = util::sqr(x);
431 
432  for(unsigned int n = 0; n < N; ++n){
433  // Uncomment to debug crashes in this function
434  // assert(type < fits.size());
435  // assert(n < sp.fits[type].size());
436  // assert(shiftBin < int(sp.fits[type][n].size()));
437 
438  const Coeffs& f = fits[n];
439 
440  corr[n] *= f.a*x_cube + f.b*x_sqr + f.c*x + f.d;
441  } // end for n
442  } // end for syst
443 
444  for(unsigned int n = 0; n < N; ++n){
445  // std::max() doesn't work with stan::math::var
446  arr[n] *= (corr[n] > 0.) ? corr[n] : 0.;
447  }
448 
449  }
T max(const caf::Proxy< T > &a, T b)
set< int >::iterator it
double corr(TGraph *g, double thresh)
T cube(T x)
More efficient cube function than pow(x,3)
Definition: MathUtil.h:26
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
OStream cout
Definition: OStream.cxx:6
std::vector< std::string > fits
assert(nhit_max >=nhit_nbins)
double T
Definition: Xdiff_gwt.C:5
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
T min(const caf::Proxy< T > &a, T b)
Spectrum ana::PredictionInterp::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.

Definition at line 452 of file PredictionInterp.cxx.

References _ShiftedComponent(), calc, febshutoff_auto::curr, samweb_client.utility::hash, and canMan::sign().

Referenced by _PredictComponentSyst(), and _ShiftedComponent().

459  {
460  return _ShiftedComponent(calc, hash, shift, flav, curr, sign, type);
461  }
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.
def sign(x)
Definition: canMan.py:197
Spectrum ana::PredictionInterp::ShiftedComponent ( osc::IOscCalcStan calc,
const TMD5 *  hash,
const SystShifts shift,
Flavors::Flavors_t  flav,
Current::Current_t  curr,
Sign::Sign_t  sign,
CoeffsType  type 
) const
Spectrum ana::PredictionInterp::ShiftSpectrum ( const Spectrum s,
CoeffsType  type,
bool  nubar,
const SystShifts shift 
) const

Definition at line 363 of file PredictionInterp.cxx.

References ana::SystShifts::ActiveSysts(), om::cout, allTimeWatchdog::endl, ana::Spectrum::GetBinnings(), ana::Spectrum::GetEigen(), ana::Spectrum::GetEigenStan(), ana::Spectrum::GetLabels(), ana::SystShifts::GetShift(), ana::SystShifts::HasAnyStan(), ana::Spectrum::HasStan(), ana::SystShifts::IsNominal(), ana::Spectrum::Livetime(), ana::Spectrum::POT(), and ShiftBins().

Referenced by _ShiftedComponent().

367  {
368  if(s.HasStan() && (!shift.IsNominal() && !shift.HasAnyStan())){
369  std::cout << "PredictionInterp used with Stan oscillations but non-Stan syst shifts. Is that what you expected?" << std::endl;
370  std::cout << " Systs set: " << std::endl;
371  for (const auto & syst : shift.ActiveSysts())
372  std::cout << " " << syst->ShortName() << " = " << shift.GetShift(syst) << std::endl;
373  }
374 
375  if(s.HasStan() || shift.HasAnyStan()){
376  // Need the second case for the NC component
377  Eigen::ArrayXstan vec = s.HasStan() ? s.GetEigenStan(s.POT()) : Eigen::ArrayXstan(s.GetEigen(s.POT()));
378 
379  ShiftBins(vec.size(), vec.data(), type, nubar, shift);
380  return Spectrum(std::move(vec), HistAxis(s.GetLabels(), s.GetBinnings()), s.POT(), s.Livetime());
381  }
382  else{
383  Eigen::ArrayXd vec = s.GetEigen(s.POT());
384  ShiftBins(vec.size(), vec.data(), type, nubar, shift);
385  return Spectrum(std::move(vec), HistAxis(s.GetLabels(), s.GetBinnings()), s.POT(), s.Livetime());
386  }
387  }
void ShiftBins(unsigned int N, T *arr, CoeffsType type, bool nubar, const SystShifts &shift) const
Helper for ShiftSpectrum.
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
const XML_Char * s
Definition: expat.h:262
Eigen::VectorXd vec
std::vector< float > Spectrum
Definition: Constants.h:570
OStream cout
Definition: OStream.cxx:6
Eigen::Array< stan::math::var, Eigen::Dynamic, 1 > ArrayXstan
Definition: StanUtils.h:7
bool ana::PredictionInterp::SplitBySign ( ) const
inline

Definition at line 116 of file PredictionInterp.h.

References fSplitBySign.

116 {return fSplitBySign;}

Member Data Documentation

Spectrum ana::PredictionInterp::fBinning
mutableprotected

Dummy spectrum to provide binning.

Definition at line 205 of file PredictionInterp.h.

Referenced by _PredictComponentSyst(), and InitFits().

ThreadLocal<std::map<Key_t, Val_t> > ana::PredictionInterp::fNomCache
mutableprotected

Definition at line 223 of file PredictionInterp.h.

Referenced by _ShiftedComponent().

osc::IOscCalc* ana::PredictionInterp::fOscOrigin
protected
std::unique_ptr<IPrediction> ana::PredictionInterp::fPredNom
protected
std::unordered_map<const ISyst*, ShiftedPreds> ana::PredictionInterp::fPreds
mutableprotected
bool ana::PredictionInterp::fSplitBySign
protected
std::unordered_map<const ISyst*, ShiftedPreds> ana::PredictionInterp::fSumPreds
mutableprotected

Definition at line 200 of file PredictionInterp.h.


The documentation for this class was generated from the following files: