PredictionScaleComp.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Cut.h"
4 
5 #include "TDirectory.h"
6 #include "TObjString.h"
7 #include "TH1.h"
8 #include "TVectorD.h"
9 #include "TParameter.h"
10 
11 #include <cassert>
12 
13 namespace ana
14 {
15  namespace
16  {
17  // We need some way to identify what components the user wants to vary in
18  // the fit. We just create N copies of this class, each one identifying one
19  // component. They're distinguished by their pointer values.
20  class DummyScaleCompSyst : public ana::ISyst
21  {
22  public:
23  DummyScaleCompSyst(const std::string& uniqifier, int i)
24  : ISyst("CompNormShift_" + std::to_string(i)+'_'+uniqifier,
25  "Component Normalization Shift " + std::to_string(i))
26  //, fN(i)
27  {}
28 
29  void Shift(double sigma, caf::SRProxy* sr,
30  double& weight) const override
31  {
32  }
33  private:
34  //int fN;
35  };
36  }
37 
38  //----------------------------------------------------------------------
40  const HistAxis& axis,
41  Cut cut,
42  std::vector< std::pair<Cut, bool> > truthcuts,
43  const SystShifts& shift,
44  const Var& wei)
45  : fComplementCut(kNoCut)
46  {
47  const std::string uniq = TString::Format("%p", (void*)this).Data();
48 
49  assert(truthcuts.size()>0 && "Please give at least one truth selection.");
50  for (int i = 0; i < int(truthcuts.size()); i++){
51  fSpectra.push_back(new Spectrum(loader, axis,
52  cut && truthcuts[i].first, shift, wei));
53  fSysts.push_back(new DummyScaleCompSyst(uniq, i));
54  fComplementCut = fComplementCut && !truthcuts[i].first;
55  fIsComponentFloated.push_back(truthcuts[i].second);
56  }
57 
58  // The idea is that if truthcuts are exhaustive, this Spectrum should wind
59  // up empty
60  fComplement = new Spectrum(loader, axis,
61  cut && fComplementCut, shift, wei);
62 
63  }
64 
65 
66  //----------------------------------------------------------------------
68  const HistAxis& axis,
69  Cut cut,
70  std::vector<Cut> truthcuts,
71  const SystShifts& shift,
72  const Var& wei)
74  {
75  const std::string uniq = TString::Format("%p", (void*)this).Data();
76 
77  assert(truthcuts.size()>0 && "Please give at least one truth selection.");
78  for (int i = 0; i < int(truthcuts.size()); i++){
79  fSpectra.push_back(new Spectrum(loader, axis,
80  cut && truthcuts[i], shift, wei));
81  fSysts.push_back(new DummyScaleCompSyst(uniq, i));
82  fComplementCut = fComplementCut && !truthcuts[i];
83  fIsComponentFloated.push_back(true);
84  }
85 
86  // The idea is that if truthcuts are exhaustive, this Spectrum should wind
87  // up empty
88  fComplement = new Spectrum(loader, axis,
89  cut && fComplementCut, shift, wei);
90  }
91 
92 
93  //----------------------------------------------------------------------
95  const HistAxis& axis1,
96  const HistAxis& axis2,
97  Cut cut,
98  std::vector<Cut> truthcuts,
99  const SystShifts& shift,
100  const Var& wei)
102  {
103  const std::string uniq = TString::Format("%p", (void*)this).Data();
104 
105  assert(truthcuts.size()>0 && "Please give at least one truth selection.");
106  for (int i = 0; i < int(truthcuts.size()); i++){
107  fSpectra.push_back(new Spectrum(loader, axis1, axis2,
108  cut && truthcuts[i], shift, wei));
109  fSysts.push_back(new DummyScaleCompSyst(uniq, i));
110  fComplementCut = fComplementCut && !truthcuts[i];
111  fIsComponentFloated.push_back(true);
112  }
113 
114  // The idea is that if truthcuts are exhaustive, this Spectrum should wind
115  // up empty
116  fComplement = new Spectrum(loader, axis1, axis2,
117  cut && fComplementCut, shift, wei);
118  }
119 
120  //----------------------------------------------------------------------
122  std::vector<Spectrum*> spectra,
123  std::vector<bool> isfloated)
124 
125  : fSpectra(spectra),
127  fComplement(complement)
128  {
129  const std::string uniq = TString::Format("%p", (void*)this).Data();
130  fIsComponentFloated=isfloated;
131  for (int i = 0; i < int(spectra.size()); i++)
132  {
133  fSysts.push_back(new DummyScaleCompSyst(uniq, i));
134  }
135  }
136 
137  //----------------------------------------------------------------------
139  {
140  Spectrum ret = *fSpectra[0];
141  for (int i = 1; i < int(fSpectra.size()); i++) ret += *fSpectra[i];
142  return ret;
143  }
144 
145  //----------------------------------------------------------------------
147  const SystShifts& shift) const
148  {
149  double pot = fSpectra[0]->POT();
150  Eigen::ArrayXd hist = fSpectra[0]->GetEigen(pot);
151  double scale = fIsComponentFloated[0] ? 1 + shift.GetShift(fSysts[0]) : 1;
152  hist *= scale;
153  for (int i = 1; i < int(fSpectra.size()); i++){
154  Eigen::ArrayXd curhist = fSpectra[i]->GetEigen(pot);
155  scale = fIsComponentFloated[i] ? 1 + shift.GetShift(fSysts[i]) : 1;
156  hist += curhist * scale;
157  }
158  return Spectrum(std::move(hist),
159  HistAxis(fSpectra[0]->GetLabels(),
160  fSpectra[0]->GetBinnings()),
161  pot, fSpectra[0]->Livetime());
162  }
163 
164  //----------------------------------------------------------------------
165  void PredictionScaleComp::SaveTo(TDirectory* dir, const std::string& name) const
166  {
167  TDirectory* tmp = gDirectory;
168 
169  dir = dir->mkdir(name.c_str()); // switch to subdir
170  dir->cd();
171 
172  TObjString("PredictionScaleComp").Write("type");
173 
174  // Save the number of MC spectra as an entry in a dummy histogram;
175  TVectorD nspectra(1);
176  nspectra[0] = fSpectra.size();
177  nspectra.Write("NSpectra");
178  // Save Vector of Spectrum*'s
179  for (int i = 0; i < int(fSpectra.size()); i++){
180  fSpectra[i]->SaveTo(dir, "Spectrum_"+std::to_string(i));
181  TParameter<bool> isfloated;
182  isfloated.SetVal(fIsComponentFloated[i]);
183  isfloated.Write(("IsComponentFloated_"+std::to_string(i)).c_str());
184  }
185  fComplement->SaveTo(dir, "Complement");
186 
187  dir->Write();
188  delete dir;
189 
190  tmp->cd();
191  }
192 
193  //----------------------------------------------------------------------
194  std::unique_ptr<PredictionScaleComp> PredictionScaleComp::LoadFrom(TDirectory* dir, const std::string& name)
195  {
196  dir = dir->GetDirectory(name.c_str()); // switch to subdir
197  assert(dir);
198 
199  Spectrum* complement = ana::LoadFrom<Spectrum>
200  (dir, "Complement").release();
201 
202  // Get number of spectra we expect
203  TVectorD* nspectra = (TVectorD*)dir->Get("NSpectra");
204  assert(nspectra->GetNoElements()==1 && "Vector of N MC Bkgds is not of size 1.");
205  // Load spectra
206  std::vector<Spectrum*> spectra;
207  std::vector<bool> isfloated;
208  for (int i = 0; i < (*nspectra)[0]; i++){
209  spectra.push_back(ana::LoadFrom<Spectrum>(
210  dir, "Spectrum_"+std::to_string(i)).release());
211  isfloated.push_back((bool)dir->Get(("IsComponentFloated_"+std::to_string(i)).c_str()));
212  }
213  delete nspectra;
214 
215  delete dir;
216 
217  return std::unique_ptr<PredictionScaleComp>(
218  new PredictionScaleComp(complement, spectra, isfloated));
219  }
220 
221 }
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2038
std::vector< Spectrum * > fSpectra
Float_t tmp
Definition: plot.C:36
std::vector< const ISyst * > fSysts
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
Spectrum Predict(osc::IOscCalc *osc) const override
Double_t scale
Definition: plot.C:25
Definition: Shift.h:6
#define pot
T GetShift(const ISyst *syst) const
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:578
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:527
double sigma(TH1F *hist, double percentile)
Base class for the various types of spectrum loader.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
TDirectory * dir
Definition: macro.C:5
static std::unique_ptr< PredictionScaleComp > LoadFrom(TDirectory *dir, const std::string &name)
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::vector< bool > fIsComponentFloated
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:117
Spectrum PredictSyst(osc::IOscCalc *osc, const SystShifts &syst) const override
PredictionScaleComp(SpectrumLoaderBase &loader, const HistAxis &axis, Cut cut, std::vector< Cut > truthcuts, const SystShifts &shift=kNoShift, const Var &wei=kUnweighted)
static constexpr Double_t sr
Definition: Munits.h:164