Loads shifted spectra from files. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-09-26/3FlavorAna/Prediction/PredictionSystNumu2017.h"

Inheritance diagram for ana::PredictionSystNumu2017:

Public Types
enum	EMode_t { kCombineSigns, kSplitBySign }

enum	CoeffsType { kNueApp, kNueSurv, kNumuSurv, kNC, kOther, kNCoeffTypes }

Public Member Functions
	PredictionSystNumu2017 (ENumu2017ExtrapType extrap, osc::IOscCalc osc, int WhichQuant=1, const std::vector< const DummyNumu2017Syst > &systs=kNumu2017Systs, const std::string fname="/nova/ana/nu_mu_ana/Ana2017/Predictions/provisional/pred_numu_reduced_v5.root")

virtual	~PredictionSystNumu2017 ()

void	SaveTo (TDirectory *dir, const std::string &name) const override

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

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

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< PredictionSystNumu2017 >	LoadFrom (TDirectory *dir, const std::string &name)

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

Protected Member Functions
	PredictionSystNumu2017 ()

void	AddNormSysts (const std::vector< const DummyNumu2017Syst * > &systs)

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< IPrediction >	fPredNom
	The nominal prediction. More...

std::unordered_map< const ISyst *, ShiftedPreds >	fPreds

std::unordered_map< const ISyst *, ShiftedPreds >	fSumPreds

osc::IOscCalc *	fOscOrigin
	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

Loads shifted spectra from files.

Definition at line 115 of file PredictionSystNumu2017.h.

Member Enumeration Documentation

enum ana::PredictionInterp::CoeffsType

inherited

Enumerator
kNueApp
kNueSurv
kNumuSurv
kNC
kOther	Taus, numu appearance.
kNCoeffTypes

Definition at line 117 of file PredictionInterp.h.

                    {
       kNueApp, kNueSurv, kNumuSurv, kNC,
       kOther, ///< Taus, numu appearance
       kNCoeffTypes
     };

enum ana::PredictionInterp::EMode_t

inherited

Enumerator
kCombineSigns
kSplitBySign

Definition at line 26 of file PredictionInterp.h.

                 {
       kCombineSigns, kSplitBySign
     };

Constructor & Destructor Documentation

ana::PredictionSystNumu2017::PredictionSystNumu2017	(	ENumu2017ExtrapType	extrap,
		osc::IOscCalc *	osc,
		int	WhichQuant = `1`,
		const std::vector< const DummyNumu2017Syst * > &	systs = `kNumu2017Systs`,
		const std::string	fname = `"/nova/ana/nu_mu_ana/Ana2017/Predictions/provisional/pred_numu_reduced_v5.root"`
	)

Definition at line 59 of file PredictionSystNumu2017.cxx.

References AddNormSysts(), osc::_IOscCalc< T >::Copy(), om::cout, allTimeWatchdog::endl, plot_validation_datamc::fname, ana::PredictionInterp::fOscOrigin, ana::PredictionInterp::fPredNom, ana::PredictionInterp::fPreds, ana::PredictionInterp::InitFits(), ana::kNuMu, ana::kNumu2017NormSystBkg, ana::kNumu2017NormSystBoth, ana::kNumu2017NormSystSig, ana::kNumuNoExtrap, ana::PredictionInterp::ShiftedPreds::preds, runNovaSAM::release, ana::PredictionInterp::ShiftedPreds::shifts, sigma(), ana::PredictionInterp::ShiftedPreds::systName, and art::to_string().

     : PredictionInterp()
   {
     fOscOrigin = osc->Copy();
 
     std::string extrapStr = "pred_";
     if (extrap == ENumu2017ExtrapType::kNuMu)         extrapStr += "xp_numu";
     if (extrap == ENumu2017ExtrapType::kNumuNoExtrap) extrapStr += "nxp";
     extrapStr += "_";
 
     std::string VarName = "numu_pred_Energy_Quant"+std::to_string(WhichQuant);
 
     std::cout << "\nLooking at VarName " << VarName << " --> Loading nominals now..." << std::endl;
 
     IPrediction *nom;
     std::string NomVar  = extrapStr+"Nominal/noShift/"+VarName;
     fPredNom = ana::LoadFromFile<IPrediction>(fname, NomVar);
     nom      = ana::LoadFromFile<IPrediction>(fname, NomVar).release();
     // --- The light level syst uses an alternative nominal.  Grab that now.
     std::string LNomVar = extrapStr+"Lightlevel/noShift/"+VarName;
     IPrediction* nomLightLevel = ana::LoadFromFile<IPrediction>(fname, LNomVar).release();
   
     for(const DummyNumu2017Syst* syst: systs){
       if(syst == &kNumu2017NormSystBoth ||
          syst == &kNumu2017NormSystSig ||
          syst == &kNumu2017NormSystBkg) continue;
 
       // Hack.  For sake of distinct names, every syst's fist 8 characters are numu2017
       std::string systName = syst->ShortName();
       systName.erase(0,8);
 
       std::cout << "Loading " << systName << "..." << std::endl;
 
       ShiftedPreds sp;
       sp.systName = systName;
 
       for(int sigma: syst->Shifts()){
         sp.shifts.push_back(sigma);
 
         if (sigma == 0){
           if (systName == "Lightlevel" || systName == "Cherenkov") {
             std::cout << "\tLooking at systName " << systName << ", so pushing back nomLightLevel " << std::endl;
             sp.preds.push_back(nomLightLevel);
           } else {
             sp.preds.push_back(nom);
           }
           continue;
         }
 
         std::string sigmaStr = "/";
         if(sigma == -2) sigmaStr += "minusTwo/";
         if(sigma == -1) sigmaStr += "minusOne/";
         if(sigma == +1) sigmaStr += "plusOne/";
         if(sigma == +2) sigmaStr += "plusTwo/"; 
 
 
         auto shCore = ana::LoadFromFile<PredictionExtrap>  (fname, extrapStr + systName + sigmaStr +VarName).release();
         sp.preds.push_back(shCore);
       } // end for sigma
 
       fPreds.emplace(syst, sp);
     } // end for syst
 
     // Calculate all the ratios etc
     InitFits();
 
     AddNormSysts(systs);
   }

