ana::PredictionInterp Class Reference

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

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-05-06/CAFAna/Prediction/PredictionInterp.h"

Inheritance diagram for ana::PredictionInterp:

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< PredictionInterp >	LoadFrom (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< 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

Implements systematic errors by interpolation between shifted templates.

Definition at line 23 of file PredictionInterp.h.

Member Enumeration Documentation

enum ana::PredictionInterp::CoeffsType

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

Enumerator
kCombineSigns
kSplitBySign

Definition at line 26 of file PredictionInterp.h.

                {
      kCombineSigns, kSplitBySign
    };

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

systs	What systematics we should be capable of interpolating
osc	The oscillation point to expand around
predGen	Construct an IPrediction from the following information.
loaders	The loaders to be passed on to the underlying prediction
shiftMC	Underlying shift. Use with care. Mostly for PredictionNumuFAHadE. Should not contain any of of systs
mode	In 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 58 of file PredictionInterp.cxx.

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

    : fOscOrigin(osc ? osc->Copy() : 0),
      fBinning(Spectrum::Uninitialized()),
      fSplitBySign(mode == kSplitBySign)
  {
    for(const ISyst* syst: systs){
      ShiftedPreds sp;
      sp.systName = syst->ShortName();

      // Genie reweight shifts aren't able to provide +/-3 sigma
      if(syst->IsGenieReweight())
        sp.shifts = {-2, -1, 0, +1, +2};
      else
        sp.shifts = {-3, -2, -1, 0, +1, +2, +3};

      for(int sigma: sp.shifts){
        SystShifts shiftHere = shiftMC;
        shiftHere.SetShift(syst, sigma);
        sp.preds.emplace_back(predGen.Generate(loaders, shiftHere).release());
      }

      fPreds.emplace(syst, std::move(sp));
    } // end for syst

    fPredNom = predGen.Generate(loaders, shiftMC);
  }

ana::PredictionInterp::~PredictionInterp ( )

virtual

Definition at line 91 of file PredictionInterp.cxx.

References FitRatios().

  {
  }

ana::PredictionInterp::PredictionInterp ( )

inline

Definition at line 123 of file PredictionInterp.h.

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

: fOscOrigin(nullptr), fBinning(Spectrum::Uninitialized()), fSplitBySign(false) {}

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 526 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().

  {
    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

protected

Templated helper for ShiftedComponent.

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

  {
    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);
  }

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().

  {
    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 115 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 120 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 110 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

Definition at line 717 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, it, nbins, ana::IPrediction::PredictComponent(), PredictComponentSyst(), PandAna.reco_validation.prod5_pid_validation::ratio(), ana::ISyst::ShortName(), std::sqrt(), ss, ana::Spectrum::ToTH1(), and submit_syst::x.

Referenced by 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(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

Definition at line 821 of file PredictionInterp.cxx.

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

Referenced by 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

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

  {
    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

Definition at line 856 of file PredictionInterp.cxx.

References ana::assert(), DebugPlotColz(), genie::utils::style::Format(), fPreds, 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

Find coefficients describing the ratios from this component.

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

  {
    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

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, 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().

  {
    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

Definition at line 240 of file PredictionInterp.cxx.

References fPreds.

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

    return allsysts;
  }

void ana::PredictionInterp::InitFits ( ) const

protected

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

  {
    // 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

protected

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

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

static

Definition at line 640 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(), PlotNusSensAna01(), and SetPrediction().

  {
    dir = dir->GetDirectory(name.c_str()); // switch to subdir
    assert(dir);

    TObjString* tag = (TObjString*)dir->Get("type");
    assert(tag);
    assert(tag->GetString() == "PredictionInterp");

    std::unique_ptr<PredictionInterp> ret(new PredictionInterp);

    LoadFromBody(dir, ret.get());

    TObjString* split_sign = (TObjString*)dir->Get("split_sign");
    // Can be missing from old files
    ret->fSplitBySign = (split_sign && split_sign->String() == "yes");

    delete dir;

    return ret;
  }

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

static

Definition at line 663 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(), string, and ana::PredictionInterp::ShiftedPreds::systName.

Referenced by ana::PredictionSystJoint2018::LoadFrom(), ana::PredictionSystJointDemo::LoadFrom(), LoadFrom(), ana::PredictionSystNue2017::LoadFrom(), ana::PredictionSystNumu2017::LoadFrom(), and 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

(

)

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

Definition at line 696 of file PredictionInterp.cxx.

References fPredNom, 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.
    malloc_trim(0);
  }

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

overridevirtual

Implements ana::IPrediction.

Definition at line 305 of file PredictionInterp.cxx.

References fPredNom.

Referenced by ana::PredictionSystJoint2018::AddNormSyst(), ana::PredictionSystJointDemo::AddNormSyst(), ana::PredictionSyst3Flavor2020::AddNormSyst(), ana::PredictionSystNue2017::AddNormSysts(), ana::PredictionSystNumu2017::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

overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 311 of file PredictionInterp.cxx.

References 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

overridevirtual

Implements ana::IPrediction.

Definition at line 317 of file PredictionInterp.cxx.

References 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

overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 326 of file PredictionInterp.cxx.

References 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

overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 570 of file PredictionInterp.cxx.

References _PredictComponentSyst().

Referenced by check_predinterp(), check_predinterp_numu(), DebugPlot(), DebugPlotColz(), and 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

overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 580 of file PredictionInterp.cxx.

References _PredictComponentSyst().

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

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

overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 335 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_nux_tester_2020(), 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

overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 347 of file PredictionInterp.cxx.

References InitFits(), ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth, and 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 33 of file IPrediction.cxx.

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

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

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

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 589 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().

  {
    InitFits();

    TDirectory* tmp = gDirectory;

    dir = dir->mkdir(name.c_str()); // switch to subdir
    dir->cd();

    TObjString("PredictionInterp").Write("type");

    fPredNom->SaveTo(dir, "pred_nom");

    for(auto& it: fPreds){
      const ShiftedPreds& sp = it.second;

      for(unsigned int i = 0; i < sp.shifts.size(); ++i){
        if(!sp.preds[i]){
          std::cout << "Can't save a PredictionInterp after MinimizeMemory()" << std::endl;
          abort();
        }

        sp.preds[i]->SaveTo(dir, TString::Format("pred_%s_%+d",
                                                 sp.systName.c_str(),
                                                 int(sp.shifts[i])).Data());
      } // end for i
    } // end for it

    ana::SaveTo(*fOscOrigin, dir, "osc_origin");

    if(!fPreds.empty()){
      TH1F hSystNames("syst_names", ";Syst names", fPreds.size(), 0, fPreds.size());
      int binIdx = 1;
      for(auto& it: fPreds){
        hSystNames.GetXaxis()->SetBinLabel(binIdx++, it.second.systName.c_str());
      }
      dir->cd();
      hSystNames.Write("syst_names");
    }

    TObjString split_sign = fSplitBySign ? "yes" : "no";
    dir->cd();
    split_sign.Write("split_sign");

    dir->Write();
    delete dir;

    tmp->cd();
  }

void ana::PredictionInterp::SetOscSeed

(

osc::IOscCalc *

oscSeed

)

Definition at line 298 of file PredictionInterp.cxx.

References osc::_IOscCalc< T >::Copy(), fOscOrigin, fPreds, 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

protected

Helper for ShiftSpectrum.

Definition at line 388 of file PredictionInterp.cxx.

References ana::PredictionInterp::Coeffs::a, ana::assert(), ana::PredictionInterp::Coeffs::b, ana::PredictionInterp::Coeffs::c, dg::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(), ana::PredictionInterp::ShiftedPreds::nCoeffs, ana::PredictionInterp::ShiftedPreds::shifts, util::sqr(), ana::PredictionInterp::ShiftedPreds::Stride(), T, and submit_syst::x.

Referenced by 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

Helper for PredictComponentSyst.

Definition at line 449 of file PredictionInterp.cxx.

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

Referenced by _PredictComponentSyst(), and _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

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

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 ShiftBins().

Referenced by _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

inline

Definition at line 116 of file PredictionInterp.h.

References fSplitBySign.

{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

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(), ana::PredictionSystNumu2017::AddNormSysts(), FitComponent(), InitFits(), LoadFromBody(), ana::PredictionSyst3Flavor2020::PredictionSyst3Flavor2020(), ana::PredictionSystJoint2018::PredictionSystJoint2018(), ana::PredictionSystJointDemo::PredictionSystJointDemo(), ana::PredictionSystNue2017::PredictionSystNue2017(), ana::PredictionSystNumu2017::PredictionSystNumu2017(), SaveTo(), and SetOscSeed().

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

protected

The nominal prediction.

Definition at line 174 of file PredictionInterp.h.

Referenced by _ShiftedComponent(), ana::PredictionSystJoint2018::AddNormSyst(), ana::PredictionSystJointDemo::AddNormSyst(), ana::PredictionSyst3Flavor2020::AddNormSyst(), ana::PredictionSystNue2017::AddNormSysts(), ana::PredictionSystNumu2017::AddNormSysts(), DebugPlot(), DebugPlotColz(), ana::PredictionSystJoint2018::GetNominalPrediction(), ana::PredictionSystJointDemo::GetNominalPrediction(), InitFits(), LoadFromBody(), MinimizeMemory(), Predict(), PredictComponent(), PredictionInterp(), ana::PredictionSyst3Flavor2020::PredictionSyst3Flavor2020(), ana::PredictionSystJoint2018::PredictionSystJoint2018(), ana::PredictionSystJointDemo::PredictionSystJointDemo(), ana::PredictionSystNue2017::PredictionSystNue2017(), ana::PredictionSystNumu2017::PredictionSystNumu2017(), and SaveTo().

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

mutableprotected

Definition at line 199 of file PredictionInterp.h.

Referenced by _PredictComponentSyst(), ana::PredictionSystJoint2018::AddNormSyst(), ana::PredictionSystJointDemo::AddNormSyst(), ana::PredictionSyst3Flavor2020::AddNormSyst(), ana::PredictionSystNue2017::AddNormSysts(), ana::PredictionSystNumu2017::AddNormSysts(), DebugPlot(), DebugPlots(), DebugPlotsColz(), GetAllSysts(), ana::PredictionSystJoint2018::GetShiftedPrediction(), ana::PredictionSystJointDemo::GetShiftedPrediction(), InitFits(), LoadFromBody(), MinimizeMemory(), PredictionInterp(), ana::PredictionSyst3Flavor2020::PredictionSyst3Flavor2020(), ana::PredictionSystJoint2018::PredictionSystJoint2018(), ana::PredictionSystJointDemo::PredictionSystJointDemo(), ana::PredictionSystNue2017::PredictionSystNue2017(), ana::PredictionSystNumu2017::PredictionSystNumu2017(), SaveTo(), SetOscSeed(), and ShiftBins().

bool ana::PredictionInterp::fSplitBySign

protected

Definition at line 225 of file PredictionInterp.h.

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

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:

• /cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-05-06/CAFAna/Prediction/PredictionInterp.h
• /cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-05-06/CAFAna/Prediction/PredictionInterp.cxx