ana::PredictionInterp Class Reference

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

#include <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 ()
virtual Spectrum	Predict (osc::IOscCalc *calc) const override
virtual Spectrum	PredictSyst (osc::IOscCalc *calc, const SystShifts &shift) const override
virtual Spectrum	PredictComponent (osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
virtual Spectrum	PredictComponentSyst (osc::IOscCalc *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
virtual void	SaveTo (TDirectory *dir) 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< std::unique_ptr< TH1 >> &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 &shortName) const
	Find coefficients describing the ratios from this component. More...
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...
Public Member Functions inherited from ana::IPrediction
virtual	~IPrediction ()
virtual Spectrum	PredictUnoscillated () const
virtual OscillatableSpectrum	ComponentCC (int from, int to) const

Static Public Member Functions
static std::unique_ptr < PredictionInterp >	LoadFrom (TDirectory *dir)

Public Attributes
std::unique_ptr< IPrediction >	fPredNom
	The nominal prediction. More...

Protected Member Functions
void	InitFits () const
void	InitFitsHelper (ShiftedPreds &sp, std::vector< std::vector< std::vector< Coeffs >>> &fits, Sign::Sign_t sign) const

Protected Attributes
std::unordered_map< const ISyst *, ShiftedPreds >	fPreds
osc::IOscCalc *	fOscOrigin
	The oscillation values we assume when evaluating the coefficients. More...
Spectrum	fBinning
	Dummy spectrum to provide binning. More...
std::map< Key_t, Val_t >	fNomCache
bool	fSplitBySign

Detailed Description

Implements systematic errors by interpolation between shifted templates.

Definition at line 22 of file PredictionInterp.h.

Member Enumeration Documentation

enum ana::PredictionInterp::CoeffsType

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

Definition at line 101 of file PredictionInterp.h.

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

enum ana::PredictionInterp::EMode_t

Enumerator
kCombineSigns
kSplitBySign

Definition at line 25 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 29 of file PredictionInterp.cxx.

    : fOscOrigin(osc ? osc->Copy() : 0),
      fBinning(0, {}, {}, 0, 0),

ana::PredictionInterp::~PredictionInterp ( )

virtual

Definition at line 60 of file PredictionInterp.cxx.

  {
    //    for(auto it: fPreds) for(IPrediction* p: it.second.preds) delete p;
    //    delete fOscOrigin;

    // It isn't really a unique ptr when we use PredictionInterpTemplates
    fPredNom.release();
  }

ana::PredictionInterp::PredictionInterp ( )

inline

Definition at line 107 of file PredictionInterp.h.

: fBinning(0, {}, {}, 0, 0) {}

Member Function Documentation

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 556 of file PredictionInterp.cxx.

  {
    DontAddDirectory guard;

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

    std::vector<TGraph*> curves(nbins);
    std::vector<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];
    }
    assert(pNom);
    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


    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(UniqueName().c_str(),
                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 677 of file PredictionInterp.cxx.

  {
    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 644 of file PredictionInterp.cxx.

  {
    const bool save = !savePattern.empty();
    const bool multiFile = savePattern.find("%s") != std::string::npos;

    if(save && !multiFile){
      new TCanvas;
      gPad->Print((savePattern+"[").c_str());
    }

    for(auto& it: fPreds){
      DebugPlot(it.first, calc, flav, curr, sign);

      if(save){
        if(multiFile){
          gPad->Print(TString::Format(savePattern.c_str(), it.second.systName.c_str()).Data());
        }
        else{
          gPad->Print(savePattern.c_str());
        }
      }
    } // end for it

    if(save && !multiFile){
      gPad->Print((savePattern+"]").c_str());
    }
  }

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 712 of file PredictionInterp.cxx.

  {
    const bool save = !savePattern.empty();
    const bool multiFile = savePattern.find("%s") != std::string::npos;

    if(save && !multiFile){
      new TCanvas;
      gPad->Print((savePattern+"[").c_str());
    }

    for(auto it: fPreds){
      new TCanvas;
      DebugPlotColz(it.first, calc, flav, curr, sign);

      if(save){
        if(multiFile){
          gPad->Print(TString::Format(savePattern.c_str(), it.second.systName.c_str()).Data());
        }
        else{
          gPad->Print(savePattern.c_str());
        }
      }
    } // end for it

    if(save && !multiFile){
      gPad->Print((savePattern+"]").c_str());
    }
  }

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 &	shortName
	)		const

Find coefficients describing the ratios from this component.

Definition at line 167 of file PredictionInterp.cxx.

  {
    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<std::unique_ptr<TH1>> ratios;
    for(auto& p: preds){
      ratios.emplace_back(Ratio(p->PredictComponent(fOscOrigin,
                                                    flav, curr, sign),
                                nom).ToTH1());

      // Check none of the ratio values is utterly crazy
      std::unique_ptr<TH1>& r = ratios.back();
      for(int i = 0; i < r->GetNbinsX()+2; ++i){
        const double y = r->GetBinContent(i);
        if(y > 500){
          std::cout << "PredictionInterp: WARNING, ratio in bin "
                    << i << " is " << y
                    << " for '" << shortName << "'. Ignoring." << std::endl;
          r->SetBinContent(i, 1);
        }
      }
    }

    return FitRatios(shifts, ratios);
  }

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

Find coefficients describing this set of shifts.

Definition at line 71 of file PredictionInterp.cxx.

  {
    assert(shifts.size() == ratios.size());

    std::vector<std::vector<Coeffs>> ret;

    const int binMax = ratios[0]->GetNbinsX();

    // Linear interpolation
    if(ratios.size() == 2){
      for(int binIdx = 0; binIdx < binMax+2; ++binIdx){
        ret.push_back({});
        const double y1 = ratios[0]->GetBinContent(binIdx);
        const double y2 = ratios[1]->GetBinContent(binIdx);
        ret.back().emplace_back(0, 0, y2-y1, y1);
      }
      return ret;
    }

    for(int binIdx = 0; binIdx < binMax+2; ++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.

      {
        const double y1 = ratios[0]->GetBinContent(binIdx);
        const double y2 = ratios[1]->GetBinContent(binIdx);
        const double y3 = ratios[2]->GetBinContent(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]->GetBinContent(binIdx);
        const double y1 = ratios[shiftIdx  ]->GetBinContent(binIdx);
        const double y2 = ratios[shiftIdx+1]->GetBinContent(binIdx);
        const double y3 = ratios[shiftIdx+2]->GetBinContent(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  ]->GetBinContent(binIdx);
        const double y1 = ratios[N+1]->GetBinContent(binIdx);
        const double y2 = ratios[N+2]->GetBinContent(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;
  }

void ana::PredictionInterp::InitFits ( ) const

protected

Definition at line 225 of file PredictionInterp.cxx.

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

    // 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 209 of file PredictionInterp.cxx.

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

static

Definition at line 479 of file PredictionInterp.cxx.

  {
    TObjString* tag = (TObjString*)dir->Get("type");
    assert(tag);
    assert(tag->GetString() == "PredictionInterp" ||
           tag->GetString() == "PredictionInterp2"); // backwards compatibility

    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");

    return ret;
  }

void ana::PredictionInterp::MinimizeMemory

(

)

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

Definition at line 531 of file PredictionInterp.cxx.

  {
    InitFits();

    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

overridevirtual

Implements ana::IPrediction.

Definition at line 259 of file PredictionInterp.cxx.

  {
    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 265 of file PredictionInterp.cxx.

  {
    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 391 of file PredictionInterp.cxx.

  {
    InitFits();

    Spectrum& ret = fBinning;
    ret.Clear();

    if(ret.POT()==0) ret.OverridePOT(1e24);

    // 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;
  }

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

overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 274 of file PredictionInterp.cxx.

  {
    InitFits();

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

void ana::PredictionInterp::SaveTo ( TDirectory *  dir ) const

overridevirtual

Reimplemented from ana::IPrediction.

Definition at line 436 of file PredictionInterp.cxx.

  {
    TDirectory* tmp = gDirectory;

    dir->cd();
    TObjString("PredictionInterp").Write("type");


    fPredNom->SaveTo(dir->mkdir("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->mkdir(TString::Format("pred_%s_%+d",
                                                       sp.systName.c_str(),
                                                       int(sp.shifts[i])).Data()));
      } // end for i
    } // end for it

    ana::SaveTo(*fOscOrigin, dir->mkdir("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());
      }
      hSystNames.Write("syst_names");
    }

    TObjString split_sign = fSplitBySign ? "yes" : "no";
    split_sign.Write("split_sign");

    tmp->cd();
  }

void ana::PredictionInterp::SetOscSeed

(

osc::IOscCalc *

oscSeed

)

Definition at line 253 of file PredictionInterp.cxx.

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

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 345 of file PredictionInterp.cxx.

  {
    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 calcs won't hash themselves.
    const bool canCache = (hash != 0);

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

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

Definition at line 287 of file PredictionInterp.cxx.

  {
    // Save time by not shifting a spectrum that is all zeros anyway
    if(s.POT() == 0 && s.Livetime() == 0) return s;

    if(nubar) assert(fSplitBySign);

    InitFits();

    // TODO histogram operations could be too slow
    TH1D* h = s.ToTH1(s.POT());

    const unsigned int N = h->GetNbinsX()+2;
    std::vector<double> corr(N,1.0);

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

      auto& fits = nubar ? sp.fitsNubar : sp.fits;

      double x = shift.GetShift(syst);

      if(x == 0) continue;

      int shiftBin = (x - sp.shifts[0])/sp.Stride();
      shiftBin = std::max(0, shiftBin);
      shiftBin = std::min(shiftBin, sp.nCoeffs-1);

      x -= sp.shifts[shiftBin];

      const double x_cube = util::cube(x);
      const double 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[type][n][shiftBin];

        corr[n] *= f.a*x_cube + f.b*x_sqr + f.c*x + f.d;
      } // end for n
    } // end for syst

    double* arr = h->GetArray();
    for(unsigned int n = 0; n < N; ++n){
      arr[n] *= std::max(corr[n], 0.);
    }

    return Spectrum(std::unique_ptr<TH1D>(h), s.GetLabels(), s.GetBinnings(), s.POT(), s.Livetime());
  }

bool ana::PredictionInterp::SplitBySign ( ) const

inline

Definition at line 100 of file PredictionInterp.h.

{return fSplitBySign;}

Member Data Documentation

Spectrum ana::PredictionInterp::fBinning

mutableprotected

Dummy spectrum to provide binning.

Definition at line 170 of file PredictionInterp.h.

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

mutableprotected

Definition at line 188 of file PredictionInterp.h.

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

protected

The oscillation values we assume when evaluating the coefficients.

Definition at line 168 of file PredictionInterp.h.

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

The nominal prediction.

Definition at line 146 of file PredictionInterp.h.

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

mutableprotected

Definition at line 165 of file PredictionInterp.h.

bool ana::PredictionInterp::fSplitBySign

protected

Definition at line 190 of file PredictionInterp.h.

---

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

• PredictionInterp.h
• PredictionInterp.cxx