ana::PredictionSystNumu2017::~PredictionSystNumu2017 ( )

virtual

Definition at line 133 of file PredictionSystNumu2017.cxx.

134 {

135 }

ana::PredictionSystNumu2017::PredictionSystNumu2017 ( )

inlineprotected

Definition at line 130 of file PredictionSystNumu2017.h.

References systs.

130 {}

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

protectedinherited

Templated helper for PredictComponentSyst.

Definition at line 526 of file PredictionInterp.cxx.

Referenced by ana::PredictionInterp::PredictComponentSyst().

   {
     InitFits();
 
     Spectrum ret = fBinning;
     ret.Clear();
 
     // Check that we're able to handle all the systs we were passed
     for(const ISyst* syst: shift.ActiveSysts()){
       if(fPreds.find(syst) == fPreds.end()){
         std::cerr << "This PredictionInterp is not set up to handle the requested systematic: " << syst->ShortName() << std::endl;
         abort();
       }
     } // end for syst
 
 
     const TMD5* hash = calc ? calc->GetParamsHash() : 0;
 
     if(curr & Current::kCC){
       if(flav & Flavors::kNuEToNuE)    ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuE,    Current::kCC, sign, kNueSurv);
       if(flav & Flavors::kNuEToNuMu)   ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuMu,   Current::kCC, sign, kOther  );
       if(flav & Flavors::kNuEToNuTau)  ret += ShiftedComponent(calc, hash, shift, Flavors::kNuEToNuTau,  Current::kCC, sign, kOther  );
 
       if(flav & Flavors::kNuMuToNuE)   ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuE,   Current::kCC, sign, kNueApp  );
       if(flav & Flavors::kNuMuToNuMu)  ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuMu,  Current::kCC, sign, kNumuSurv);
       if(flav & Flavors::kNuMuToNuTau) ret += ShiftedComponent(calc, hash, shift, Flavors::kNuMuToNuTau, Current::kCC, sign, kOther   );
     }
     if(curr & Current::kNC){
       assert(flav == Flavors::kAll); // Don't know how to calculate anything else
 
       ret += ShiftedComponent(calc, hash, shift, Flavors::kAll, Current::kNC, sign, kNC);
     }
 
     delete hash;
 
     return ret;
   }

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

protectedinherited

Templated helper for ShiftedComponent.

Definition at line 474 of file PredictionInterp.cxx.

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

Referenced by ana::PredictionInterp::ShiftedComponent().

   {
     static_assert(std::is_same<T, double>::value || std::is_same<T, stan::math::var>::value,
                   "PredictionInterp::ShiftedComponent() can only be called using doubles or stan::math::vars");
 
     if(fSplitBySign && sign == Sign::kBoth){
       return (ShiftedComponent(calc, hash, shift, flav, curr, Sign::kAntiNu, type)+
               ShiftedComponent(calc, hash, shift, flav, curr, Sign::kNu,     type));
     }
 
     // Must be the base case of the recursion to use the cache. Otherwise we
     // can cache systematically shifted versions of our children, which is
     // wrong. Also, some calculators won't hash themselves.
     // Moreover, caching is not going to work with stan::math::vars
     // since they get reset every time Stan's log_prob() is called.
     const bool canCache = (hash != 0) && !std::is_same<T, stan::math::var>::value;
 
     const Key_t key = {flav, curr, sign};
     auto it = fNomCache->find(key);
 
     // Should the interpolation use the nubar fits?
     const bool nubar = (fSplitBySign && sign == Sign::kAntiNu);
 
     // We have the nominal for this exact combination of flav, curr, sign, calc
     // stored.  Shift it and return.
     if(canCache && it != fNomCache->end() && it->second.hash == *hash){
       return ShiftSpectrum(it->second.nom, type, nubar, shift);
     }
 
     // We need to compute the nominal again for whatever reason
     const Spectrum nom = fPredNom->PredictComponent(calc, flav, curr, sign);
 
     if(canCache){
       // Insert into the cache if not already there, or update if there but
       // with old oscillation parameters.
       if(it == fNomCache->end())
         fNomCache->emplace(key, Val_t({*hash, nom}));
       else
         it->second = {*hash, nom};
     }
 
     return ShiftSpectrum(nom, type, nubar, shift);
   }

void ana::PredictionSystNumu2017::AddNormSysts ( const std::vector< const DummyNumu2017Syst * > & systs )

protected

Definition at line 138 of file PredictionSystNumu2017.cxx.

