OscillatableSpectrum.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Core/OscCurve.h"
5 #include "CAFAna/Core/Ratio.h"
6 #include "Utilities/func/Stan.h"
8 
10 
11 #include "OscLib/IOscCalc.h"
12 
13 #include "TDirectory.h"
14 #include "TH2.h"
15 #include "TMD5.h"
16 #include "TObjString.h"
17 
18 #include <cassert>
19 #include <memory>
20 
21 namespace ana
22 {
23  // Duplicate here because we can't include Vars.h
24  const Var kTrueE([](const caf::SRProxy* sr)
25  {
26  assert(sr->mc.nnu == 1);
27  return sr->mc.nu[0].E;
28  });
29 
30  const Cut kHasNu([](const caf::SRProxy* sr)
31  {
32  if(sr->mc.nnu == 0) return false;
33  assert(sr->mc.nnu == 1);
34  return true;
35  });
36 
37  //----------------------------------------------------------------------
41  const Var& var,
42  const Cut& cut,
43  const SystShifts& shift,
44  const Var& wei)
45  : ReweightableSpectrum(loader,
46  HistAxis(label, bins, var),
47  HistAxis("True Energy (GeV)", kTrueEnergyBins, kTrueE),
48  kHasNu && cut, shift, wei)
49  {
50  }
51 
52  //----------------------------------------------------------------------
54  const HistAxis& axis,
55  const Cut& cut,
56  const SystShifts& shift,
57  const Var& wei)
58  : ReweightableSpectrum(loader, axis,
59  HistAxis("True Energy (GeV)", kTrueEnergyBins, kTrueE),
60  kHasNu && cut, shift, wei)
61  {
62  }
63 
64  //----------------------------------------------------------------------
66  const HistAxis& recoAxis,
67  double pot, double livetime)
68  : ReweightableSpectrum(std::move(mat), recoAxis,
69  HistAxis("True Energy (GeV)", kTrueEnergyBins, kTrueE),
70  pot, livetime)
71  {
72  }
73 
74  //----------------------------------------------------------------------
76  {
77  }
78 
79  //----------------------------------------------------------------------
82  {
83  if(rhs.fCache->hash){
84  fCache->spect = rhs.fCache->spect;
85  fCache->hash = std::make_unique<TMD5>(*rhs.fCache->hash);
86  }
87 
88  assert( rhs.fReferences.empty() ); // Copying with pending loads is unexpected
89  }
90 
91  //----------------------------------------------------------------------
94  {
95  if(rhs.fCache->hash){
96  fCache->spect = std::move(rhs.fCache->spect);
97  fCache->hash = std::move(rhs.fCache->hash);
98  }
99 
100  assert( rhs.fReferences.empty() ); // Copying with pending loads is unexpected
101  }
102 
103  //----------------------------------------------------------------------
105  {
106  if(this == &rhs) return *this;
107 
109 
110  if(rhs.fCache->hash){
111  fCache->spect = rhs.fCache->spect;
112  fCache->hash = std::make_unique<TMD5>(*rhs.fCache->hash);
113  }
114  else{
115  fCache->hash.reset();
116  }
117 
118  assert( rhs.fReferences.empty() ); // Copying with pending loads is unexpected
119  assert( fReferences.empty() ); // Copying with pending loads is unexpected
120 
121  return *this;
122  }
123 
124  //----------------------------------------------------------------------
126  {
127  if(this == &rhs) return *this;
128 
130 
131  if(rhs.fCache->hash){
132  fCache->spect = std::move(rhs.fCache->spect);
133  fCache->hash = std::move(rhs.fCache->hash);
134  }
135  else{
136  fCache->hash.reset();
137  }
138 
139  assert( rhs.fReferences.empty() ); // Copying with pending loads is unexpected
140  assert( fReferences.empty() ); // Copying with pending loads is unexpected
141 
142  return *this;
143  }
144 
145  //----------------------------------------------------------------------
146  template<class T> Spectrum OscillatableSpectrum::
147  _Oscillated(osc::_IOscCalc<T>* calc, int from, int to) const
148  {
149  TMD5* hash = calc->GetParamsHash();
150  if(hash && fCache->hash && *hash == *fCache->hash){
151  delete hash;
152  return fCache->spect;
153  }
154 
155  const OscCurve curve(calc, from, to);
156  const Spectrum ret = WeightedBy(curve);
157  if(hash){
158  fCache->spect = ret;
159  fCache->hash.reset(hash);
160  }
161 
162  return ret;
163  }
164 
165  //----------------------------------------------------------------------
167  int from, int to) const
168  {
169  return _Oscillated(calc, from, to);
170  }
171 
172  //----------------------------------------------------------------------
174  int from, int to) const
175  {
176  return _Oscillated(calc, from, to);
177  }
178 
179  //----------------------------------------------------------------------
181  {
183 
184  // invalidate
185  fCache->hash.reset(nullptr);
186 
187  return *this;
188  }
189 
190  //----------------------------------------------------------------------
192  {
193  OscillatableSpectrum ret = *this;
194  ret += rhs;
195  return ret;
196  }
197 
198  //----------------------------------------------------------------------
200  {
202 
203  // invalidate
204  fCache->hash.reset(nullptr);
205 
206  return *this;
207  }
208 
209  //----------------------------------------------------------------------
211  {
212  OscillatableSpectrum ret = *this;
213  ret -= rhs;
214  return ret;
215  }
216 
217  //----------------------------------------------------------------------
218  void OscillatableSpectrum::SaveTo(TDirectory* dir, const std::string& name) const
219  {
220  _SaveTo(dir, name, "OscillatableSpectrum");
221  }
222 
223  //----------------------------------------------------------------------
224  std::unique_ptr<OscillatableSpectrum> OscillatableSpectrum::LoadFrom(TDirectory* dir, const std::string& name)
225  {
226  dir = dir->GetDirectory(name.c_str()); // switch to subdir
227  assert(dir);
228 
229  DontAddDirectory guard;
230 
231  TObjString* tag = (TObjString*)dir->Get("type");
232  assert(tag);
233  assert(tag->GetString() == "OscillatableSpectrum");
234  delete tag;
235 
236  TH2D* spect = (TH2D*)dir->Get("hist");
237  assert(spect);
238  TH1* hPot = (TH1*)dir->Get("pot");
239  assert(hPot);
240  TH1* hLivetime = (TH1*)dir->Get("livetime");
241  assert(hLivetime);
242 
243  std::vector<std::string> labels;
244  std::vector<Binning> bins;
245 
246  for(int i = 0; ; ++i){
247  const std::string subname = TString::Format("bins%d", i).Data();
248  TDirectory* subdir = dir->GetDirectory(subname.c_str());
249  if(!subdir) break;
250  delete subdir;
251  bins.push_back(*Binning::LoadFrom(dir, subname));
252  TObjString* label = (TObjString*)dir->Get(TString::Format("label%d", i));
253  labels.push_back(label ? label->GetString().Data() : "");
254  delete label;
255  }
256 
257  delete dir;
258 
259  auto ret = std::make_unique<OscillatableSpectrum>(kNullLoader,
260  HistAxis(labels, bins),
261  kNoCut);
262 
263  // ROOT histogram storage is row-major, but Eigen is column-major by
264  // default
265  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen:: Dynamic, Eigen::RowMajor> MatRowMajor;
266  ret->fMat = Eigen::Map<MatRowMajor>(spect->GetArray(),
267  ret->fMat.rows(),
268  ret->fMat.cols());
269 
270  delete spect;
271 
272  ret->fPOT = hPot->Integral(0, -1);
273  ret->fLivetime = hLivetime->Integral(0, -1);
274 
275  delete hPot;
276  delete hLivetime;
277 
278  return ret;
279  }
280 }
static std::unique_ptr< OscillatableSpectrum > LoadFrom(TDirectory *dir, const std::string &name)
static std::unique_ptr< Binning > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Binning.cxx:230
const XML_Char * name
Definition: expat.h:151
ReweightableSpectrum & operator+=(const ReweightableSpectrum &rhs)
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
subdir
Definition: cvnie.py:7
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
General interface to oscillation calculators.
Definition: StanTypedefs.h:23
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
ThreadLocal< OscCache > fCache
Spectrum with the value of a second variable, allowing for reweighting
var_value< double > var
Definition: StanTypedefs.h:14
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
OscillatableSpectrum()
Constructor for Uninitialized()
osc::OscCalcDumb calc
OscillatableSpectrum & operator-=(const OscillatableSpectrum &rhs)
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
Definition: Vars.cxx:85
ReweightableSpectrum & operator-=(const ReweightableSpectrum &rhs)
Spectrum WeightedBy(const Ratio &weights) const
Reco spectrum with truth weights applied.
#define pot
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
caf::StandardRecord * sr
Float_t mat
Definition: plot.C:39
const Binning kTrueEnergyBins
Default true-energy bin edges.
Definition: Binning.h:14
OscillatableSpectrum & operator+=(const OscillatableSpectrum &rhs)
loader
Definition: demo0.py:10
OscillatableSpectrum operator+(const OscillatableSpectrum &rhs) const
virtual TMD5 * GetParamsHash() const
Use to check two calculators are in the same state.
Definition: IOscCalc.h:39
Base class for the various types of spectrum loader.
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::set< ReweightableSpectrum ** > fReferences
const Cut cut
Definition: exporter_fd.C:30
Transition probability for any one channel as a function of energy.
Definition: OscCurve.h:11
void _SaveTo(TDirectory *dir, const std::string &name, const std::string &type) const
Spectrum Oscillated(osc::IOscCalc *calc, int from, int to) const
ReweightableSpectrum & operator=(const ReweightableSpectrum &rhs)
const Cut kHasNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return true;})
Definition: TruthCuts.h:56
TDirectory * dir
Definition: macro.C:5
double livetime
Definition: saveFDMCHists.C:21
assert(nhit_max >=nhit_nbins)
Spectrum _Oscillated(osc::_IOscCalc< T > *calc, int from, int to) const
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
OscillatableSpectrum & operator=(const OscillatableSpectrum &rhs)
Assignment operator.
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Spectrum with true energy information, allowing it to be oscillated
OscillatableSpectrum operator-(const OscillatableSpectrum &rhs) const
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
void SaveTo(TDirectory *dir, const std::string &name) const
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
enum BeamMode string