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

Compare a data spectrum to MC using a covariance matrix in the chi^2. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-08-14/CAFAna/Experiment/CovMxExperiment.h"

Inheritance diagram for ana::CovMxExperiment:
ana::SingleSampleExperiment ana::IExperiment

Classes

struct  Component
 Helper class used to specify components that are separately treated with their own covariance matrices. More...
 

Public Member Functions

 CovMxExperiment (const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, const std::map< Component, TMatrixD * > &covmx, double cosmicScaleError=0)
 
 CovMxExperiment (const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, TMatrixD *const covmx, double cosmicScaleError=0)
 Simplified version for use with a single covariance matrix. More...
 
 CovMxExperiment (const IPrediction *pred, const Spectrum &data, const std::map< Component, TMatrixD * > &covmx)
 No-cosmics version (see first constructor for explanation of parameters) More...
 
 CovMxExperiment (const IPrediction *pred, const Spectrum &data, TMatrixD *const covmx)
 Single-matrix variant of no-cosmics constructor. More...
 
virtual double ChiSq (osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
 
 CovMxExperiment (const SingleSampleExperiment &)=delete
 
CovMxExperimentoperator= (const SingleSampleExperiment &)=delete
 
virtual void Derivative (osc::IOscCalc *calc, const SystShifts &shift, std::unordered_map< const ISyst *, double > &dch) const override
 
virtual stan::math::var LogLikelihood (osc::_IOscCalcAdjustable< stan::math::var > *osc, const SystShifts &syst=SystShifts::Nominal()) const override
 
virtual stan::math::var LogLikelihood (osc::IOscCalcAdjustableStan *osc, const SystShifts &syst=kNoShift) const
 
virtual void SaveTo (TDirectory *dir, const std::string &name) const override
 

Static Public Member Functions

static std::unique_ptr< SingleSampleExperimentLoadFrom (TDirectory *dir, const std::string &name)
 

Protected Member Functions

void MakeInvCovMxs (const std::map< Component, TMatrixD * > &covmx)
 
TH1D * PredHistIncCosmics (osc::IOscCalc *calc, const SystShifts &syst) const
 

Protected Attributes

std::map< Component, std::unique_ptr< TMatrixD > > fCovMxInv
 
const IPredictionfMC
 
Spectrum fData
 
TH1D * fCosmic
 
double fCosmicScaleError
 
bool fPoissonError
 
std::vector< std::pair< double, double > > fBinErrors
 

Detailed Description

Compare a data spectrum to MC using a covariance matrix in the chi^2.

Definition at line 22 of file CovMxExperiment.h.

Constructor & Destructor Documentation

ana::CovMxExperiment::CovMxExperiment ( const IPrediction pred,
const Spectrum data,
const Spectrum cosmic,
const std::map< Component, TMatrixD * > &  covmx,
double  cosmicScaleError = 0 
)
inline
Parameters
predSource of oscillated MC beam predictions
dataData spectrum to compare to
cosmicCosmic ray background component
covmxBin-to-bin covariance matrix for each component you want to separately treat. See below for alternate constructor if you have just a single covariance matrix (or pass a map with the single key Component(Flavors::kAll,Current::kBoth,Sign::kBoth)).
cosmicScaleErrorfractional uncertainty on cosmic normalization

Definition at line 49 of file CovMxExperiment.h.

References MakeInvCovMxs().

Referenced by CovMxExperiment().

54  : SingleSampleExperiment(pred, data, cosmic, cosmicScaleError)
55  {
56  this->MakeInvCovMxs(covmx);
57  };
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
const XML_Char const XML_Char * data
Definition: expat.h:268
Sum up livetimes from individual cosmic triggers.
void MakeInvCovMxs(const std::map< Component, TMatrixD * > &covmx)
ana::CovMxExperiment::CovMxExperiment ( const IPrediction pred,
const Spectrum data,
const Spectrum cosmic,
TMatrixD *const  covmx,
double  cosmicScaleError = 0 
)
inline

Simplified version for use with a single covariance matrix.

Definition at line 60 of file CovMxExperiment.h.

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

65  : SingleSampleExperiment(pred, data, cosmic, cosmicScaleError)
66  {
68  };
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
const XML_Char const XML_Char * data
Definition: expat.h:268
Interactions of both types.
Definition: IPrediction.h:42
Sum up livetimes from individual cosmic triggers.
void MakeInvCovMxs(const std::map< Component, TMatrixD * > &covmx)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
ana::CovMxExperiment::CovMxExperiment ( const IPrediction pred,
const Spectrum data,
const std::map< Component, TMatrixD * > &  covmx 
)
inline

No-cosmics version (see first constructor for explanation of parameters)

Definition at line 71 of file CovMxExperiment.h.

References MakeInvCovMxs().

75  {
76  this->MakeInvCovMxs(covmx);
77  };
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
const XML_Char const XML_Char * data
Definition: expat.h:268
void MakeInvCovMxs(const std::map< Component, TMatrixD * > &covmx)
ana::CovMxExperiment::CovMxExperiment ( const IPrediction pred,
const Spectrum data,
TMatrixD *const  covmx 
)
inline

Single-matrix variant of no-cosmics constructor.

Definition at line 80 of file CovMxExperiment.h.

References ChiSq(), CovMxExperiment(), ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth, MakeInvCovMxs(), ana::SystShifts::Nominal(), and operator=().

84  {
86  };
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
const XML_Char const XML_Char * data
Definition: expat.h:268
Interactions of both types.
Definition: IPrediction.h:42
void MakeInvCovMxs(const std::map< Component, TMatrixD * > &covmx)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
ana::CovMxExperiment::CovMxExperiment ( const SingleSampleExperiment )
delete

Member Function Documentation

double ana::CovMxExperiment::ChiSq ( osc::IOscCalcAdjustable osc,
const SystShifts syst = SystShifts::Nominal() 
) const
overridevirtual

Reimplemented from ana::SingleSampleExperiment.

Definition at line 23 of file CovMxExperiment.cxx.

References bin, chi2(), ana::Chi2CovMx(), ana::HistCache::Delete(), release_diff::diff, e, stan::math::fabs(), ana::SingleSampleExperiment::fCosmic, ana::SingleSampleExperiment::fCosmicScaleError, fCovMxInv, ana::SingleSampleExperiment::fData, ana::SingleSampleExperiment::fMC, ana::SystShifts::GetShift(), ana::SystShifts::IsNominal(), ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth, ana::kCosmicBkgScaleSyst, pot, ana::Spectrum::POT(), ana::IPrediction::Predict(), ana::IPrediction::PredictComponent(), scale, and ana::Spectrum::ToTH1().

Referenced by CovMxExperiment().

24  {
25  // the covariance matrix replaces the systematic shift approach.
26  if (!syst.IsNominal())
27  throw std::runtime_error("ExperimentWithCovMx::ChiSq(): Systematic shifts cannot be specified when fitting against covariance matrices");
28 
29  const double pot = fData.POT();
30  TH1D* hdata = fData.ToTH1(pot);
31 
32  std::map<Component, TH1D*> predictions;
33  for (const auto & mxPair : this->fCovMxInv)
34  {
35  const auto & comp = mxPair.first;
36  if (comp.flavor == Flavors::kAll && comp.current == Current::kBoth && comp.sign == Sign::kBoth)
37  {
38  // if you specified one matrix as applying to the whole prediction,
39  // why did you give me more of them??
40  if (this->fCovMxInv.size() > 1)
41  throw std::runtime_error("CovMxExperiment::ChiSq(): one component corresponds to the whole prediction but you also provided others?!");
42 
43  predictions.emplace(comp, fMC->Predict(calc).ToTH1(pot));
44  }
45  else
46  predictions.emplace(mxPair.first, fMC->PredictComponent(calc, comp.flavor, comp.current, comp.sign).ToTH1(pot));
47  }
48 
49  // will want the total prediction in a moment,
50  // but since we assume the components add to the total
51  // in the following calculations, might as well ensure it here.
52  TH1D * totalPred = fMC->Predict(calc).ToTH1(pot);
53  TH1D diff(*totalPred);
54  for (const auto & predPair : predictions)
55  diff.Add(predPair.second, -1);
56  for (int bin = 0; bin <= totalPred->GetNbinsX() + 1; bin++)
57  {
58  if ( std::fabs(diff.GetBinContent(bin)) > 1e-4 )
59  throw std::runtime_error("CovMxExperiment::ChiSq(): components don't sum to total prediction!");
60  }
61 
62  // add the cosmics into the prediction.
63  if(fCosmic)
64  {
65  if(fCosmicScaleError != 0)
66  {
67  const double scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
68  totalPred->Add(fCosmic, scale);
69  }
70  else
71  totalPred->Add(fCosmic);
72 
73  }
74 
75  // first, compute the statistical uncertainty covariance matrix
76  // (which is diagonal). we'll need this repeatedly below.
77  // use the *total* prediction to get the sigma^2 here because
78  // that's the expected fluctuation size
79  TMatrixD stat_mx(totalPred->GetNbinsX(), totalPred->GetNbinsX());
80  for (int binNum = 1; binNum <= totalPred->GetNbinsX(); binNum++)
81  stat_mx[binNum-1][binNum-1] = totalPred->GetBinContent(binNum);
82 
83 
84  // We need to include the statistical covariance matrix
85  // to do the chi^2 calculation or it'll be way off.
86  // we need to recalculate the sum for *every* new configuration of
87  // the IPrediction, unfortunately--
88  // and since we have to then invert it, this could be very costly.
89  // StackExchange to the rescue: somebody there wanted to solve a similar problem
90  // (where we already have the inverse of a matrix A
91  // and want to compute (A+B)^{-1] with some other invertible B):
92  // http://math.stackexchange.com/a/17780
93  // method has apparently been discovered in various forms by multiple people,
94  // but this version is due to K. Miller, Math.Mag. 54, 67 (1981)
95  //
96 
97 
98  // Strictly speaking, it's not possible to divide the prediction up
99  // into pieces and then use the covariance matrices for those pieces
100  // to compute a chi^2 because it's not possible to know which was which
101  // in the data. However, we'll cheat a little and do the following:
102  // for each component, we'll subtract the rest of the prediction
103  // and compare background-subtracted data to that component.
104  // Because we have to include the statistical covariance matrix every time,
105  // this over-counts the statistical uncertainty by a factor of
106  // roughly (N-1), where N is the number of components.
107  // Probably this isn't exactly the right statistical thing to do,
108  // but for a first approximation, we'll take an average of the chi^2s
109  // weighted by how many events they contribute to the prediction.
110 
111  double chi2 = 0;
112  for (auto it_comp = this->fCovMxInv.begin(); it_comp != this->fCovMxInv.end(); ++it_comp)
113  {
114  const auto & comp = it_comp->first;
115  const auto & compCovMxInv = it_comp->second;
116  const auto & predComp = *(predictions.at(comp));
117  TH1D hdata_compSub = *hdata - ( *totalPred - predComp );
118 
119 
120  // now recursively update the matrix
121  // by calculating the inverse of (A+B_k),
122  // where B_k is one of the N rank-1 matrices
123  // (when N is the rank of B) that B can be decomposed into.
124  // our case is even easier because we know B is diagonal
125  // (it's the statistical uncertainty covariance matrix).
126  TMatrixD Bk(stat_mx.GetNrows(), stat_mx.GetNcols());
127  TMatrixD CkInv(*compCovMxInv); // here we're setting what Miller calls C_1 = A
128 
129  // the rank of the stat mx is just the number of rows:
130  // it's diagonal and the bins are orthogonal.
131  // note that we're counting from zero instead of 1
132  // for convenience in the TMatrixDs
133  for (int k = 0; k < stat_mx.GetNrows(); k++)
134  {
135  // reset the last iteration so we can re-use
136  if (k > 0)
137  Bk[k-1][k-1] = 0;
138 
139  // set current iteration B_k
140  Bk[k][k] = stat_mx[k][k];
141 
142  // can short-circuit the trace in Miller's algorithm
143  // since our Bk come from a diagonal matrix:
144  // they're single entries on the diagonal
145  double gk = 1. / (1 + CkInv[k][k]*Bk[k][k]);
146 
147  CkInv -= gk * CkInv * Bk * CkInv;
148  }
149 
150  // event-weighted chi^2
151  chi2 += predComp.Integral() * Chi2CovMx(predictions.at(comp), &hdata_compSub, &CkInv);
152 
153  }
154 
155  // divide out the total number of events to the get the weighted chi^2
156  chi2 /= totalPred->Integral();
157 
158 
159  // finally, clean up all the temp hists we got from Spectrum::ToTH1()
160  HistCache::Delete(hdata);
161  HistCache::Delete(totalPred);
162  for (auto & predPair : predictions)
163  HistCache::Delete(predPair.second);
164 
165 
166  return chi2;
167  }
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:564
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
double chi2()
Double_t scale
Definition: plot.C:25
Interactions of both types.
Definition: IPrediction.h:42
static void Delete(TH1D *&h)
Definition: HistCache.cxx:220
#define pot
osc::OscCalcDumb calc
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:263
std::map< Component, std::unique_ptr< TMatrixD > > fCovMxInv
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
double Chi2CovMx(const TVectorD *e, const TVectorD *o, const TMatrixD *covmxinv)
Definition: Utilities.cxx:176
void ana::SingleSampleExperiment::Derivative ( osc::IOscCalc calc,
const SystShifts shift,
std::unordered_map< const ISyst *, double > &  dch 
) const
overridevirtualinherited

Reimplemented from ana::IExperiment.

Definition at line 111 of file SingleSampleExperiment.cxx.

References ana::HistCache::Delete(), ana::IPrediction::Derivative(), exit(), ana::SingleSampleExperiment::fBinErrors, ana::SingleSampleExperiment::fCosmic, ana::SingleSampleExperiment::fCosmicScaleError, ana::SingleSampleExperiment::fData, ana::SingleSampleExperiment::fMC, ana::SingleSampleExperiment::fPoissonError, ana::SystShifts::GetShift(), chisquared::hpred, MECModelEnuComparisons::i, it, ana::LogLikelihoodDerivative(), pot, ana::Spectrum::POT(), ana::SingleSampleExperiment::PredHistIncCosmics(), sigma(), and ana::Spectrum::ToTH1().

Referenced by ana::SingleSampleExperiment::ChiSq().

114  {
115 
116  const double pot = fData.POT();
117 
118  std::unordered_map<const ISyst*, std::vector<double>> dp;
119  for(auto it: dchi) dp[it.first] = {};
120  fMC->Derivative(calc, shift, pot, dp);
121 
122  if(dp.empty()){ // prediction doesn't implement derivatives
123  dchi.clear(); // pass on that info to our caller
124  return;
125  }
126 
127  TH1D* hpred = PredHistIncCosmics(calc, shift);
128  TH1D* hdata = fData.ToTH1(pot);
129 
130  for(auto& it: dchi){
131  if(it.first != &kCosmicBkgScaleSyst){
132  it.second += LogLikelihoodDerivative(hpred, hdata, dp[it.first]);
133  }
134  else{
135  const unsigned int N = fCosmic->GetNbinsX()+2;
136  const double* ca = fCosmic->GetArray();
137  std::vector<double> cosErr(N);
138  if(fPoissonError){
139  double sigma = shift.GetShift(&kCosmicBkgScaleSyst);
140  for(unsigned int i = 1; i<= fBinErrors.size(); i++){
141  double err = 0;
142  if(sigma<0) err = fBinErrors[i-1].second;
143  if(sigma>0) err = fBinErrors[i-1].first;
144  cosErr[i] = ca[i]*err;
145  }
146  }
147  else {
148  for(unsigned int i = 1; i <= N; ++i) cosErr[i] = ca[i]*fCosmicScaleError;
149  }
150  it.second += LogLikelihoodDerivative(hpred, hdata, cosErr);
151  }
152  }
153  exit(0);
154  HistCache::Delete(hpred);
155  HistCache::Delete(hdata);
156  }
set< int >::iterator it
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:564
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
TH1D * PredHistIncCosmics(osc::IOscCalc *calc, const SystShifts &syst) const
virtual void Derivative(osc::IOscCalc *calc, const SystShifts &shift, double pot, std::unordered_map< const ISyst *, std::vector< double >> &dchi) const
Definition: IPrediction.h:92
double LogLikelihoodDerivative(double e, double o, double dedx)
Definition: Utilities.cxx:154
static void Delete(TH1D *&h)
Definition: HistCache.cxx:220
#define pot
std::vector< std::pair< double, double > > fBinErrors
double POT() const
Definition: Spectrum.h:263
double sigma(TH1F *hist, double percentile)
exit(0)
std::unique_ptr< SingleSampleExperiment > ana::SingleSampleExperiment::LoadFrom ( TDirectory *  dir,
const std::string &  name 
)
staticinherited

Definition at line 196 of file SingleSampleExperiment.cxx.

References ana::assert(), dir, ana::Spectrum::LoadFrom(), ana::LoadFrom< IPrediction >(), mc, runNovaSAM::release, and runNovaSAM::ret.

Referenced by ana::LoadFrom< IExperiment >().

197  {
198  dir = dir->GetDirectory(name.c_str()); // switch to subdir
199  assert(dir);
200 
201  TObjString* ptag = (TObjString*)dir->Get("type");
202  assert(ptag);
203  assert(ptag->GetString() == "SingleSampleExperiment");
204  delete ptag;
205 
206  const IPrediction* mc = ana::LoadFrom<IPrediction>(dir, "mc").release();
207  const std::unique_ptr<Spectrum> data = Spectrum::LoadFrom(dir, "data");
208 
209  TH1D* cosmic = 0;
210  if(dir->Get("cosmic")) cosmic = (TH1D*)dir->Get("cosmic");
211 
212  delete dir;
213 
214  auto ret = std::make_unique<SingleSampleExperiment>(mc, *data);
215  if(cosmic) ret->fCosmic = cosmic;
216  return ret;
217  }
const XML_Char * name
Definition: expat.h:151
const XML_Char const XML_Char * data
Definition: expat.h:268
Sum up livetimes from individual cosmic triggers.
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:1071
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:33
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
virtual stan::math::var ana::IExperiment::LogLikelihood ( osc::IOscCalcAdjustableStan osc,
const SystShifts syst = kNoShift 
) const
inlinevirtualinherited

Reimplemented in test::GaussQuadExperiment.

Definition at line 38 of file IExperiment.h.

References ana::assert(), dir, and ana::IExperiment::SaveTo().

Referenced by ana::StanFitter::log_prob().

40  {
41  assert(false && "unimplemented");
42  return 0;
43  };
assert(nhit_max >=nhit_nbins)
stan::math::var ana::SingleSampleExperiment::LogLikelihood ( osc::_IOscCalcAdjustable< stan::math::var > *  osc,
const SystShifts syst = SystShifts::Nominal() 
) const
overridevirtualinherited

Definition at line 159 of file SingleSampleExperiment.cxx.

References ana::HistCache::Delete(), ana::SingleSampleExperiment::fData, ana::SingleSampleExperiment::fMC, ana::SystShifts::IsNominal(), ana::LogLikelihood(), ana::Spectrum::POT(), plot_validation_datamc::pred, ana::IPrediction::Predict(), ana::IPrediction::PredictSyst(), and ana::Spectrum::ToTH1().

Referenced by test_stanfit_systpulls().

161  {
162  auto pred = syst.IsNominal()
163  ? fMC->Predict(osc).ToBins(fData.POT())
164  : fMC->PredictSyst(osc, syst).ToBins(fData.POT());
165  TH1D* hdata = fData.ToTH1(fData.POT());
166 
167  // fully-qualified so that we get the one in StanUtils.h
168  auto ll = ana::LogLikelihood(pred, hdata) / -2.; // LogLikelihood(), confusingly, returns chi2=-2*LL
169 
170  HistCache::Delete(hdata);
171 
172  return ll;
173  }
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:564
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
stan::math::var LogLikelihood(const std::vector< stan::math::var > &exp, const TH1 *obs)
Variant that handles the prediction in the form of Stan vars.
Definition: StanUtils.cxx:11
static void Delete(TH1D *&h)
Definition: HistCache.cxx:220
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:85
double POT() const
Definition: Spectrum.h:263
void ana::CovMxExperiment::MakeInvCovMxs ( const std::map< Component, TMatrixD * > &  covmx)
inlineprotected

Definition at line 96 of file CovMxExperiment.h.

References fCovMxInv, and ana::SymmMxInverse().

Referenced by CovMxExperiment().

97  {
98  for (const auto & mxPair : covmx)
99  this->fCovMxInv.emplace( mxPair.first, SymmMxInverse(*(mxPair.second)) );
100 
101  }
std::unique_ptr< TMatrixD > SymmMxInverse(const TMatrixD &mx)
Definition: Utilities.cxx:219
std::map< Component, std::unique_ptr< TMatrixD > > fCovMxInv
CovMxExperiment& ana::CovMxExperiment::operator= ( const SingleSampleExperiment )
delete

Referenced by CovMxExperiment().

TH1D * ana::SingleSampleExperiment::PredHistIncCosmics ( osc::IOscCalc calc,
const SystShifts syst 
) const
protectedinherited

Definition at line 63 of file SingleSampleExperiment.cxx.

References ana::HistCache::Copy(), ana::HistCache::Delete(), ana::SingleSampleExperiment::fBinErrors, ana::SingleSampleExperiment::fCosmic, ana::SingleSampleExperiment::fCosmicScaleError, ana::SingleSampleExperiment::fData, ana::SingleSampleExperiment::fMC, ana::SingleSampleExperiment::fPoissonError, ana::SystShifts::GetShift(), chisquared::hpred, MECModelEnuComparisons::i, ana::Spectrum::POT(), plot_validation_datamc::pred, ana::IPrediction::PredictSyst(), scale, ana::SystShifts::SetShift(), sigma(), and ana::Spectrum::ToTH1().

Referenced by ana::SingleSampleExperiment::ChiSq(), ana::SingleSampleExperiment::Derivative(), and ana::SingleSampleExperiment::~SingleSampleExperiment().

65  {
66  SystShifts systNoCosmic = syst;
67  systNoCosmic.SetShift(&kCosmicBkgScaleSyst, 0);
68 
69  const Spectrum pred = fMC->PredictSyst(calc, systNoCosmic);
70  TH1D* hpred = pred.ToTH1(fData.POT());
71  if(fCosmic){
72  if(fCosmicScaleError != 0 && !fPoissonError){
73  const double scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
74  hpred->Add(fCosmic, scale);
75  }
76  if(fPoissonError){
77  TH1D* hcosm = HistCache::Copy(fCosmic);
78  double sigma = syst.GetShift(&kCosmicBkgScaleSyst);
79  for(int i = 1; i<= hcosm->GetNbinsX(); i++){
80  if (sigma>0) hcosm->SetBinContent(i, hcosm->GetBinContent(i)*(1+sigma*fBinErrors[i-1].first));
81  if (sigma<0) hcosm->SetBinContent(i, hcosm->GetBinContent(i)*(1+sigma*fBinErrors[i-1].second));
82  }
83  hpred->Add(hcosm);
84  HistCache::Delete(hcosm);
85  }
86  else{
87  hpred->Add(fCosmic);
88  }
89  }
90  return hpred;
91  }
static TH1D * Copy(const TH1D *h)
Definition: HistCache.cxx:138
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
Double_t scale
Definition: plot.C:25
static void Delete(TH1D *&h)
Definition: HistCache.cxx:220
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:85
std::vector< std::pair< double, double > > fBinErrors
double POT() const
Definition: Spectrum.h:263
double sigma(TH1F *hist, double percentile)
void ana::SingleSampleExperiment::SaveTo ( TDirectory *  dir,
const std::string &  name 
) const
overridevirtualinherited

Reimplemented from ana::IExperiment.

Definition at line 175 of file SingleSampleExperiment.cxx.

References dir, ana::SingleSampleExperiment::fCosmic, ana::SingleSampleExperiment::fData, ana::SingleSampleExperiment::fMC, ana::IPrediction::SaveTo(), ana::Spectrum::SaveTo(), and tmp.

176  {
177  TDirectory* tmp = dir;
178 
179  dir = dir->mkdir(name.c_str()); // switch to subdir
180  dir->cd();
181 
182  TObjString("SingleSampleExperiment").Write("type");
183 
184  fMC->SaveTo(dir, "mc");
185  fData.SaveTo(dir, "data");
186 
187  if(fCosmic) fCosmic->Write("cosmic");
188 
189  dir->Write();
190  delete dir;
191 
192  tmp->cd();
193  }
const XML_Char * name
Definition: expat.h:151
Float_t tmp
Definition: plot.C:36
virtual void SaveTo(TDirectory *dir, const std::string &name) const
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:1040
TDirectory * dir
Definition: macro.C:5

Member Data Documentation

std::vector< std::pair <double, double> > ana::SingleSampleExperiment::fBinErrors
protectedinherited
TH1D* ana::SingleSampleExperiment::fCosmic
protectedinherited
double ana::SingleSampleExperiment::fCosmicScaleError
protectedinherited
std::map<Component, std::unique_ptr<TMatrixD> > ana::CovMxExperiment::fCovMxInv
protected

Definition at line 103 of file CovMxExperiment.h.

Referenced by ChiSq(), and MakeInvCovMxs().

Spectrum ana::SingleSampleExperiment::fData
protectedinherited
const IPrediction* ana::SingleSampleExperiment::fMC
protectedinherited
bool ana::SingleSampleExperiment::fPoissonError
protectedinherited

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