Referenced by PredictionSystNumu2017().

   {
     TH1* h = Predict(fOscOrigin).ToTH1(1);
     const int nBins = h->GetNbinsX();
     delete h;
 
 
     // Got this from Alex. TODO add reference
     // pythag == sum in quadrature
     // leaving Sig and bkg hooks in place, but they are not used anywhere
     const double normBothSyst = 3.5/100;
     const double normBkgSyst = 3.5/100;
     const double normSigSyst = 3.5/100;
 
     // TODO find a way to write these three norm systs as a loop
     //    if(std::count(systs.begin(), systs.end(), &kNumu2017NormSystBoth)){
     if(std::count_if(systs.begin(), systs.end(),
                      [](const ISyst* s){return s->ShortName() == "numu2017NormBoth";})){
       ShiftedPreds normsp;
       normsp.systName = "NormBoth";
       normsp.shifts = {0};
       normsp.preds = {fPredNom.get()};
       for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
         normsp.fits[coeff].resize(nBins+2);
         for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
           normsp.fits[coeff][binIdx].push_back(Coeffs(0, 0, normBothSyst, 1));
         }
       }
       normsp.FillRemaps();
       fPreds.emplace(&kNumu2017NormSystBoth, normsp);
     }
 
 
     //    if(std::count(systs.begin(), systs.end(), &kNumu2017NormSystBkg)){
     if(std::count_if(systs.begin(), systs.end(),
                      [](const ISyst* s){return s->ShortName() == "numu2017normBkg";})){
       ShiftedPreds normbksp;
       normbksp.systName = "NumuNormBkg";
       normbksp.shifts = {0};
       normbksp.preds = {fPredNom.get()};
       for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
         normbksp.fits[coeff].resize(nBins+2);
         for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
           normbksp.fits[coeff][binIdx].push_back(Coeffs(0, 0, (coeff != kNumuSurv) ? normBkgSyst : 0, 1));
         }
       }
       normbksp.FillRemaps();
       fPreds.emplace(&kNumu2017NormSystBkg, normbksp);
     }
 
 
     //    if(std::count(systs.begin(), systs.end(), &kNumu2017NormSystSig)){
     if(std::count_if(systs.begin(), systs.end(),
                      [](const ISyst* s){return s->ShortName() == "numu2017normSig";})){
       ShiftedPreds normsigsp;
       normsigsp.systName = "NumuNormSig";
       normsigsp.shifts ={0};
       normsigsp.preds = {fPredNom.get()};
       for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
         normsigsp.fits[coeff].resize(nBins+2);
         for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
           normsigsp.fits[coeff][binIdx].push_back(Coeffs(0, 0, (coeff == kNumuSurv) ? normSigSyst : 0, 1));
         }
       }
       normsigsp.FillRemaps();
       fPreds.emplace(&kNumu2017NormSystSig, normsigsp);
     }
   }

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 128 of file IPrediction.cxx.

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

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

   {
     std::cout << "WARNING! ComponentCC is unimplemented in IPrediction" << std::endl; abort();
   }

Spectrum ana::IPrediction::ComponentNC ( ) const

virtualinherited

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

Definition at line 138 of file IPrediction.cxx.

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

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

   {
     std::cout << "WARNING! ComponentNC is unimplemented in IPrediction" << std::endl; abort();
   }

Spectrum ana::IPrediction::ComponentNCAnti ( ) const

virtualinherited

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

Definition at line 143 of file IPrediction.cxx.

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

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

   {
     std::cout << "WARNING! ComponentNCAnti is unimplemented in IPrediction" << std::endl; abort();
   }

Spectrum ana::IPrediction::ComponentNCTotal ( ) const

virtualinherited

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

Definition at line 133 of file IPrediction.cxx.

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

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

   {
     std::cout << "WARNING! ComponentNCTotal is unimplemented in IPrediction" << std::endl; abort();
   }

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

inherited

Definition at line 719 of file PredictionInterp.cxx.

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

Referenced by ana::PredictionInterp::DebugPlots().

   {
     InitFits();
 
     auto it = fPreds.find(syst);
     if(it == fPreds.end()){
       std::cout << "PredictionInterp::DebugPlots(): "
                 << syst->ShortName() << " not found" << std::endl;
       return;
     }
 
     std::unique_ptr<TH1> nom(fPredNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
     const int nbins = nom->GetNbinsX();
 
     TGraph* curves[nbins];
     TGraph* points[nbins];
 
     for(int i = 0; i <= 80; ++i){
       const double x = .1*i-4;
       const SystShifts ss(it->first, x);
       std::unique_ptr<TH1> h(PredictComponentSyst(calc, ss, flav, curr, sign).ToTH1(18e20));
 
       for(int bin = 0; bin < nbins; ++bin){
         if(i == 0){
           curves[bin] = new TGraph;
           points[bin] = new TGraph;
         }
 
         const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
 
         if(!std::isnan(ratio)) curves[bin]->SetPoint(curves[bin]->GetN(), x, ratio);
         else curves[bin]->SetPoint(curves[bin]->GetN(), x, 1);
       } // end for bin
     } // end for i (x)
 
     // As elswhere, to allow BirksC etc that need a different nominal to plot
     // right.
     IPrediction* pNom = 0;
     for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
       if(it->second.shifts[shiftIdx] == 0) pNom = it->second.preds[shiftIdx];;
     }
     if(pNom){ // if not, probably MinimizeMemory() was called
       std::unique_ptr<TH1> hnom(pNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
 
       for(unsigned int shiftIdx = 0; shiftIdx < it->second.shifts.size(); ++shiftIdx){
         if(!it->second.preds[shiftIdx]) continue; // Probably MinimizeMemory()
         std::unique_ptr<TH1> h;
         h = std::move(std::unique_ptr<TH1>(it->second.preds[shiftIdx]->PredictComponent(calc, flav, curr, sign).ToTH1(18e20)));
 
         for(int bin = 0; bin < nbins; ++bin){
           const double ratio = h->GetBinContent(bin+1)/hnom->GetBinContent(bin+1);
           if(!std::isnan(ratio)) points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], ratio);
           else points[bin]->SetPoint(points[bin]->GetN(), it->second.shifts[shiftIdx], 1);
         }
       } // end for shiftIdx
     } // end if pNom
 
 
     int nx = int(sqrt(nbins));
     int ny = int(sqrt(nbins));
     if(nx*ny < nbins) ++nx;
     if(nx*ny < nbins) ++ny;
 
     TCanvas* c = new TCanvas;
     c->Divide(nx, ny);
 
     for(int bin = 0; bin < nbins; ++bin){
       c->cd(bin+1);
       (new TH2F("",
                 TString::Format("%s %g < %s < %g;Shift;Ratio",
                                 it->second.systName.c_str(),
                                 nom->GetXaxis()->GetBinLowEdge(bin+1),
                                 nom->GetXaxis()->GetTitle(),
                                 nom->GetXaxis()->GetBinUpEdge(bin+1)),
                 100, -4, +4, 100, .5, 1.5))->Draw();
       curves[bin]->Draw("l same");
       points[bin]->SetMarkerStyle(kFullDotMedium);
       points[bin]->Draw("p same");
     } // end for bin
 
     c->cd(0);
   }

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

