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