ReweightableSpectrum.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Core/Ratio.h"
5 #include "CAFAna/Core/Stan.h"
6 
7 #include "TDirectory.h"
8 #include "TH2.h"
9 #include "TObjString.h"
10 
11 #include <cassert>
12 #include <iostream>
13 #include <memory>
14 
15 namespace ana
16 {
17  //----------------------------------------------------------------------
19  const LabelsAndBins& recoAxis,
20  const LabelsAndBins& trueAxis,
21  double pot, double livetime)
22  : fMat(mat), fPOT(pot), fLivetime(livetime),
23  fAxisX(recoAxis), fAxisY(trueAxis)
24  {
25  }
26 
27  //----------------------------------------------------------------------
29  {
30  for(ReweightableSpectrum** ref: fReferences) *ref = 0;
31  }
32 
33  //----------------------------------------------------------------------
35  : fAxisX(rhs.fAxisX), fAxisY(rhs.fAxisY)
36  {
37  DontAddDirectory guard;
38 
39  fMat = rhs.fMat;
40  fPOT = rhs.fPOT;
41  fLivetime = rhs.fLivetime;
42 
43  assert( rhs.fReferences.empty() ); // Copying with pending loads is unexpected
44  }
45 
46  //----------------------------------------------------------------------
48  {
49  if(this == &rhs) return *this;
50 
51  DontAddDirectory guard;
52 
53  fAxisX = rhs.fAxisX;
54  fAxisY = rhs.fAxisY;
55 
56  fMat = rhs.fMat;
57  fPOT = rhs.fPOT;
58  fLivetime = rhs.fLivetime;
59 
60  assert( fReferences.empty() ); // Copying with pending loads is unexpected
61 
62  return *this;
63  }
64 
65  //----------------------------------------------------------------------
66  TH2D* ReweightableSpectrum::ToTH2(double pot) const
67  {
68  // Could have a file temporarily open
69  DontAddDirectory guard;
70 
71  TH2D* ret = MakeTH2D(UniqueName().c_str(), "", fAxisX.GetBins1D(), fAxisY.GetBins1D());
72 
73  for(int i = 0; i < fMat.rows(); ++i){
74  for(int j = 0; j < fMat.cols(); ++j){
75  ret->SetBinContent(j, i, fMat(i, j));
76  }
77  }
78 
79  if(fPOT){
80  ret->Scale(pot/fPOT);
81  }
82  else{
83  // How did it get events with no POT?
84  assert(ret->Integral() == 0);
85  }
86 
87  ret->GetXaxis()->SetTitle(fAxisX.GetLabel1D().c_str());
88  ret->GetYaxis()->SetTitle(fAxisY.GetLabel1D().c_str());
89 
90  return ret;
91  }
92 
93  //----------------------------------------------------------------------
94  void ReweightableSpectrum::Fill(double x, double y, double w)
95  {
97  }
98 
99  /// Helper for \ref Unweighted
100  inline Eigen::ArrayXd ProjectionX(const Eigen::MatrixXd& mat)
101  {
102  return Eigen::RowVectorXd::Ones(mat.rows()) * mat;
103  }
104 
105  /// Helper for \ref WeightingVariable
106  inline Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd& mat)
107  {
108  return mat * Eigen::VectorXd::Ones(mat.cols());
109  }
110 
111  //----------------------------------------------------------------------
113  {
115  }
116 
117  //----------------------------------------------------------------------
119  {
121  }
122 
123  //----------------------------------------------------------------------
125  {
126  if(!ws.HasStan()){
127  const Eigen::VectorXd& vec = ws.GetEigen();
128 
129  return Spectrum(Hist::Adopt(Eigen::ArrayXd(vec.transpose() * fMat)),
130  fAxisX, fPOT, fLivetime);
131  }
132  else{
133  const Eigen::VectorXstan& vec = ws.GetEigenStan();
134 
135  return Spectrum(Hist::AdoptStan(vec.transpose() * fMat),
136  fAxisX, fPOT, fLivetime);
137  }
138  }
139 
140  //----------------------------------------------------------------------
142  {
143  // This is a big component of what extrapolations do, so it should be fast
144 
145  const Ratio ratio(target, WeightingVariable());
146  // We want to multiply all the rows by this ratio, so left-multiply
147 
148  fMat = ratio.GetEigen().matrix().asDiagonal() * fMat;
149  }
150 
151  //----------------------------------------------------------------------
153  {
154  // This is a big component of what extrapolations do, so it should be fast
155 
156  const Ratio ratio(target, UnWeighted());
157  // We want to multiply all the columns by this ratio
158 
159  fMat *= ratio.GetEigen().matrix().asDiagonal();
160  }
161 
162  //----------------------------------------------------------------------
164  {
165  // In this case it would be OK to have no POT/livetime
166  if(rhs.fMat.sum() == 0) return *this;
167 
168 
169  if((!fPOT && !fLivetime) || (!rhs.fPOT && !rhs.fLivetime)){
170  std::cout << "Error: can't sum ReweightableSpectrum with no POT or livetime."
171  << fPOT << " " << rhs.fPOT
172  << std::endl;
173  // abort();
174  return *this;
175  }
176 
177  if(!fLivetime && !rhs.fPOT){
178  std::cout << "Error: can't sum ReweightableSpectrum with POT ("
179  << fPOT << ") but no livetime and ReweightableSpectrum with livetime ("
180  << rhs.fLivetime << " sec) but no POT." << std::endl;
181  abort();
182  }
183 
184  if(!fPOT && !rhs.fLivetime){
185  std::cout << "Error: can't sum ReweightableSpectrum with livetime ("
186  << fLivetime << " sec) but no POT and ReweightableSpectrum with POT ("
187  << rhs.fPOT << ") but no livetime." << std::endl;
188  abort();
189  }
190 
191  // And now there are still a bunch of good cases to consider
192 
193  if(fPOT && rhs.fPOT){
194  // Scale by POT when possible
195  fMat += rhs.fMat * sign*fPOT/rhs.fPOT;
196 
197  if(fLivetime && rhs.fLivetime){
198  // If POT/livetime ratios match, keep regular lifetime, otherwise zero
199  // it out.
200  if(AlmostEqual(fLivetime*rhs.fPOT, rhs.fLivetime*fPOT))
201  fLivetime = 0;
202  }
203  if(!fLivetime && rhs.fLivetime){
204  // If the RHS has a livetime and we don't, copy it in (suitably scaled)
205  fLivetime = rhs.fLivetime * fPOT/rhs.fPOT;
206  }
207  // Otherwise, keep our own livetime (if any)
208 
209  return *this;
210  }
211 
212  if(fLivetime && rhs.fLivetime){
213  // Scale by livetime, the only thing in common
214  fMat += rhs.fMat * sign*fLivetime/rhs.fLivetime;
215 
216  if(!fPOT && rhs.fPOT){
217  // If the RHS has a POT and we don't, copy it in (suitably scaled)
218  fPOT = rhs.fPOT * fLivetime/rhs.fLivetime;
219  }
220  // Otherwise, keep our own POT (if any)
221 
222  return *this;
223  }
224 
225  // That should have been all the cases. I definitely want to know what
226  // happened if it wasn't.
227  std::cout << "ReweightableSpectrum::operator+=(). How did we get here? "
228  << fPOT << " " << fLivetime << " "
229  << rhs.fPOT << " " << rhs.fLivetime << std::endl;
230  abort();
231  }
232 
233  //----------------------------------------------------------------------
235  {
236  return PlusEqualsHelper(rhs, +1);
237  }
238 
239  //----------------------------------------------------------------------
241  {
242  ReweightableSpectrum ret = *this;
243  ret += rhs;
244  return ret;
245  }
246 
247  //----------------------------------------------------------------------
249  {
250  return PlusEqualsHelper(rhs, -1);
251  }
252 
253  //----------------------------------------------------------------------
255  {
256  ReweightableSpectrum ret = *this;
257  ret -= rhs;
258  return ret;
259  }
260 
261  //----------------------------------------------------------------------
263  {
264  fMat.setZero();
265  }
266 
267  //----------------------------------------------------------------------
269  {
270  fReferences.erase(ref);
271  }
272 
273  //----------------------------------------------------------------------
275  {
276  fReferences.insert(ref);
277  }
278 
279  //----------------------------------------------------------------------
280  void ReweightableSpectrum::SaveTo(TDirectory* dir, const std::string& name) const
281  {
282  _SaveTo(dir, name, "ReweightableSpectrum");
283  }
284 
285  //----------------------------------------------------------------------
287  const std::string& name,
288  const std::string& type) const
289  {
290  TDirectory* tmp = gDirectory;
291 
292  dir = dir->mkdir(name.c_str()); // switch to sbudir
293  dir->cd();
294 
295  TObjString(type.c_str()).Write("type");
296 
297  TH2* h = ToTH2(fPOT);
298  h->Write("hist");
299 
300  TH1D hPot("", "", 1, 0, 1);
301  hPot.Fill(.5, fPOT);
302  hPot.Write("pot");
303  TH1D hLivetime("", "", 1, 0, 1);
304  hLivetime.Fill(.5, fLivetime);
305  hLivetime.Write("livetime");
306 
307  for(unsigned int i = 0; i < fAxisX.NDimensions(); ++i){
308  TObjString(fAxisX.GetLabels()[i].c_str()).Write(TString::Format("label%d", i).Data());
309  fAxisX.GetBinnings()[i].SaveTo(dir, TString::Format("bins%d", i).Data());
310  }
311 
312  for(unsigned int i = 0; i < fAxisY.NDimensions(); ++i){
313  TObjString(fAxisY.GetLabels()[i].c_str()).Write(TString::Format("labely%d", i).Data());
314  fAxisY.GetBinnings()[i].SaveTo(dir, TString::Format("binsy%d", i).Data());
315  }
316 
317  dir->Write();
318  delete dir;
319 
320  delete h;
321 
322  tmp->cd();
323  }
324 
325  //----------------------------------------------------------------------
326  std::unique_ptr<ReweightableSpectrum> ReweightableSpectrum::LoadFrom(TDirectory* dir, const std::string& name)
327  {
328  dir = dir->GetDirectory(name.c_str()); // switch to subdir
329  assert(dir);
330 
331  TObjString* tag = (TObjString*)dir->Get("type");
332  assert(tag);
333  assert(tag->GetString() == "ReweightableSpectrum");
334  delete tag;
335 
336  TH2D* spect = (TH2D*)dir->Get("hist");
337  assert(spect);
338  TH1* hPot = (TH1*)dir->Get("pot");
339  assert(hPot);
340  TH1* hLivetime = (TH1*)dir->Get("livetime");
341  assert(hLivetime);
342 
343  std::vector<std::string> labels, labelsy;
344  std::vector<Binning> bins, binsy;
345 
346  for(int i = 0; ; ++i){
347  const std::string subname = TString::Format("bins%d", i).Data();
348  TDirectory* subdir = dir->GetDirectory(subname.c_str());
349  if(!subdir) break;
350  delete subdir;
351  bins.push_back(*Binning::LoadFrom(dir, subname.c_str()));
352  TObjString* label = (TObjString*)dir->Get(TString::Format("label%d", i));
353  labels.push_back(label ? label->GetString().Data() : "");
354  }
355 
356  for(int i = 0; ; ++i){
357  const std::string subname = TString::Format("binsy%d", i).Data();
358  TDirectory* subdir = dir->GetDirectory(subname.c_str());
359  if(!subdir) break;
360  delete subdir;
361  binsy.push_back(*Binning::LoadFrom(dir, subname.c_str()));
362  TObjString* labely = (TObjString*)dir->Get(TString::Format("labely%d", i));
363  labelsy.push_back(labely ? labely->GetString().Data() : "");
364  }
365 
366  // Backwards compatibility
367  if(labelsy.empty()) labelsy.push_back(spect->GetYaxis()->GetTitle());
368  if(binsy.empty()) binsy.push_back(Binning::FromTAxis(spect->GetYaxis()));
369 
370  const LabelsAndBins xax(labels, bins);
371  const LabelsAndBins yax(labelsy, binsy);
372 
373  // ROOT histogram storage is row-major, but Eigen is column-major by
374  // default
375  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen:: Dynamic, Eigen::RowMajor> MatRowMajor;
376  auto ret = std::make_unique<ReweightableSpectrum>(
377  Eigen::Map<MatRowMajor>(spect->GetArray(),
378  yax.GetBins1D().NBins()+2,
379  xax.GetBins1D().NBins()+2),
380  xax, yax, hPot->Integral(0, -1), hLivetime->Integral(0, -1));
381 
382  delete spect;
383 
384  delete hPot;
385  delete hLivetime;
386 
387  delete dir;
388 
389  return ret;
390  }
391 }
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)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const XML_Char * target
Definition: expat.h:268
const Eigen::ArrayXd & GetEigen() const
Definition: Ratio.h:44
Eigen::ArrayXd ProjectionX(const Eigen::MatrixXd &mat)
Helper for Unweighted.
const SpillVar fPOT([](const caf::SRSpillProxy *spill){return spill->spillpot;})
void Fill(double x, double y, double w=1)
subdir
Definition: cvnie.py:7
static Binning FromTAxis(const TAxis *ax)
Definition: Binning.cxx:173
int FindBin(double x) const
Definition: Binning.cxx:155
Eigen::Matrix< stan::math::var, Eigen::Dynamic, 1 > VectorXstan
Definition: Hist.h:14
Spectrum with the value of a second variable, allowing for reweighting
std::vector< double > Spectrum
Definition: Constants.h:743
Float_t tmp
Definition: plot.C:36
static std::unique_ptr< ReweightableSpectrum > LoadFrom(TDirectory *dir, const std::string &name)
void ReweightToRecoSpectrum(const Spectrum &target)
Recale bins so that Unweighted will return target.
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
ReweightableSpectrum operator-(const ReweightableSpectrum &rhs) const
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
void RemoveLoader(ReweightableSpectrum **)
const std::string & GetLabel1D() const
bool AlmostEqual(double a, double b, double eps)
Definition: UtilsExt.cxx:40
ReweightableSpectrum & operator-=(const ReweightableSpectrum &rhs)
const Binning & GetBins1D() const
Appropriate binning and labelling for that 1D Var.
Spectrum WeightedBy(const Ratio &weights) const
Reco spectrum with truth weights applied.
#define pot
const Eigen::ArrayXstan & GetEigenStan() const
Definition: Ratio.h:45
const double j
Definition: BetheBloch.cxx:29
Float_t mat
Definition: plot.C:39
TFile * fMat
Definition: makeCovMatrix.C:18
Represent the ratio between two spectra.
Definition: Ratio.h:8
OStream cout
Definition: OStream.cxx:6
ReweightableSpectrum operator+(const ReweightableSpectrum &rhs) const
TH2D * MakeTH2D(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Utility function to avoid 4-way combinatorial explosion on the bin types.
Definition: UtilsExt.cxx:86
const std::vector< Binning > & GetBinnings() const
Definition: LabelsAndBins.h:69
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::set< ReweightableSpectrum ** > fReferences
unsigned int NDimensions() const
Definition: LabelsAndBins.h:66
static Hist AdoptStan(Eigen::ArrayXstan &&v)
Definition: Hist.cxx:140
void _SaveTo(TDirectory *dir, const std::string &name, const std::string &type) const
ReweightableSpectrum & PlusEqualsHelper(const ReweightableSpectrum &rhs, int sign)
ReweightableSpectrum & operator=(const ReweightableSpectrum &rhs)
TDirectory * dir
Definition: macro.C:5
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
double livetime
Definition: saveFDMCHists.C:21
int NBins() const
Definition: Binning.h:29
void SaveTo(TDirectory *dir, const std::string &name) const
static Hist Adopt(Eigen::ArrayXd &&v)
Definition: Hist.cxx:149
void AddLoader(ReweightableSpectrum **)
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
void ReweightToTrueSpectrum(const Spectrum &target)
Rescale bins so that WeightingVariable will return target.
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t w
Definition: plot.C:20
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:28
bool HasStan() const
Definition: Ratio.h:43
const std::vector< std::string > & GetLabels() const
Definition: LabelsAndBins.h:68
def sign(x)
Definition: canMan.py:197
TH2D * ToTH2(double pot) const
gm Write()
enum BeamMode string