inherited

Definition at line 824 of file PredictionInterp.cxx.

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

Referenced by ana::PredictionInterp::DebugPlotsColz().

   {
     InitFits();
 
     std::unique_ptr<TH1> nom(fPredNom->PredictComponent(calc, flav, curr, sign).ToTH1(18e20));
     const int nbins = nom->GetNbinsX();
 
     TH2* h2 = new TH2F("", (syst->LatexName()+";;#sigma").c_str(),
                        nbins, nom->GetXaxis()->GetXmin(), nom->GetXaxis()->GetXmax(),
                        80, -4, +4);
     h2->GetXaxis()->SetTitle(nom->GetXaxis()->GetTitle());
 
     for(int i = 1; i <= 80; ++i){
       const double y = h2->GetYaxis()->GetBinCenter(i);
       const SystShifts ss(syst, y);
       std::unique_ptr<TH1> h(PredictComponentSyst(calc, ss, flav, curr, sign).ToTH1(18e20));
 
       for(int bin = 0; bin < nbins; ++bin){
         const double ratio = h->GetBinContent(bin+1)/nom->GetBinContent(bin+1);
 
         if(!isnan(ratio) && !isinf(ratio))
           h2->Fill(h2->GetXaxis()->GetBinCenter(bin), y, ratio);
       } // end for bin
     } // end for i (x)
 
     h2->Draw("colz");
     h2->SetMinimum(0.5);
     h2->SetMaximum(1.5);
   }

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

inherited

Definition at line 807 of file PredictionInterp.cxx.

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

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

   {
     for(auto& it: fPreds){
       DebugPlot(it.first, calc, flav, curr, sign);
 
       if(!savePattern.empty()){
         assert(savePattern.find("%s") != std::string::npos);
         gPad->Print(TString::Format(savePattern.c_str(), it.second.systName.c_str()).Data());
       }
     } // end for it
   }

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

inherited

Definition at line 859 of file PredictionInterp.cxx.

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

   {
     InitFits();
 
     for(auto& it: fPreds){
       new TCanvas;
       DebugPlotColz(it.first, calc, flav, curr, sign);
 
       if(!savePattern.empty()){
         assert(savePattern.find("%s") != std::string::npos);
         gPad->Print(TString::Format(savePattern.c_str(), it.second.systName.c_str()).Data());
       }
     } // end for it
   }

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

inherited

Find coefficients describing the ratios from this component.

Definition at line 197 of file PredictionInterp.cxx.

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

Referenced by ana::PredictionInterp::FitRatios(), and ana::PredictionInterp::InitFitsHelper().

   {
     IPrediction* pNom = 0;
     for(unsigned int i = 0; i < shifts.size(); ++i){
       if(shifts[i] == 0) {pNom = preds[i]; };
     }
     assert(pNom);
 
     // Do it this way rather than via fPredNom so that systematics evaluated
     // relative to some alternate nominal (eg Birks C where the appropriate
     // nominal is no-rock) can work.
     const Spectrum nom = pNom->PredictComponent(fOscOrigin,
                                                 flav, curr, sign);
 
     std::vector<Eigen::ArrayXd> ratios;
     ratios.reserve(preds.size());
     for(auto& p: preds){
       ratios.emplace_back(Ratio(p->PredictComponent(fOscOrigin,
                                                     flav, curr, sign),
                                 nom).GetEigen());
 
       // Check none of the ratio values is crazy
       Eigen::ArrayXd& r = ratios.back();
       for(int i = 0; i < r.size(); ++i){
         if(r[i] > 500){
           std::cout << "PredictionInterp: WARNING, ratio in bin "
                     << i << " for " << shifts[&p-&preds.front()]
                     << " sigma shift of " << systName
                     << " is " << r[i] << ". Ignoring." << std::endl;
           r[i] = 1;
         }
       }
     }
 
     return FitRatios(shifts, ratios);
   }

std::vector< std::vector< PredictionInterp::Coeffs > > ana::PredictionInterp::FitRatios	(	const std::vector< double > &	shifts,
		const std::vector< Eigen::ArrayXd > &	ratios
	)		const

inherited

Find coefficients describing this set of shifts.

