Hist.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Hist.h"
2 
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Core/UtilsExt.h"
5 
6 #include "CAFAna/Core/Stan.h"
7 
8 #include "TH1.h"
9 #include "THnSparse.h"
10 
11 #include "TDirectory.h"
12 
13 namespace
14 {
15  namespace util{template<class T> T sqr(const T& x){return x*x;}}
16 }
17 
19 {
20 public:
22  {
23  fEnabled = (getenv("CAFANA_STAT_ERRS") != 0);
24  if(fEnabled) std::cout << "Propagation of statistical uncertainties enabled" << std::endl;
25  }
26 
27  operator bool() const {return fEnabled;}
28 
29 protected:
30  bool fEnabled;
31 } gStatErrs;
32 
33 namespace ana
34 {
35  //----------------------------------------------------------------------
36  Hist::Hist() : fType(kUninitialized), fSqrtErrs(false)
37  {
38  }
39 
40  //----------------------------------------------------------------------
42  {
43  Hist ret;
44  ret.fType = kDense;
45  ret.fData = Eigen::ArrayXd::Zero(nbins+2);
46  ret.fSqrtErrs = true;
47  return ret;
48  }
49 
50  //----------------------------------------------------------------------
52  {
53  Hist ret;
54  ret.fType = kSparse;
55  ret.fDataSparse = Eigen::SparseVector<double>(nbins+2);
56  ret.fSqrtErrs = true;
57  return ret;
58  }
59 
60  //----------------------------------------------------------------------
61  Hist::Hist(const Hist& rhs)
62  : Hist()
63  {
64  assert(rhs.Initialized());
65 
66  fType = rhs.fType;
67  // Only one of these will actually have contents
69  fDataStan = rhs.fDataStan;
70  fData = rhs.fData;
71  fSumSq = rhs.fSumSq;
72  fSqrtErrs = rhs.fSqrtErrs;
73  }
74 
75  //----------------------------------------------------------------------
77  : Hist()
78  {
79  assert(rhs.Initialized());
80 
81  fType = rhs.fType;
82  std::swap(fDataSparse, rhs.fDataSparse);
83  std::swap(fDataStan, rhs.fDataStan);
84  std::swap(fData, rhs.fData);
85  std::swap(fSumSq, rhs.fSumSq);
86  fSqrtErrs = rhs.fSqrtErrs;
87  }
88 
89  //----------------------------------------------------------------------
90  Hist& Hist::operator=(const Hist& rhs)
91  {
92  if(this == &rhs) return *this;
93 
94  DontAddDirectory guard;
95 
96  assert(rhs.Initialized());
97 
98  fType = rhs.fType;
100  fDataStan = rhs.fDataStan;
101  fData = rhs.fData;
102  fSumSq = rhs.fSumSq;
103  fSqrtErrs = rhs.fSqrtErrs;
104 
105  return *this;
106  }
107 
108  //----------------------------------------------------------------------
110  {
111  if(this == &rhs) return *this;
112 
113  assert(rhs.Initialized());
114 
115  fType = rhs.fType;
116  std::swap(fDataSparse, rhs.fDataSparse);
117  std::swap(fDataStan, rhs.fDataStan);
118  std::swap(fData, rhs.fData);
119  std::swap(fSumSq, rhs.fSumSq);
120  fSqrtErrs = rhs.fSqrtErrs;
121 
122  return *this;
123  }
124 
125  //----------------------------------------------------------------------
127  {
128  }
129 
130  //----------------------------------------------------------------------
131  Hist Hist::AdoptSparse(Eigen::SparseVector<double>&& v)
132  {
133  Hist ret;
134  ret.fType = kSparse;
135  ret.fDataSparse = std::move(v);
136  return ret;
137  }
138 
139  //----------------------------------------------------------------------
141  {
142  Hist ret;
143  ret.fType = kDenseStan;
144  ret.fDataStan = std::move(v);
145  return ret;
146  }
147 
148  //----------------------------------------------------------------------
149  Hist Hist::Adopt(Eigen::ArrayXd&& v)
150  {
151  Hist ret;
152  ret.fType = kDense;
153  ret.fData = std::move(v);
154  return ret;
155  }
156 
157  //----------------------------------------------------------------------
159  {
160  Hist ret;
161 
162  TH1D* h = (TH1D*)dir->Get("hist");
163  THnSparseD* hSparse = (THnSparseD*)dir->Get("hist_sparse");
164  assert(bool(h) != bool(hSparse));
165 
166  if(h){
167  ret.fType = kDense;
168  ret.fData = Eigen::Map<Eigen::ArrayXd>(h->GetArray(), h->GetNbinsX()+2);
169 
170  if(gStatErrs){
171  const int N = h->GetNbinsX()+2;
172  ret.fSumSq.resize(N);
173  for(int i = 0; i < N; ++i) ret.fSumSq[i] = util::sqr(h->GetBinError(i));
174  }
175  }
176  if(hSparse){
177  ret.fType = kSparse;
178 
179  // Number of bins there would be in a dense version
180  int nbins = 1;
181  const int N = hSparse->GetNdimensions();
182  for(int d = 0; d < N; ++d)
183  nbins *= hSparse->GetAxis(d)->GetNbins()+2;
184 
185  ret.fDataSparse.resize(nbins);
186  ret.fDataSparse.setZero();
187  for(int i = 0; i < hSparse->GetNbins(); ++i){
188  int idx[N];
189  const double y = hSparse->GetBinContent(i, idx);
190  int idx1d = idx[0];
191  int mult = 1;
192  for(int d = 1; d < N; ++d){
193  mult *= hSparse->GetAxis(d-1)->GetNbins()+2;
194  idx1d += mult*idx[d];
195  }
196  ret.fDataSparse.coeffRef(idx1d) = y;
197  }
198  }
199 
200  delete h;
201  delete hSparse;
202 
203  return ret;
204  }
205 
206  //----------------------------------------------------------------------
207  TH1D* Hist::ToTH1(const Binning& bins) const
208  {
209  assert(Initialized());
210 
211  TH1D* ret = MakeTH1D(UniqueName().c_str(), "", bins);
212 
213  for(int i = 0; i < bins.NBins()+2; ++i){
214  switch(fType){
215  case kDense: ret->SetBinContent(i, fData[i]); break;
216  case kDenseStan: ret->SetBinContent(i, fDataStan[i].val()); break;
217  // Interface requires returning literally a TH1D in any case
218  case kSparse: ret->SetBinContent(i, fDataSparse.coeff(i)); break;
219  default:
220  abort(); // unreachable
221  }
222  if(fSumSq.size() > 0 || fSqrtErrs) ret->SetBinError(i, GetBinError(i));
223  }
224  return ret;
225  }
226 
227  //----------------------------------------------------------------------
228  int Hist::GetNbinsX() const
229  {
230  assert(Initialized());
231 
232  switch(fType){
233  case kSparse: return fDataSparse.size()-2;
234  case kDenseStan: return fDataStan .size()-2;
235  case kDense: return fData .size()-2;
236  default: abort(); // unreachable
237  }
238  }
239 
240  //----------------------------------------------------------------------
241  double Hist::GetBinError(int i) const
242  {
243  assert(Initialized());
244 
245  if(fSqrtErrs) return sqrt(GetBinContent(i));
246  if(fSumSq.size() > 0) return sqrt(fSumSq[i]);
247 
248  return 0;
249  }
250 
251  //----------------------------------------------------------------------
252  double Hist::Integral() const
253  {
254  assert(Initialized());
255 
256  switch(fType){
257  case kSparse: return fDataSparse.sum();
258  case kDenseStan: return fDataStan .sum().val();
259  case kDense: return fData .sum();
260  default: abort(); // unreachable
261  }
262  }
263 
264  //----------------------------------------------------------------------
265  void Hist::Fill(const Binning& bins, double x, double w)
266  {
267  assert(Initialized());
268 
269  if(w != 1) fSqrtErrs = false;
270 
271  if(fType == kDense && gStatErrs && fSumSq.size() == 0) fSumSq.resize(fData.size());
272 
273  switch(fType){
274  case kSparse:
275  fDataSparse.coeffRef(bins.FindBin(x)) += w;
276  break;
277  case kDenseStan:
278  std::cout << "Hist::Fill() not supported for stan vars" << std::endl;
279  abort();
280  case kDense:
281  {
282  const int bin = bins.FindBin(x);
283  fData[bin] += w;
284  if(gStatErrs){
285  if(fSumSq.size() == 0) fSumSq.resize(fData.size());
286  fSumSq[bin] += w*w;
287  }
288  }
289  break;
290  default:
291  abort(); // unreachable
292  }
293  }
294 
295  //----------------------------------------------------------------------
296  void Hist::Scale(double s)
297  {
298  assert(Initialized());
299 
300  switch(fType){
301  case kSparse: fDataSparse *= s; break;
302  case kDenseStan: fDataStan *= s; break;
303  case kDense: fData *= s; break;
304  default:
305  abort(); // unreachable
306  }
307 
308  if(s != 1) fSqrtErrs = false;
309  fSumSq *= s;
310  }
311 
312  //----------------------------------------------------------------------
314  {
315  assert(Initialized());
316  switch (fType){
317  case kSparse:
318  fDataStan = fDataSparse * s;
319  fDataSparse.resize(0);
320  break;
321 
322  case kDense:
323  fDataStan = fData * s;
324  fData.resize(0);
325  break;
326 
327  case kDenseStan:
328  fDataStan *= s;
329  break;
330 
331  default:
332  abort(); // unreachable
333  }
334 
335  fType = kDenseStan;
336 
337  if(s != 1) fSqrtErrs = false;
338  fSumSq *= s.val();
339  }
340 
341  //----------------------------------------------------------------------
343  {
344  fSumSq.resize(0);
345  fSqrtErrs = true;
346  }
347 
348  //----------------------------------------------------------------------
349  double Hist::GetBinContent(int i) const
350  {
351  assert(Initialized());
352 
353  switch(fType){
354  case kSparse: return fDataSparse.coeff(i);
355  case kDenseStan: return fDataStan[i].val();
356  case kDense: return fData[i];
357  default:
358  abort(); // unreachable
359  }
360  }
361 
362  //----------------------------------------------------------------------
363  void Hist::SetBinContent(int i, double x)
364  {
365  assert(Initialized());
366 
367  switch(fType){
368  case kSparse: fDataSparse.coeffRef(i) = x; break;
369  case kDenseStan:
370  std::cout << "Hist::SetBinContent() not implemented for stan vars" << std::endl;
371  abort();
372  case kDense: fData[i] = x; break;
373  default:
374  abort(); // unreachable
375  }
376 
377  fSqrtErrs = false;
378  if(fSumSq.size() > 0) fSumSq[i] = 0;
379  }
380 
381  //----------------------------------------------------------------------
382  void Hist::Reset()
383  {
384  switch(fType){
385  case kSparse: fDataSparse.setZero(); break;
386  case kDenseStan: fDataStan .setZero(); break;
387  case kDense: fData .setZero(); break;
388  default: ; // OK?
389  }
390 
391  fSqrtErrs = true;
392  fSumSq.resize(0);
393  }
394 
395  //----------------------------------------------------------------------
396  void Hist::Add(const Eigen::SparseVector<double>& rhs, double scale)
397  {
398  switch(fType){
399  case kSparse: fDataSparse += rhs * scale; break;
400  case kDenseStan: fDataStan += rhs * scale; break;
401  case kDense: fData += rhs * scale; break;
402  default: abort(); // unreachable
403  }
404  }
405 
406  //----------------------------------------------------------------------
407  void Hist::Add(const Eigen::ArrayXstan& rhs, double scale)
408  {
409  switch(fType){
410  case kSparse:
411  fType = kDenseStan;
412  fDataStan = rhs * scale;
414  fDataSparse.resize(0);
415  break;
416 
417  case kDenseStan:
418  fDataStan += rhs * scale;
419  break;
420 
421  case kDense:
422  fType = kDenseStan;
423  fDataStan = fData + rhs * scale;
424  fData.resize(0);
425  break;
426 
427  default: abort(); // unreachable
428  }
429  }
430 
431  //----------------------------------------------------------------------
432  void Hist::Add(const Eigen::ArrayXd& rhs, double scale)
433  {
434  switch(fType){
435  case kSparse:
436  fType = kDense;
437  fData = rhs * scale;
438  fData += fDataSparse;
439  fDataSparse.resize(0);
440  break;
441 
442  case kDenseStan:
443  fDataStan += rhs * scale;
444  break;
445 
446  case kDense:
447  fData += rhs * scale;
448  break;
449 
450  default: abort(); // unreachable
451  }
452  }
453 
454  //----------------------------------------------------------------------
455  void Hist::Add(const Hist& rhs, double scale)
456  {
457  assert(Initialized());
458  assert(rhs.Initialized());
459 
460  switch(rhs.fType){
461  case kSparse: Add(rhs.fDataSparse, scale); break;
462  case kDenseStan: Add(rhs.fDataStan, scale); break;
463  case kDense: Add(rhs.fData, scale); break;
464  default: abort(); // unreachable
465  }
466 
467  if(scale != 1 || !rhs.fSqrtErrs) fSqrtErrs = false;
468 
469  if(fSumSq.size() > 0){
470  if(rhs.fSumSq.size() > 0) fSumSq += scale * rhs.fSumSq;
471  // otherwise nothing to add in
472  }
473  else{
474  fSumSq = scale * rhs.fSumSq;
475  }
476  }
477 
478  //----------------------------------------------------------------------
479  void Hist::Multiply(const Hist& rhs)
480  {
481  assert(Initialized());
482  assert(rhs.Initialized());
483 
484  if(fType == kSparse || rhs.fType == kSparse){
485  std::cout << "Hist::Multiply() not implemented for sparse vectors" << std::endl;
486  abort();
487  }
488 
489  if(fType == kDenseStan){
490  if(rhs.fType == kDenseStan){
491  fDataStan *= rhs.fDataStan;
492  }
493  else{
494  fDataStan *= rhs.fData;
495  }
496  }
497  else{
498  if(rhs.fType == kDenseStan){
499  fType = kDenseStan;
500  fDataStan = fData * rhs.fDataStan;
501  fData.resize(0);
502  }
503  else{
504  fData *= rhs.fData;
505  }
506  }
507 
508  if(fType == kDense && rhs.fType == kDense){
509  if(fSumSq.size() > 0){
510  if(rhs.fSumSq.size() > 0){
511  fSumSq = fSumSq * util::sqr(rhs.fData) + rhs.fSumSq * util::sqr(fData);
512  }
513  else{
514  fSumSq *= util::sqr(rhs.fData);
515  }
516  }
517  else if(rhs.fSumSq.size() > 0){
518  fSumSq = rhs.fSumSq * util::sqr(fData);
519  }
520  }
521  else{
522  // Didn't bother to implement error prop for multiplying stan/sparse
523  // hists
524  fSumSq.resize(0);
525  }
526 
527  fSqrtErrs = false;
528  }
529 
530  //----------------------------------------------------------------------
531  void Hist::Divide(const Hist& rhs)
532  {
533  assert(Initialized());
534  assert(rhs.Initialized());
535 
536  if(fType == kSparse || rhs.fType == kSparse){
537  std::cout << "Hist::Divide() not implemented for sparse vectors" << std::endl;
538  abort();
539  }
540 
541  if(fType == kDenseStan){
542  if(rhs.fType == kDenseStan){
543  fDataStan /= rhs.fDataStan;
544  }
545  else{
546  fDataStan /= rhs.fData;
547  }
548  }
549  else{
550  if(rhs.fType == kDenseStan){
551  fType = kDenseStan;
552  fDataStan = fData / rhs.fDataStan;
553  fData.resize(0);
554  }
555  else{
556  fData /= rhs.fData;
557  }
558  }
559 
560  if(fType == kDense && rhs.fType == kDense){
561  if(fSumSq.size() > 0){
562  if(rhs.fSumSq.size() > 0){
564  }
565  else{
566  fSumSq /= util::sqr(rhs.fData);
567  }
568  }
569  else if(rhs.fSumSq.size() > 0){
571  }
572  }
573  else{
574  // Didn't bother to implement error prop for dividing stan/sparse
575  // hists
576  fSumSq.resize(0);
577  }
578 
579  fSqrtErrs = false;
580  }
581 
582  //----------------------------------------------------------------------
583  void Hist::Write(const Binning& bins) const
584  {
585  assert(Initialized());
586 
587  if(fType == kDenseStan){
588  std::cout << "Hist::Write() not implemented (impossible?) for stan vars" << std::endl;
589  abort();
590  }
591  if(fType == kDense){
592  TH1D* h = ToTH1(bins);
593  h->Write("hist");
594  delete h;
595  }
596  if(fType == kSparse){
597  const int n = bins.NBins();
598  const double x0 = bins.IsSimple() ? bins.Min() : 0;
599  const double x1 = bins.IsSimple() ? bins.Max() : bins.NBins();
600  THnSparseD* h = new THnSparseD("", "", 1, &n, &x0, &x1);
601 
602  for(Eigen::SparseVector<double>::InnerIterator it(fDataSparse); it; ++it){
603  const int idx = it.index();
604  h->SetBinContent(&idx, it.value());
605  if(fSqrtErrs) h->SetBinError(&idx, sqrt(it.value()));
606  }
607 
608  h->Write("hist_sparse");
609  delete h;
610  }
611  }
612 }
class StatErrorsEnabled gStatErrs
~Hist()
Definition: Hist.cxx:126
void SetBinContent(int i, double x)
Definition: Hist.cxx:363
Eigen::ArrayXd fSumSq
Accumulate errors, if enabled.
Definition: Hist.h:89
Hist()
Definition: Hist.cxx:36
void Write(const Binning &bins) const
Definition: Hist.cxx:583
Filter events based on their run/event numbers.
void Add(const Hist &rhs, double scale=1)
Definition: Hist.cxx:455
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
set< int >::iterator it
double Integral() const
Definition: Hist.cxx:252
void ResetErrors()
Definition: Hist.cxx:342
Float_t x1[n_points_granero]
Definition: compare.C:5
int FindBin(double x) const
Definition: Binning.cxx:180
TH1D * MakeTH1D(const char *name, const char *title, const Binning &bins)
Utility function to avoid need to switch on bins.IsSimple()
Definition: UtilsExt.cxx:74
T sqrt(T number)
Definition: d0nt_math.hpp:156
Definition: Hist.h:29
static Hist FromDirectory(TDirectory *dir)
Definition: Hist.cxx:158
int GetNbinsX() const
Definition: Hist.cxx:228
TH1D * ToTH1(const Binning &bins) const
Definition: Hist.cxx:207
Eigen::SparseVector< double > fDataSparse
Definition: Hist.h:86
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
bool IsSimple() const
Definition: Binning.h:33
bool fSqrtErrs
Special case when filled with unweighted data.
Definition: Hist.h:90
static Hist ZeroSparse(int nbins)
Definition: Hist.cxx:51
const XML_Char * s
Definition: expat.h:262
bool Initialized() const
Definition: Hist.h:39
const int nbins
Definition: cellShifts.C:15
Double_t scale
Definition: plot.C:25
double Min() const
Definition: Binning.h:30
static Hist AdoptSparse(Eigen::SparseVector< double > &&v)
Definition: Hist.cxx:131
std::string getenv(std::string const &name)
Float_t d
Definition: plot.C:236
void Scale(double s)
Definition: Hist.cxx:296
void Reset()
Definition: Hist.cxx:382
Eigen::ArrayXstan fDataStan
Definition: Hist.h:87
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
void Fill(const Binning &bins, double x, double w)
Definition: Hist.cxx:265
double Max() const
Definition: Binning.h:31
const Binning bins
Definition: NumuCC_CPiBin.h:8
static Hist AdoptStan(Eigen::ArrayXstan &&v)
Definition: Hist.cxx:140
Eigen::ArrayXd fData
Definition: Hist.h:88
EType fType
Definition: Hist.h:84
void Divide(const Hist &rhs)
Definition: Hist.cxx:531
Eigen::Array< stan::math::var, Eigen::Dynamic, 1 > ArrayXstan
Definition: StanUtils.h:7
TDirectory * dir
Definition: macro.C:5
void Zero()
int NBins() const
Definition: Binning.h:29
double GetBinContent(int i) const
Definition: Hist.cxx:349
static Hist Zero(int nbins)
Definition: Hist.cxx:41
static Hist Adopt(Eigen::ArrayXd &&v)
Definition: Hist.cxx:149
assert(nhit_max >=nhit_nbins)
double T
Definition: Xdiff_gwt.C:5
double val() const
Definition: var.hpp:294
Hist & operator=(const Hist &)
Definition: Hist.cxx:90
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
double GetBinError(int i) const
Definition: Hist.cxx:241
Float_t w
Definition: plot.C:20
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void Multiply(const Hist &rhs)
Definition: Hist.cxx:479