SingleSampleExperiment.cxx
Go to the documentation of this file.
2 
7 
9 
10 #include "TDirectory.h"
11 #include "TObjString.h"
12 #include "TH1.h"
13 #include "TFeldmanCousins.h"
14 
15 namespace ana
16 {
18 
19  //----------------------------------------------------------------------
21  const Spectrum& data,
22  const Spectrum& cosmic,
23  double cosmicScaleError,
24  bool PoissonError)
25  : fMC(pred), fData(data),
26  fCosmic(cosmic.ToTH1(data.Livetime(), kLivetime)),
27  fCosmicScaleError(cosmicScaleError),
28  fPoissonError(PoissonError)
29  {
30  // Fill the vector of errors in case of PoissonError = true
31  if(fPoissonError){
32  TH1D* cosmicErrHist = cosmic.ToTH1(cosmic.Livetime(), kLivetime);
33 
34  TFeldmanCousins fc(0.6827);
35  for(int i = 1; i<= cosmicErrHist->GetNbinsX(); i++){
36  double y = cosmicErrHist->GetBinContent(i);
37  double errup = 0;
38  double errlow = 0;
39  if(y!=0){
40  if ( y < 50 ) errlow = (y - fc.CalculateLowerLimit(y,0))/y;
41  else errlow = (sqrt(y))/y;
42  if ( y < 30 ) errup = (fc.CalculateUpperLimit(y,0)-y)/y;
43  else errup = (sqrt(y))/y;
44  }
45  fBinErrors.push_back({errup, errlow});
46  }
47 
48  HistCache::Delete(cosmicErrHist);
49  }
50 
51  }
52 
53  //----------------------------------------------------------------------
55  {
57  }
58 
59  //----------------------------------------------------------------------
62  const SystShifts& syst) const
63  {
64  SystShifts systNoCosmic = syst;
65  systNoCosmic.SetShift(&kCosmicBkgScaleSyst, 0);
66 
67  const Spectrum pred = fMC->PredictSyst(calc, systNoCosmic);
68  TH1D* hpred = pred.ToTH1(fData.POT());
69  if(fCosmic){
70  if(fCosmicScaleError != 0 && !fPoissonError){
71  const double scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
72  hpred->Add(fCosmic, scale);
73  }
74  if(fPoissonError){
75  TH1D* hcosm = HistCache::Copy(fCosmic);
76  double sigma = syst.GetShift(&kCosmicBkgScaleSyst);
77  for(int i = 1; i<= hcosm->GetNbinsX(); i++){
78  if (sigma>0) hcosm->SetBinContent(i, hcosm->GetBinContent(i)*(1+sigma*fBinErrors[i-1].first));
79  if (sigma<0) hcosm->SetBinContent(i, hcosm->GetBinContent(i)*(1+sigma*fBinErrors[i-1].second));
80  }
81  hpred->Add(hcosm);
82  HistCache::Delete(hcosm);
83  }
84  else{
85  hpred->Add(fCosmic);
86  }
87  }
88  return hpred;
89  }
90 
91  //----------------------------------------------------------------------
93  const SystShifts& syst) const
94  {
95  TH1D* hpred = PredHistIncCosmics(calc, syst);
96  TH1D* hdata = fData.ToTH1(fData.POT());
97 
98  // full namespace qualification to avoid degeneracy with method inherited from IExperiment
99  const double ll = ana::LogLikelihood(hpred, hdata);
100 
101  HistCache::Delete(hpred);
102  HistCache::Delete(hdata);
103 
104  return ll;
105  }
106 
107  //----------------------------------------------------------------------
110  const SystShifts& shift,
111  std::unordered_map<const ISyst*, double>& dchi) const
112  {
113 
114  const double pot = fData.POT();
115 
116  std::unordered_map<const ISyst*, std::vector<double>> dp;
117  for(auto it: dchi) dp[it.first] = {};
118  fMC->Derivative(calc, shift, pot, dp);
119 
120  if(dp.empty()){ // prediction doesn't implement derivatives
121  dchi.clear(); // pass on that info to our caller
122  return;
123  }
124 
125  TH1D* hpred = PredHistIncCosmics(calc, shift);
126  TH1D* hdata = fData.ToTH1(pot);
127 
128  for(auto& it: dchi){
129  if(it.first != &kCosmicBkgScaleSyst){
130  it.second += LogLikelihoodDerivative(hpred, hdata, dp[it.first]);
131  }
132  else{
133  const unsigned int N = fCosmic->GetNbinsX()+2;
134  const double* ca = fCosmic->GetArray();
135  std::vector<double> cosErr(N);
136  if(fPoissonError){
137  double sigma = shift.GetShift(&kCosmicBkgScaleSyst);
138  for(unsigned int i = 1; i<= fBinErrors.size(); i++){
139  double err = 0;
140  if(sigma<0) err = fBinErrors[i-1].second;
141  if(sigma>0) err = fBinErrors[i-1].first;
142  cosErr[i] = ca[i]*err;
143  }
144  }
145  else {
146  for(unsigned int i = 1; i <= N; ++i) cosErr[i] = ca[i]*fCosmicScaleError;
147  }
148  it.second += LogLikelihoodDerivative(hpred, hdata, cosErr);
149  }
150  }
151  exit(0);
152  HistCache::Delete(hpred);
153  HistCache::Delete(hdata);
154  }
155 
156  //----------------------------------------------------------------------
158  const SystShifts &syst) const
159  {
160  auto pred = syst.IsNominal()
161  ? fMC->Predict(osc).ToBins(fData.POT())
162  : fMC->PredictSyst(osc, syst).ToBins(fData.POT());
163  TH1D* hdata = fData.ToTH1(fData.POT());
164 
165  // fully-qualified so that we get the one in StanUtils.h
166  auto ll = ana::LogLikelihood(pred, hdata) / -2.; // LogLikelihood(), confusingly, returns chi2=-2*LL
167 
168  HistCache::Delete(hdata);
169 
170  return ll;
171  }
172  //----------------------------------------------------------------------
173  void SingleSampleExperiment::SaveTo(TDirectory* dir) const
174  {
175  TDirectory* tmp = dir;
176 
177  dir->cd();
178  TObjString("SingleSampleExperiment").Write("type");
179 
180  fMC->SaveTo(dir->mkdir("mc"));
181  fData.SaveTo(dir->mkdir("data"));
182 
183  if(fCosmic) fCosmic->Write("cosmic");
184 
185  tmp->cd();
186  }
187 
188  //----------------------------------------------------------------------
189  std::unique_ptr<SingleSampleExperiment> SingleSampleExperiment::LoadFrom(TDirectory* dir)
190  {
191  TObjString* ptag = (TObjString*)dir->Get("type");
192  assert(ptag);
193  assert(ptag->GetString() == "SingleSampleExperiment");
194 
195  assert(dir->GetDirectory("mc"));
196  assert(dir->GetDirectory("data"));
197 
198 
199  const IPrediction* mc = ana::LoadFrom<IPrediction>(dir->GetDirectory("mc")).release();
200  const std::unique_ptr<Spectrum> data = Spectrum::LoadFrom(dir->GetDirectory("data"));
201 
202  TH1D* cosmic = 0;
203  if(dir->Get("cosmic")) cosmic = (TH1D*)dir->Get("cosmic");
204 
205  auto ret = std::make_unique<SingleSampleExperiment>(mc, *data);
206  if(cosmic) ret->fCosmic = cosmic;
207  return ret;
208  }
209 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: CutFlow_header.h:5
set< int >::iterator it
static TH1D * Copy(const TH1D *h)
Definition: HistCache.cxx:138
bool IsNominal() const
Definition: SystShifts.h:43
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 void SaveTo(TDirectory *dir) const
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
virtual stan::math::var LogLikelihood(osc::_IOscCalculatorAdjustable< stan::math::var > *osc, const SystShifts &syst=SystShifts::Nominal()) const override
T sqrt(T number)
Definition: d0nt_math.hpp:156
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
virtual double ChiSq(osc::IOscCalculatorAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
virtual void Derivative(osc::IOscCalculator *calc, const SystShifts &shift, double pot, std::unordered_map< const ISyst *, std::vector< double >> &dchi) const
Definition: IPrediction.h:92
Float_t tmp
Definition: plot.C:36
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
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
virtual Spectrum PredictSyst(osc::IOscCalculator *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:98
const XML_Char const XML_Char * data
Definition: expat.h:268
double LogLikelihoodDerivative(double e, double o, double dedx)
Definition: Utilities.cxx:154
Double_t scale
Definition: plot.C:25
Sum up livetimes from individual cosmic triggers.
static void Delete(TH1D *&h)
Definition: HistCache.cxx:220
#define pot
T GetShift(const ISyst *syst) const
std::vector< std::pair< double, double > > fBinErrors
double POT() const
Definition: Spectrum.h:263
osc::OscCalculatorDumb calc
Oscillation probability calculators.
Definition: Calcs.h:5
virtual void SaveTo(TDirectory *dir) const override
double sigma(TH1F *hist, double percentile)
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
virtual void Derivative(osc::IOscCalculator *calc, const SystShifts &shift, std::unordered_map< const ISyst *, double > &dch) const override
TH1D * PredHistIncCosmics(osc::IOscCalculator *calc, const SystShifts &syst) const
TDirectory * dir
Definition: macro.C:5
string syst
Definition: plotSysts.py:176
exit(0)
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< Spectrum > LoadFrom(TDirectory *dir)
Definition: Spectrum.cxx:1066
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:70
static std::unique_ptr< SingleSampleExperiment > LoadFrom(TDirectory *dir)
void SaveTo(TDirectory *dir) const
Definition: Spectrum.cxx:1040