Definition at line 97 of file PredictionInterp.cxx.

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

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

   {
     if(ratios.size() < 2){
       std::cout << "PredictionInterp::FitRatios(): ratios.size() = " << ratios.size() << " - how did that happen?" << std::endl;
 
       abort();
     }
 
     assert(shifts.size() == ratios.size());
 
     std::vector<std::vector<Coeffs>> ret;
 
     const int binMax = ratios[0].size();
 
     for(int binIdx = 0; binIdx < binMax; ++binIdx){
       ret.push_back({});
 
       // This is cubic interpolation. For each adjacent set of four points we
       // determine coefficients for a cubic which will be the curve between the
       // center two. We constrain the function to match the two center points
       // and to have the right mean gradient at them. This causes this patch to
       // match smoothly with the next one along. The resulting function is
       // continuous and first and second differentiable. At the ends of the
       // range we fit a quadratic instead with only one constraint on the
       // slope. The coordinate conventions are that point y1 sits at x=0 and y2
       // at x=1. The matrices are simply the inverses of writing out the
       // constraints expressed above.
 
       // Special-case for linear interpolation
       if(ratios.size() == 2){
         const double y0 = ratios[0][binIdx];
         const double y1 = ratios[1][binIdx];
 
         ret.back().emplace_back(0, 0, y1-y0, y0);
         continue;
       }
 
       {
         const double y1 = ratios[0][binIdx];
         const double y2 = ratios[1][binIdx];
         const double y3 = ratios[2][binIdx];
         const double v[3] = {y1, y2, (y3-y1)/2};
         const double m[9] = { 1, -1,  1,
                              -2,  2, -1,
                               1,  0,  0};
         const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
         ret.back().emplace_back(0, res(0), res(1), res(2));
       }
 
       // We're assuming here that the shifts are separated by exactly 1 sigma.
       for(unsigned int shiftIdx = 1; shiftIdx < ratios.size()-2; ++shiftIdx){
         const double y0 = ratios[shiftIdx-1][binIdx];
         const double y1 = ratios[shiftIdx  ][binIdx];
         const double y2 = ratios[shiftIdx+1][binIdx];
         const double y3 = ratios[shiftIdx+2][binIdx];
 
         const double v[4] = {y1, y2, (y2-y0)/2, (y3-y1)/2};
         const double m[16] = { 2, -2,  1,  1,
                               -3,  3, -2, -1,
                                0,  0,  1,  0,
                                1,  0,  0,  0};
         const TVectorD res = TMatrixD(4, 4, m) * TVectorD(4, v);
         ret.back().emplace_back(res(0), res(1), res(2), res(3));
       } // end for shiftIdx
 
       {
         const int N = ratios.size()-3;
         const double y0 = ratios[N  ][binIdx];
         const double y1 = ratios[N+1][binIdx];
         const double y2 = ratios[N+2][binIdx];
         const double v[3] = {y1, y2, (y2-y0)/2};
         const double m[9] = {-1,  1, -1,
                               0,  0,  1,
                               1,  0,  0};
         const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
         ret.back().emplace_back(0, res(0), res(1), res(2));
       }
     } // end for binIdx
 
     double stride = -1;
     for(unsigned int i = 0; i < shifts.size()-1; ++i){
       const double newStride = shifts[i+1]-shifts[i];
       assert((stride < 0 || fabs(stride-newStride) < 1e-3) &&
              "Variably-spaced syst templates are unsupported");
       stride = newStride;
     }
 
     // If the stride is actually not 1, need to rescale all the coefficients
     for(std::vector<Coeffs>& cs: ret)
       for(Coeffs& c: cs){
         c = Coeffs(c.a/util::cube(stride),
                    c.b/util::sqr(stride),
                    c.c/stride,
                    c.d);}
     return ret;
   }

std::vector< const ISyst * > ana::PredictionInterp::GetAllSysts ( ) const

inherited

Definition at line 240 of file PredictionInterp.cxx.

References ana::PredictionInterp::fPreds.

   {
     std::vector<const ISyst*> allsysts;
     for (const auto &p : fPreds)
       allsysts.push_back(p.first);
 
     return allsysts;
   }

void ana::PredictionInterp::InitFits ( ) const

protectedinherited

Definition at line 266 of file PredictionInterp.cxx.

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

   {
     // No systs
     if(fPreds.empty()){
       if(fBinning.POT() > 0 || fBinning.Livetime() > 0) return;
     }
     // Already initialized
     else if(!fPreds.empty() && !fPreds.begin()->second.fits.empty()) return;
 
     for(auto& it: fPreds){
       ShiftedPreds& sp = it.second;
 
       if(fSplitBySign){
         InitFitsHelper(sp, sp.fits, Sign::kNu);
         InitFitsHelper(sp, sp.fitsNubar, Sign::kAntiNu);
       }
       else{
         InitFitsHelper(sp, sp.fits, Sign::kBoth);
       }
       sp.nCoeffs = sp.fits[0][0].size();
 
       // Copy the outputs into the remapped indexing order. TODO this is very
       // ugly. Best would be to generate things in this order natively.
       sp.FillRemaps();
     }
 
     // Predict something, anything, so that we can know what binning to use
     fBinning = fPredNom->Predict(fOscOrigin);
     fBinning.Clear();
   }

void ana::PredictionInterp::InitFitsHelper	(	ShiftedPreds &	sp,
		std::vector< std::vector< std::vector< Coeffs >>> &	fits,
		Sign::Sign_t	sign
	)		const

protectedinherited

Definition at line 250 of file PredictionInterp.cxx.

Referenced by ana::PredictionInterp::InitFits().

   {
     fits.resize(kNCoeffTypes);
 
     fits[kNueApp]   = FitComponent(sp.shifts, sp.preds, Flavors::kNuMuToNuE,  Current::kCC, sign, sp.systName);
     fits[kNueSurv]  = FitComponent(sp.shifts, sp.preds, Flavors::kNuEToNuE,   Current::kCC, sign, sp.systName);
     fits[kNumuSurv] = FitComponent(sp.shifts, sp.preds, Flavors::kNuMuToNuMu, Current::kCC, sign, sp.systName);
 
     fits[kNC]       = FitComponent(sp.shifts, sp.preds, Flavors::kAll, Current::kNC, sign, sp.systName);
 
     fits[kOther] = FitComponent(sp.shifts, sp.preds, Flavors::kNuEToNuMu | Flavors::kAllNuTau, Current::kCC, sign, sp.systName);
   }

std::unique_ptr< PredictionSystNumu2017 > ana::PredictionSystNumu2017::LoadFrom	(	TDirectory *	dir,
		const std::string &	name
	)

static

