SingleSampleExperiment.cxx
Go to the documentation of this file.
2 
7 
8 #include "OscLib/IOscCalc.h"
9 
10 #include "Utilities/func/Stan.h"
11 
12 #include "TDirectory.h"
13 #include "TObjString.h"
14 #include "TH1.h"
15 #include "TFeldmanCousins.h"
16 
17 namespace ana
18 {
19  REGISTER_LOADFROM("SingleSampleExperiment", IExperiment, SingleSampleExperiment);
20 
22 
23  //----------------------------------------------------------------------
25  const Spectrum& data,
26  const Spectrum& cosmic,
27  double cosmicScaleError,
28  bool PoissonError)
29  : fMC(pred), fData(data),
30  fCosmic(cosmic),
31  fCosmicScaleError(cosmicScaleError),
32  fPoissonError(PoissonError)
33  {
34  // Fill the vector of errors in case of PoissonError = true
35  if(fPoissonError){
36  // Note this is scaled to the cosmic's intrinsic stats level
37  Eigen::ArrayXd cosmicErrHist = cosmic.GetEigen(cosmic.Livetime(), kLivetime);
38 
39  TFeldmanCousins fc(0.6827);
40  for(int i = 1; i<= cosmicErrHist.size()-2; i++){
41  double y = cosmicErrHist[i];
42  double errup = 0;
43  double errlow = 0;
44  if(y!=0){
45  if ( y < 50 ) errlow = (y - fc.CalculateLowerLimit(y,0))/y;
46  else errlow = (sqrt(y))/y;
47  if ( y < 30 ) errup = (fc.CalculateUpperLimit(y,0)-y)/y;
48  else errup = (sqrt(y))/y;
49  }
50  fBinErrors.push_back({errup, errlow});
51  }
52  }
53 
54  }
55 
56  //----------------------------------------------------------------------
58  {
59  }
60 
61  //----------------------------------------------------------------------
62  template <typename T>
63  Eigen::Array<T, Eigen::Dynamic, 1>
65  const SystShifts& syst) const
66  {
67  using EigenArray = Eigen::Array<T, Eigen::Dynamic, 1>;
68 
69  SystShifts systNoCosmic = syst;
70  systNoCosmic.RemoveShift(&kCosmicBkgScaleSyst);
71 
72  Spectrum specPred = fMC->PredictSyst(calc, systNoCosmic);
73  EigenArray apred;
74  if (specPred.HasStan())
75  {
76  if constexpr (std::is_same_v<T, stan::math::var>)
77  apred = specPred.GetEigenStan(fData.POT());
78  }
79  else
80  apred = specPred.GetEigen(fData.POT());
81 
82  if(fCosmic.has_value()){
83  // note the conversion here from double to stan::math::var if T is stan::math::var
84  EigenArray acosm = fCosmic->GetEigen(fData.Livetime(), kLivetime);
85  if(fCosmicScaleError != 0 && !fPoissonError){
86  const T scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
87  acosm *= scale;
88  }
89  if(fPoissonError){
90  T sigma = syst.GetShift(&kCosmicBkgScaleSyst);
91  for(int i = 1; i<= acosm.size()-2; i++){
92  if (sigma>0) acosm[i] *= (1+sigma*fBinErrors[i-1].first);
93  if (sigma<0) acosm[i] *= (1+sigma*fBinErrors[i-1].second);
94  }
95  }
96 
97  apred += acosm;
98  }
99 
100  return apred;
101  }
102 
103  template Eigen::ArrayXd
105  const SystShifts& syst) const;
106  template Eigen::ArrayXstan
108  const SystShifts& syst) const;
109 
110  //----------------------------------------------------------------------
112  const SystShifts& syst) const
113  {
114  Eigen::ArrayXd pred = PredHistIncCosmics(calc, syst);
115  Eigen::ArrayXd data = fData.GetEigen(fData.POT());
116 
117  // full namespace qualification to avoid degeneracy with method inherited from IExperiment
118  return ana::LogLikelihood(pred, data);
119  }
120 
121  //----------------------------------------------------------------------
123  const SystShifts &syst) const
124  {
126  Eigen::ArrayXd data = fData.GetEigen(fData.POT());
127 
128  // fully-qualified so that we get the one in StanUtils.h
129  return ana::LogLikelihood(pred, data) / -2.; // LogLikelihood(), confusingly, returns chi2=-2*LL
130  }
131 
132  //----------------------------------------------------------------------
133  void SingleSampleExperiment::SaveTo(TDirectory* dir, const std::string& name) const
134  {
135  TDirectory* tmp = dir;
136 
137  dir = dir->mkdir(name.c_str()); // switch to subdir
138  dir->cd();
139 
140  TObjString("SingleSampleExperiment").Write("type");
141 
142  fMC->SaveTo(dir, "mc");
143  fData.SaveTo(dir, "data");
144 
145  if(fCosmic.has_value()) fCosmic->SaveTo(dir, "cosmic");
146 
147  dir->Write();
148  delete dir;
149 
150  tmp->cd();
151  }
152 
153  //----------------------------------------------------------------------
154  std::unique_ptr<SingleSampleExperiment> SingleSampleExperiment::LoadFrom(TDirectory* dir, const std::string& name)
155  {
156  dir = dir->GetDirectory(name.c_str()); // switch to subdir
157  assert(dir);
158 
159  TObjString* ptag = (TObjString*)dir->Get("type");
160  assert(ptag);
161  assert(ptag->GetString() == "SingleSampleExperiment");
162  delete ptag;
163 
165  const std::unique_ptr<Spectrum> data = Spectrum::LoadFrom(dir, "data");
166 
167  auto ret = std::make_unique<SingleSampleExperiment>(mc, *data);
168 
169  if(dir->Get("cosmic")){
170  ret->fCosmic = *ana::LoadFrom<Spectrum>(dir, "cosmic");
171  }
172 
173  delete dir;
174 
175  return ret;
176  }
177 }
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::optional< Spectrum > fCosmic
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:188
Float_t tmp
Definition: plot.C:36
virtual void SaveTo(TDirectory *dir, const std::string &name) const
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
osc::OscCalcDumb calc
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const XML_Char const XML_Char * data
Definition: expat.h:268
Double_t scale
Definition: plot.C:25
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Sum up livetimes from individual cosmic triggers.
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
T GetShift(const ISyst *syst) const
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:49
std::vector< std::pair< double, double > > fBinErrors
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:18
double POT() const
Definition: Spectrum.h:219
Oscillation probability calculators.
Definition: Calcs.h:5
double sigma(TH1F *hist, double percentile)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Eigen::Array< stan::math::var, Eigen::Dynamic, 1 > ArrayXstan
Definition: StanUtils.h:7
TDirectory * dir
Definition: macro.C:5
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
assert(nhit_max >=nhit_nbins)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Dummy syst to communicate with SingleSampleExperiment.
static std::unique_ptr< SingleSampleExperiment > LoadFrom(TDirectory *dir, const std::string &name)
double T
Definition: Xdiff_gwt.C:5
Eigen::Array< T, Eigen::Dynamic, 1 > PredHistIncCosmics(osc::_IOscCalc< T > *calc, const SystShifts &syst) const
virtual stan::math::var LogLikelihood(osc::IOscCalcAdjustableStan *osc, const SystShifts &syst=kNoShift) const override
void RemoveShift(const ISyst *syst)
Definition: SystShifts.cxx:73
bool HasStan() const
Definition: Spectrum.h:186
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:222
const Eigen::ArrayXstan & GetEigenStan() const
Definition: Spectrum.h:189