SingleSampleExperiment.cxx
Go to the documentation of this file.
2 
6 
8 
9 #include "TDirectory.h"
10 #include "TObjString.h"
11 #include "TH1.h"
12 
13 namespace ana
14 {
16 
17  //----------------------------------------------------------------------
19  const Spectrum& data,
20  const Spectrum& cosmic,
21  double cosmicScaleError)
22  : fMC(pred), fData(data),
23  fCosmic(cosmic.ToTH1(data.Livetime(), kLivetime)),
24  fCosmicScaleError(cosmicScaleError)
25  {
26  }
27 
28  //----------------------------------------------------------------------
30  const Spectrum& data,
31  const TH1D* cosmic,
32  double cosmicScaleError,
33  double fromlivetime,
34  double tolivetime)
35  : fMC(pred), fData(data), fCosmic(HistCache::Copy(cosmic)),
36  fCosmicScaleError(cosmicScaleError)
37  {
38  //default for backward capacity
39  if(fromlivetime ==1 && tolivetime ==1){
40 // std::cout<< "Single Experiment Warning: using TH1D for cosmic, please double check POT/Livetime... And switch to Spectrum if possible."<<std::endl;
41  }
42  //other cases, always positive livetimes
43  else if(fromlivetime > 0 && tolivetime > 0){
44 
45  std::cout<<"Single Experiment rescaling cosmic Livetime from: "<<
46  fromlivetime<<" to: "<<
47  tolivetime<<
48  "Please use Spectrum for cosmic next time."<<std::endl;
49 
50  fCosmic->Scale(tolivetime/fromlivetime);
51  }
52  // others not allowed
53  else{
54  std::cout<<"Single Experiment Warning: please double check cosmic POT/Livetime setting."<<std::endl;
55  abort();
56  }
57 
58  }
59 
60  //----------------------------------------------------------------------
62  {
64  }
65 
66  //----------------------------------------------------------------------
69  const SystShifts& syst) const
70  {
71  SystShifts systNoCosmic = syst;
72  systNoCosmic.SetShift(&kCosmicBkgScaleSyst, 0);
73 
74  const Spectrum pred = fMC->PredictSyst(calc, systNoCosmic);
75 
76  TH1D* hpred = pred.ToTH1(fData.POT());
77 
78  if(fCosmic){
79  if(fCosmicScaleError != 0){
80  const double scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
81  hpred->Add(fCosmic, scale);
82  }
83  else{
84  hpred->Add(fCosmic);
85  }
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  const double ll = LogLikelihood(hpred, hdata);
99 
100  HistCache::Delete(hpred);
101  HistCache::Delete(hdata);
102 
103  return ll;
104  }
105 
106  //----------------------------------------------------------------------
109  const SystShifts& shift,
110  std::unordered_map<const ISyst*, double>& dchi) const
111  {
112  const double pot = fData.POT();
113 
114  std::unordered_map<const ISyst*, std::vector<double>> dp;
115  for(auto it: dchi) dp[it.first] = {};
116  fMC->Derivative(calc, shift, pot, dp);
117 
118  if(dp.empty()){ // prediction doesn't implement derivatives
119  dchi.clear(); // pass on that info to our caller
120  return;
121  }
122 
123  TH1D* hpred = PredHistIncCosmics(calc, shift);
124  TH1D* hdata = fData.ToTH1(pot);
125 
126  for(auto& it: dchi){
127  if(it.first != &kCosmicBkgScaleSyst){
128  it.second += LogLikelihoodDerivative(hpred, hdata, dp[it.first]);
129  }
130  else{
131  const unsigned int N = fCosmic->GetNbinsX()+2;
132  const double* ca = fCosmic->GetArray();
133  std::vector<double> cosErr(N);
134  for(unsigned int i = 0; i < N; ++i) cosErr[i] = ca[i]*fCosmicScaleError;
135  it.second += LogLikelihoodDerivative(hpred, hdata, cosErr);
136  }
137  }
138 
139  HistCache::Delete(hpred);
140  HistCache::Delete(hdata);
141  }
142 
143  //----------------------------------------------------------------------
144  void SingleSampleExperiment::SaveTo(TDirectory* dir) const
145  {
146  TDirectory* tmp = dir;
147 
148  dir->cd();
149  TObjString("SingleSampleExperiment").Write("type");
150 
151  fMC->SaveTo(dir->mkdir("mc"));
152  fData.SaveTo(dir->mkdir("data"));
153 
154  if(fCosmic) fCosmic->Write("cosmic");
155 
156  tmp->cd();
157  }
158 
159  //----------------------------------------------------------------------
160  std::unique_ptr<SingleSampleExperiment> SingleSampleExperiment::LoadFrom(TDirectory* dir)
161  {
162  TObjString* ptag = (TObjString*)dir->Get("type");
163  assert(ptag);
164  assert(ptag->GetString() == "SingleSampleExperiment");
165 
166  assert(dir->GetDirectory("mc"));
167  assert(dir->GetDirectory("data"));
168 
169 
170  const IPrediction* mc = ana::LoadFrom<IPrediction>(dir->GetDirectory("mc")).release();
171  const std::unique_ptr<Spectrum> data = Spectrum::LoadFrom(dir->GetDirectory("data"));
172 
173  TH1D* cosmic = 0;
174  if(dir->Get("cosmic")) cosmic = (TH1D*)dir->Get("cosmic");
175 
176  auto ret = std::make_unique<SingleSampleExperiment>(mc, *data);
177  if(cosmic) ret->fCosmic = cosmic;
178  return ret;
179  }
180 }
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
set< int >::iterator it
General interface to oscillation calculators.
Definition: FwdDeclare.h: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:553
virtual void SaveTo(TDirectory *dir) const
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
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
Helper for Spectrum.
Definition: HistCache.h:27
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
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
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0)
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
General interface to any calculator that lets you set the parameters.
Sum up livetimes from individual cosmic triggers.
static void Delete(TH1D *&h)
Definition: HistCache.cxx:215
#define pot
T GetShift(const ISyst *syst) const
double POT() const
Definition: Spectrum.h:263
virtual void SaveTo(TDirectory *dir) const override
OStream cout
Definition: OStream.cxx:6
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
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:1055
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:1029