Definition at line 208 of file PredictionSystNumu2017.cxx.

References ana::assert(), dir, ana::kNumu2017NormSystBkg, ana::kNumu2017NormSystBoth, ana::kNumu2017NormSystSig, ana::PredictionInterp::LoadFromBody(), runNovaSAM::ret, and getGoodRuns4SAM::tag.

   {
     dir = dir->GetDirectory(name.c_str()); // switch to subdir
     assert(dir);
 
     TObjString* tag = (TObjString*)dir->Get("type");
     assert(tag);
 
     assert(tag->GetString() == "PredictionSystNumu2017");
     delete tag;
 
     std::unique_ptr<PredictionSystNumu2017> ret(new PredictionSystNumu2017);
     LoadFromBody(dir, ret.get(), {&kNumu2017NormSystBoth, &kNumu2017NormSystSig, &kNumu2017NormSystBkg});
 
     // Just enable all norm systs unconditionally. Way too hard to do otherwise
     ret->AddNormSysts({&kNumu2017NormSystBoth, &kNumu2017NormSystSig, &kNumu2017NormSystBkg});
 
     delete dir;
 
     return ret;
   }

void ana::PredictionInterp::LoadFromBody	(	TDirectory *	dir,
		PredictionInterp *	ret,
		std::vector< const ISyst * >	veto = `{}`
	)

staticinherited

Definition at line 663 of file PredictionInterp.cxx.

References dir, genie::utils::style::Format(), ana::PredictionInterp::fOscOrigin, ana::PredictionInterp::fPredNom, ana::PredictionInterp::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(), ana::PredictionInterp::LoadFrom(), ana::PredictionSystNue2017::LoadFrom(), LoadFrom(), and ana::PredictionInterp::PredictionInterp().

   {
     ret->fPredNom = ana::LoadFrom<IPrediction>(dir, "pred_nom");
 
     TH1* hSystNames = (TH1*)dir->Get("syst_names");
     if(hSystNames){
       for(int systIdx = 0; systIdx < hSystNames->GetNbinsX(); ++systIdx){
         ShiftedPreds sp;
         sp.systName = hSystNames->GetXaxis()->GetBinLabel(systIdx+1);
 
         const ISyst* syst = Registry<ISyst>::ShortNameToPtr(sp.systName);
 
         if(std::find(veto.begin(), veto.end(), syst) != veto.end()) continue;
 
         for(int shift = -3; shift <= +3; ++shift){
           const std::string subname = TString::Format("pred_%s_%+d", sp.systName.c_str(), shift).Data();
           TDirectory *preddir = dir->GetDirectory(subname.c_str());
           if(!preddir) continue; // Can happen for genie systs
           delete preddir;
 
           sp.shifts.push_back(shift);
           sp.preds.emplace_back(ana::LoadFrom<IPrediction>(dir, subname).release());
         } // end for shift
 
         ret->fPreds.emplace(syst, std::move(sp));
       } // end for systIdx
     } // end if hSystNames
 
     ret->fOscOrigin = ana::LoadFrom<osc::IOscCalc>(dir, "osc_origin").release();
   }

void ana::PredictionInterp::MinimizeMemory ( )

inherited

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

Definition at line 696 of file PredictionInterp.cxx.

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

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

   {
     std::set<IPrediction*> todel;
     for(auto& it: fPreds){
       std::vector<IPrediction*>& preds = it.second.preds;
       for(unsigned int i = 0; i < preds.size(); ++i){
         if(preds[i] != fPredNom.get()){
           todel.insert(preds[i]);
           preds[i] = 0;
         }
       }
     }
 
     for(IPrediction* p: todel) delete p;
 
     // We probably just freed up a lot of memory, but malloc by default hangs
     // on to all of it as cache.
     #ifndef DARWINBUILD
     malloc_trim(0);
     #endif
   }

Spectrum ana::PredictionInterp::Predict ( osc::IOscCalc * calc ) const

overridevirtualinherited

Implements ana::IPrediction.

Definition at line 305 of file PredictionInterp.cxx.

References ana::PredictionInterp::fPredNom.

Referenced by ana::PredictionSystJoint2018::AddNormSyst(), ana::PredictionSystJointDemo::AddNormSyst(), ana::PredictionSyst3Flavor2020::AddNormSyst(), ana::PredictionSystNue2017::AddNormSysts(), AddNormSysts(), drawSystsShiftingNDdata_updatedAna(), GetNuePrediction2017(), GetNuePrediction2018(), GetNuePrediction2019(), ana::GetNuePrediction2020(), getSensitivity(), makeSystTable(), makeSystTable_reduced(), saveContours_oscpar(), saveContours_systs(), syst_plot_test(), test_stanfit_systpulls(), TestPred(), and Toy_analyses().

   {
     return fPredNom->Predict(calc);
   }

Spectrum ana::PredictionInterp::Predict ( osc::IOscCalcStan * calc ) const

overridevirtualinherited

Reimplemented from ana::IPrediction.

Definition at line 311 of file PredictionInterp.cxx.

References ana::PredictionInterp::fPredNom.

   {
     return fPredNom->Predict(calc);
   }

Spectrum ana::PredictionInterp::PredictComponent	(	osc::IOscCalc *	calc,
		Flavors::Flavors_t	flav,
		Current::Current_t	curr,
		Sign::Sign_t	sign
	)		const

overridevirtualinherited

Implements ana::IPrediction.

Definition at line 317 of file PredictionInterp.cxx.

References ana::PredictionInterp::fPredNom.

Referenced by syst_plot_test().

   {
     return fPredNom->PredictComponent(calc, flav, curr, sign);
   }

Spectrum ana::PredictionInterp::PredictComponent	(	osc::IOscCalcStan *	calc,
		Flavors::Flavors_t	flav,
		Current::Current_t	curr,
		Sign::Sign_t	sign
	)		const

overridevirtualinherited

Reimplemented from ana::IPrediction.

Definition at line 326 of file PredictionInterp.cxx.

References ana::PredictionInterp::fPredNom.

   {
     return fPredNom->PredictComponent(calc, flav, curr, sign);
   }

Spectrum ana::PredictionInterp::PredictComponentSyst	(	osc::IOscCalc *	calc,
		const SystShifts &	shift,
		Flavors::Flavors_t	flav,
		Current::Current_t	curr,
		Sign::Sign_t	sign
	)		const

overridevirtualinherited

Reimplemented from ana::IPrediction.

Definition at line 570 of file PredictionInterp.cxx.

References ana::PredictionInterp::_PredictComponentSyst().

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

   {
     return _PredictComponentSyst(calc, shift, flav, curr, sign);
   }

Spectrum ana::PredictionInterp::PredictComponentSyst	(	osc::IOscCalcStan *	calc,
		const SystShifts &	shift,
		Flavors::Flavors_t	flav,
		Current::Current_t	curr,
		Sign::Sign_t	sign
	)		const

overridevirtualinherited

Reimplemented from ana::IPrediction.

Definition at line 580 of file PredictionInterp.cxx.

References ana::PredictionInterp::_PredictComponentSyst().

   {
     return _PredictComponentSyst(calc, shift, flav, curr, sign);
   }

Spectrum ana::PredictionInterp::PredictSyst	(	osc::IOscCalc *	calc,
		const SystShifts &	shift
	)		const

overridevirtualinherited

Reimplemented from ana::IPrediction.

Definition at line 335 of file PredictionInterp.cxx.

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

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

   {
     InitFits();
 
     return PredictComponentSyst(calc, shift,
                                 Flavors::kAll,
                                 Current::kBoth,
                                 Sign::kBoth);
   }

Spectrum ana::PredictionInterp::PredictSyst	(	osc::IOscCalcStan *	calc,
		const SystShifts &	shift
	)		const

overridevirtualinherited

Reimplemented from ana::IPrediction.

Definition at line 347 of file PredictionInterp.cxx.

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

   {
     InitFits();
 
     return PredictComponentSyst(calc,
                                 shift,
                                 Flavors::kAll,
                                 Current::kBoth,
                                 Sign::kBoth);
   }

Spectrum ana::IPrediction::PredictUnoscillated ( ) const

virtualinherited

Reimplemented in ana::PredictionSterile.

Definition at line 56 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().

   {
     // Default implementation
     osc::NoOscillations noosc;
     return Predict(&noosc);
   }

void ana::PredictionSystNumu2017::SaveTo	(	TDirectory *	dir,
		const std::string &	name
	)		const

overridevirtual

Reimplemented from ana::PredictionInterp.

Definition at line 231 of file PredictionSystNumu2017.cxx.

References dir, ana::PredictionInterp::SaveTo(), and tmp.

   {
     PredictionInterp::SaveTo(dir, name);
 
     TDirectory* tmp = gDirectory;
 
     dir = (TDirectory*)dir->GetDirectory(name.c_str());
     dir->cd();
     // Overwrite "IPrediction"
     TObjString("PredictionSystNumu2017").Write("type");
 
     dir->Write();
     delete dir;
 
     tmp->cd();
   }

void ana::PredictionInterp::SetOscSeed ( osc::IOscCalc * oscSeed )

inherited

Definition at line 298 of file PredictionInterp.cxx.

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

                                                        {
     fOscOrigin = oscSeed->Copy();
     for(auto& it: fPreds) it.second.fits.clear();
     InitFits();
   }

template<typename T >

void ana::PredictionInterp::ShiftBins	(	unsigned int	N,
		T *	arr,
		CoeffsType	type,
		bool	nubar,
		const SystShifts &	shift
	)		const

protectedinherited

Helper for ShiftSpectrum.

Definition at line 388 of file PredictionInterp.cxx.

Referenced by ana::PredictionInterp::ShiftSpectrum().

   {
     static_assert(std::is_same_v<T, double> ||
                   std::is_same_v<T, stan::math::var>,
                   "PredictionInterp::ShiftBins() can only be called using doubles or stan::math::vars");
     if(nubar) assert(fSplitBySign);
 
 
     if(!std::is_same_v<T, stan::math::var> && shift.HasAnyStan()){
       std::cout << "PredictionInterp: stan shifts on non-stan spectrum, something is wrong" << std::endl;
       abort();
     }
 
     T corr[N];
     for(unsigned int i = 0; i < N; ++i) corr[i] = 1;
 
     for(auto& it: fPreds){
       const ISyst* syst = it.first;
       const ShiftedPreds& sp = it.second;
 
       T x = shift.GetShift<T>(syst);
 
       // need to actually do the calculation for the autodiff version
       // to make sure the gradient is computed correctly
       if(x == 0 && (!std::is_same_v<T, stan::math::var> || !shift.HasStan(syst))) continue;
 
       int shiftBin = (util::GetValAs<double>(x) - sp.shifts[0])/sp.Stride();
       shiftBin = std::max(0, shiftBin);
       shiftBin = std::min(shiftBin, sp.nCoeffs - 1);
 
       const Coeffs* fits = nubar ? &sp.fitsNubarRemap[type][shiftBin].front() : &sp.fitsRemap[type][shiftBin].front();
 
       x -= sp.shifts[shiftBin];
 
       const T x_cube = util::cube(x);
       const T x_sqr = util::sqr(x);
 
       for(unsigned int n = 0; n < N; ++n){
         // Uncomment to debug crashes in this function
         // assert(type < fits.size());
         // assert(n < sp.fits[type].size());
         // assert(shiftBin < int(sp.fits[type][n].size()));
 
         const Coeffs& f = fits[n];
 
         corr[n] *= f.a*x_cube + f.b*x_sqr + f.c*x + f.d;
      } // end for n
     } // end for syst
 
     for(unsigned int n = 0; n < N; ++n){
       // std::max() doesn't work with stan::math::var
       arr[n] *= (corr[n] > 0.) ? corr[n] : 0.;
     }
 
   }

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

inherited

Helper for PredictComponentSyst.

Definition at line 449 of file PredictionInterp.cxx.

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

Referenced by ana::PredictionInterp::_PredictComponentSyst(), and ana::PredictionInterp::_ShiftedComponent().

   {
     return _ShiftedComponent(calc, hash, shift, flav, curr, sign, type);
   }

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

inherited

Spectrum ana::PredictionInterp::ShiftSpectrum	(	const Spectrum &	s,
		CoeffsType	type,
		bool	nubar,
		const SystShifts &	shift
	)		const

inherited

Definition at line 360 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 ana::PredictionInterp::ShiftBins().

Referenced by ana::PredictionInterp::_ShiftedComponent().

   {
     if(s.HasStan() && (!shift.IsNominal() && !shift.HasAnyStan())){
       std::cout << "PredictionInterp used with Stan oscillations but non-Stan syst shifts. Is that what you expected?" << std::endl;
       std::cout << "  Systs set: " << std::endl;
       for (const auto & syst : shift.ActiveSysts())
         std::cout << "      " << syst->ShortName() << " = " << shift.GetShift(syst) << std::endl;
     }
 
     if(s.HasStan() || shift.HasAnyStan()){
       // Need the second case for the NC component
       Eigen::ArrayXstan vec = s.HasStan() ? s.GetEigenStan(s.POT()) : Eigen::ArrayXstan(s.GetEigen(s.POT()));
 
       ShiftBins(vec.size(), vec.data(), type, nubar, shift);
       return Spectrum(std::move(vec), HistAxis(s.GetLabels(), s.GetBinnings()), s.POT(), s.Livetime());
     }
     else{
       Eigen::ArrayXd vec = s.GetEigen(s.POT());
       ShiftBins(vec.size(), vec.data(), type, nubar, shift);
       return Spectrum(std::move(vec), HistAxis(s.GetLabels(), s.GetBinnings()), s.POT(), s.Livetime());
     }
   }

bool ana::PredictionInterp::SplitBySign ( ) const

inlineinherited

Definition at line 116 of file PredictionInterp.h.

References ana::PredictionInterp::fSplitBySign.

116 {return fSplitBySign;}

ana::PredictionInterp::fSplitBySign

bool fSplitBySign

Definition: PredictionInterp.h:225

Member Data Documentation

Spectrum ana::PredictionInterp::fBinning

mutableprotectedinherited

Dummy spectrum to provide binning.

Definition at line 205 of file PredictionInterp.h.

Referenced by ana::PredictionInterp::_PredictComponentSyst(), and ana::PredictionInterp::InitFits().

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

mutableprotectedinherited

Definition at line 223 of file PredictionInterp.h.

Referenced by ana::PredictionInterp::_ShiftedComponent().

osc::IOscCalc* ana::PredictionInterp::fOscOrigin

protectedinherited

The oscillation values we assume when evaluating the coefficients.

Definition at line 203 of file PredictionInterp.h.

Referenced by ana::PredictionSystJoint2018::AddNormSyst(), ana::PredictionSystJointDemo::AddNormSyst(), ana::PredictionSyst3Flavor2020::AddNormSyst(), ana::PredictionSystNue2017::AddNormSysts(), AddNormSysts(), ana::PredictionInterp::FitComponent(), ana::PredictionInterp::InitFits(), ana::PredictionInterp::LoadFromBody(), ana::PredictionSyst3Flavor2020::PredictionSyst3Flavor2020(), ana::PredictionSystJoint2018::PredictionSystJoint2018(), ana::PredictionSystJointDemo::PredictionSystJointDemo(), ana::PredictionSystNue2017::PredictionSystNue2017(), PredictionSystNumu2017(), ana::PredictionInterp::SaveTo(), and ana::PredictionInterp::SetOscSeed().

std::unique_ptr<IPrediction> ana::PredictionInterp::fPredNom

protectedinherited

The nominal prediction.

Definition at line 174 of file PredictionInterp.h.

std::unordered_map<const ISyst*, ShiftedPreds> ana::PredictionInterp::fPreds

mutableprotectedinherited

Definition at line 199 of file PredictionInterp.h.

bool ana::PredictionInterp::fSplitBySign

protectedinherited

Definition at line 225 of file PredictionInterp.h.

Referenced by ana::PredictionInterp::_ShiftedComponent(), ana::PredictionSystJoint2018::AddNormSyst(), ana::PredictionSystJointDemo::AddNormSyst(), ana::PredictionSyst3Flavor2020::AddNormSyst(), ana::PredictionInterp::InitFits(), ana::PredictionSyst3Flavor2020::PredictionSyst3Flavor2020(), ana::PredictionSystJoint2018::PredictionSystJoint2018(), ana::PredictionSystJointDemo::PredictionSystJointDemo(), ana::PredictionInterp::SaveTo(), ana::PredictionInterp::ShiftBins(), and ana::PredictionInterp::SplitBySign().

std::unordered_map<const ISyst*, ShiftedPreds> ana::PredictionInterp::fSumPreds

mutableprotectedinherited

Definition at line 200 of file PredictionInterp.h.

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

/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-09-26/3FlavorAna/Prediction/PredictionSystNumu2017.h
/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-09-26/3FlavorAna/Prediction/PredictionSystNumu2017.cxx

Public Types

Public Member Functions

Static Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation