ATLAS Offline Software
Flex2DHisto.icc
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "RVersion.h"
6 
7 //____________________________________________________________________
8 template <class T>
9 inline Flex2DHisto<T> * Flex2DHisto<T>::create( unsigned nbinsx, const double& xmin, const double& xmax,
10  unsigned nbinsy, const double& ymin, const double& ymax )
11 {
12  return new(LWPools::acquire(sizeof(Flex2DHisto<T>)+extraAllocSize(nbinsx,nbinsy))) Flex2DHisto<T>(nbinsx,xmin,xmax,nbinsy,ymin,ymax);
13 }
14 
15 //____________________________________________________________________
16 template <class T>
17 template <class TFloat>
18 inline Flex2DHisto<T> * Flex2DHisto<T>::create( unsigned nbinsx, const TFloat* xbins,
19  unsigned nbinsy, const double& ymin, const double& ymax )
20 {
21  return new(LWPools::acquire(sizeof(Flex2DHisto<T>)+extraAllocSize(nbinsx,nbinsy))) Flex2DHisto<T>(nbinsx,xbins,nbinsy,ymin,ymax);
22 }
23 
24 //____________________________________________________________________
25 template <class T>
26 template <class TFloat>
27 inline Flex2DHisto<T> * Flex2DHisto<T>::create( unsigned nbinsx, const double& xmin, const double& xmax,
28  unsigned nbinsy, const TFloat* ybins )
29 {
30  return new(LWPools::acquire(sizeof(Flex2DHisto<T>)+extraAllocSize(nbinsx,nbinsy))) Flex2DHisto<T>(nbinsx,xmin,xmax,nbinsy,ybins);
31 }
32 
33 //____________________________________________________________________
34 template <class T>
35 template <class TFloat>
36 inline Flex2DHisto<T> * Flex2DHisto<T>::create( unsigned nbinsx, const TFloat*xbins,
37  unsigned nbinsy, const TFloat*ybins )
38 {
39  return new(LWPools::acquire(sizeof(Flex2DHisto<T>)+extraAllocSize(nbinsx,nbinsy))) Flex2DHisto<T>(nbinsx,xbins,nbinsy,ybins);
40 }
41 
42 //____________________________________________________________________
43 template <class T>
44 inline void Flex2DHisto<T>::destroy(Flex2DHisto<T> *h)
45 {
46  if (h) {
47  unsigned nx = h->getNBinsX();
48  unsigned ny = h->getNBinsY();
49  h->~Flex2DHisto<T>();
50  LWPools::release(reinterpret_cast<char*>(h),sizeof(Flex2DHisto<T>)+extraAllocSize(nx,ny));
51  }
52 }
53 
54 //____________________________________________________________________
55 template <class T>
56 inline Flex2DHisto<T>::Flex2DHisto( unsigned nbinsx, const double& xmin, const double& xmax,
57  unsigned nbinsy, const double& ymin, const double& ymax )
58  : m_invDeltaX(nbinsx/(xmax-xmin)),
59  m_invDeltaY(nbinsy/(ymax-ymin)),
60  m_sumW(0),
61  m_sumW2(0),
62  m_sumWX(0),
63  m_sumWX2(0),
64  m_sumWY(0),
65  m_sumWY2(0),
66  m_sumWXY(0),
67  m_nEntries(0),
68  m_nbinsxPlus1(nbinsx+1),
69  m_nbinsyPlus1(nbinsy+1),
70  m_nbinsxPlus2(nbinsx+2),
71  m_xmin(xmin),
72  m_xmax(xmax),
73  m_ymin(ymin),
74  m_ymax(ymax),
75  m_varBinsX(0),
76  m_varBinsY(0),
77  m_flexArray((nbinsx+2)*(nbinsy+2))
78 {
79  assert(xmin<xmax);
80  assert(ymin<ymax);
81  assert(nbinsx>0);
82  assert(nbinsy>0);
83  assert(getNBinsX()==nbinsx);
84  assert(getNBinsY()==nbinsy);
85 }
86 
87 //____________________________________________________________________
88 template <class T>
89 template <class TFloat>
90 inline Flex2DHisto<T>::Flex2DHisto( unsigned nbinsx, const TFloat*xbins,
91  unsigned nbinsy, const double& ymin, const double& ymax )
92  : m_invDeltaX(nbinsx/(xbins[nbinsx]-xbins[0])),
93  m_invDeltaY(nbinsy/(ymax-ymin)),
94  m_sumW(0),
95  m_sumW2(0),
96  m_sumWX(0),
97  m_sumWX2(0),
98  m_sumWY(0),
99  m_sumWY2(0),
100  m_sumWXY(0),
101  m_nEntries(0),
102  m_nbinsxPlus1(nbinsx+1),
103  m_nbinsyPlus1(nbinsy+1),
104  m_nbinsxPlus2(nbinsx+2),
105  m_xmin(xbins[0]),
106  m_xmax(xbins[nbinsx]),
107  m_ymin(ymin),
108  m_ymax(ymax),
109  m_varBinsX(LWPools::acquire<float>(nbinsx+1)),
110  m_varBinsY(0),
111  m_flexArray((nbinsx+2)*(nbinsy+2))
112 {
113  assert(m_xmin<m_xmax);
114  assert(m_ymin<m_ymax);
115  assert(nbinsx>0);
116  assert(nbinsy>0);
117 #ifndef NDEBUG
118  for (unsigned i = 0; i<nbinsx;++i)
119  assert(xbins[i]<xbins[i+1]);
120 #endif
121  for (unsigned i = 0; i<m_nbinsxPlus1;++i)
122  m_varBinsX[i]=xbins[i];
123 }
124 
125 //____________________________________________________________________
126 template <class T>
127 template <class TFloat>
128 inline Flex2DHisto<T>::Flex2DHisto( unsigned nbinsx, const double& xmin, const double& xmax,
129  unsigned nbinsy, const TFloat* ybins )
130  : m_invDeltaX(nbinsx/(xmax-xmin)),
131  m_invDeltaY(nbinsy/(ybins[nbinsy]-ybins[0])),
132  m_sumW(0),
133  m_sumW2(0),
134  m_sumWX(0),
135  m_sumWX2(0),
136  m_sumWY(0),
137  m_sumWY2(0),
138  m_sumWXY(0),
139  m_nEntries(0),
140  m_nbinsxPlus1(nbinsx+1),
141  m_nbinsyPlus1(nbinsy+1),
142  m_nbinsxPlus2(nbinsx+2),
143  m_xmin(xmin),
144  m_xmax(xmax),
145  m_ymin(ybins[0]),
146  m_ymax(ybins[nbinsy]),
147  m_varBinsX(0),
148  m_varBinsY(LWPools::acquire<float>(nbinsy+1)),
149  m_flexArray((nbinsx+2)*(nbinsy+2))
150 {
151  assert(m_xmin<m_xmax);
152  assert(m_ymin<m_ymax);
153  assert(nbinsx>0);
154  assert(nbinsy>0);
155 #ifndef NDEBUG
156  for (unsigned i = 0; i<nbinsy;++i)
157  assert(ybins[i]<ybins[i+1]);
158 #endif
159  for (unsigned i = 0; i<m_nbinsyPlus1;++i)
160  m_varBinsY[i]=ybins[i];
161 }
162 //____________________________________________________________________
163 template <class T>
164 template <class TFloat>
165 inline Flex2DHisto<T>::Flex2DHisto( unsigned nbinsx, const TFloat* xbins,
166  unsigned nbinsy, const TFloat* ybins )
167  : m_invDeltaX(nbinsx/(xbins[nbinsx]-xbins[0])),
168  m_invDeltaY(nbinsy/(ybins[nbinsy]-ybins[0])),
169  m_sumW(0),
170  m_sumW2(0),
171  m_sumWX(0),
172  m_sumWX2(0),
173  m_sumWY(0),
174  m_sumWY2(0),
175  m_sumWXY(0),
176  m_nEntries(0),
177  m_nbinsxPlus1(nbinsx+1),
178  m_nbinsyPlus1(nbinsy+1),
179  m_nbinsxPlus2(nbinsx+2),
180  m_xmin(xbins[0]),
181  m_xmax(xbins[nbinsx]),
182  m_ymin(ybins[0]),
183  m_ymax(ybins[nbinsy]),
184  m_varBinsX(LWPools::acquire<float>(nbinsx+1)),
185  m_varBinsY(LWPools::acquire<float>(nbinsy+1)),
186  m_flexArray((nbinsx+2)*(nbinsy+2))
187 {
188  assert(m_xmin<m_xmax);
189  assert(m_ymin<m_ymax);
190  assert(nbinsx>0);
191  assert(nbinsy>0);
192 #ifndef NDEBUG
193  for (unsigned i = 0; i<nbinsx;++i)
194  assert(xbins[i]<xbins[i+1]);
195  for (unsigned i = 0; i<nbinsy;++i)
196  assert(ybins[i]<ybins[i+1]);
197 #endif
198  for (unsigned i = 0; i<m_nbinsxPlus1;++i)
199  m_varBinsX[i]=xbins[i];
200  for (unsigned i = 0; i<m_nbinsyPlus1;++i)
201  m_varBinsY[i]=ybins[i];
202 }
203 
204 //____________________________________________________________________
205 template <class T>
206 inline Flex2DHisto<T>::~Flex2DHisto()
207 {
208  if (m_varBinsX)
209  LWPools::release(m_varBinsX,m_nbinsxPlus1);
210  if (m_varBinsY)
211  LWPools::release(m_varBinsY,m_nbinsyPlus1);
212 }
213 
214 //____________________________________________________________________
215 template <class T>
216 inline unsigned Flex2DHisto<T>::valueToXBin(const double& x) const
217 {
218  return LWBinUtils::valueToBin( x, m_varBinsX, m_invDeltaX, m_xmin, m_xmax, m_nbinsxPlus1 );
219 }
220 
221 //____________________________________________________________________
222 template <class T>
223 inline unsigned Flex2DHisto<T>::valueToYBin(const double& y) const
224 {
225  return LWBinUtils::valueToBin( y, m_varBinsY, m_invDeltaY,m_ymin, m_ymax, m_nbinsyPlus1 );
226 }
227 
228 
229 //____________________________________________________________________
230 template <class T>
231 inline double Flex2DHisto<T>::getBinCenterX(int bin) const
232 {
233  return LWBinUtils::getBinCenter( bin, m_varBinsX,m_invDeltaX, m_xmin,m_nbinsxPlus1);
234 }
235 
236 //____________________________________________________________________
237 template <class T>
238 inline double Flex2DHisto<T>::getBinCenterY(int bin) const
239 {
240  return LWBinUtils::getBinCenter( bin, m_varBinsY,m_invDeltaY, m_ymin,m_nbinsyPlus1);
241 }
242 
243 //____________________________________________________________________
244 template <class T>
245 inline void Flex2DHisto<T>::fill(const double& x, const double& y)
246 {
247  unsigned binx(valueToXBin(x)), biny(valueToYBin(y));
248 #ifndef NDEBUG
249  if (binx==USHRT_MAX||biny==USHRT_MAX)
250  return;
251 #endif
252  m_flexArray.fill(internal_bin(binx,biny));
253  //Update stats (sums not for over/under flow):
254  ++m_nEntries;
255  if (binx>0&&binx<m_nbinsxPlus1
256  &&biny>0&&biny<m_nbinsyPlus1) {
257  ++m_sumW;
258  ++m_sumW2;
259  m_sumWX += x;
260  m_sumWX2 += x*x;
261  m_sumWY += y;
262  m_sumWY2 += y*y;
263  m_sumWXY += x*y;
264  }
265 }
266 
267 //____________________________________________________________________
268 template <class T>
269 inline void Flex2DHisto<T>::fill(const double& x, const double& y, const double& w)
270 {
271 #ifndef NDEBUG
272  if (w!=w) {
273  std::cout<<"LWHisto: Saw NaN in fill weight"<<std::endl;
274  return;
275  }
276 #endif
277  unsigned binx(valueToXBin(x)), biny(valueToYBin(y));
278 #ifndef NDEBUG
279  if (binx==USHRT_MAX||biny==USHRT_MAX)
280  return;
281 #endif
282  m_flexArray.fill(internal_bin(binx,biny),w);
283  //Update stats (sums not for over/under flow):
284  ++m_nEntries;
285  if (binx>0&&binx<m_nbinsxPlus1
286  &&biny>0&&biny<m_nbinsyPlus1) {
287  //NB: root used to use fabs(w) instead of w in the formulas below:
288  m_sumW += w;
289  m_sumW2 += w*w;
290  const double tmpwy(w*y);
291  m_sumWY += tmpwy;
292  m_sumWY2 += tmpwy*y;
293  const double tmpwx(w*x);
294  m_sumWX += tmpwx;
295  m_sumWX2 += tmpwx*x;
296  m_sumWXY += tmpwx*y;
297  }
298 }
299 
300 //____________________________________________________________________
301 template <class T>
302 inline unsigned Flex2DHisto<T>::getEntries() const
303 {
304  return m_nEntries;
305 }
306 
307 //____________________________________________________________________
308 template <class T>
309 inline void Flex2DHisto<T>::setEntries(unsigned n)
310 {
311  m_nEntries = n;
312 }
313 
314 //____________________________________________________________________
315 template <class T>
316 inline double Flex2DHisto<T>::getSumW() const
317 {
318  return m_sumW;
319 }
320 
321 //____________________________________________________________________
322 template <class T>
323 inline double Flex2DHisto<T>::getSumW2() const
324 {
325  return m_sumW2;
326 }
327 
328 //____________________________________________________________________
329 template <class T>
330 inline double Flex2DHisto<T>::getSumWX() const
331 {
332  return m_sumWX;
333 }
334 
335 //____________________________________________________________________
336 template <class T>
337 inline double Flex2DHisto<T>::getSumWX2() const
338 {
339  return m_sumWX2;
340 }
341 
342 //____________________________________________________________________
343 template <class T>
344 inline double Flex2DHisto<T>::getSumWY() const
345 {
346  return m_sumWY;
347 }
348 
349 //____________________________________________________________________
350 template <class T>
351 inline double Flex2DHisto<T>::getSumWY2() const
352 {
353  return m_sumWY2;
354 }
355 
356 //____________________________________________________________________
357 template <class T>
358 inline double Flex2DHisto<T>::getSumWXY() const
359 {
360  return m_sumWXY;
361 }
362 
363 //____________________________________________________________________
364 template <class T>
365 inline void Flex2DHisto<T>::setSums( const double& sumW,
366  const double& sumW2,
367  const double& sumWX,
368  const double& sumWX2,
369  const double& sumWY,
370  const double& sumWY2,
371  const double& sumWXY )
372 {
373  m_sumW = sumW;
374  m_sumW2 = sumW2;
375  m_sumWX = sumWX;
376  m_sumWX2 = sumWX2;
377  m_sumWY = sumWY;
378  m_sumWY2 = sumWY2;
379  m_sumWXY = sumWXY;
380 }
381 
382 //____________________________________________________________________
383 template <class T>
384 inline void Flex2DHisto<T>::copyContents(T*__restrict__ cont, double*__restrict__ err) const
385 {
386  m_flexArray.copyContents(cont,err);
387 }
388 
389 //____________________________________________________________________
390 template <class T>
391 inline double Flex2DHisto<T>::Integral() const
392 {
393  double overflowbins( getBinContent(0,0)
394  +getBinContent(m_nbinsxPlus1,0)
395  +getBinContent(0,m_nbinsyPlus1)
396  +getBinContent(m_nbinsxPlus1,m_nbinsyPlus1));
397  for (unsigned ybin=1;ybin<m_nbinsyPlus1;++ybin) {
398  overflowbins += getBinContent(0,ybin);
399  overflowbins += getBinContent(m_nbinsxPlus1,ybin);
400  }
401  for (unsigned xbin=1;xbin<m_nbinsxPlus1;++xbin) {
402  overflowbins += getBinContent(xbin,0);
403  overflowbins += getBinContent(xbin,m_nbinsyPlus1);
404  }
405  return m_flexArray.Integral()-overflowbins;
406 }
407 
408 //____________________________________________________________________
409 template <class T>
410 inline bool Flex2DHisto<T>::holdsSeparateSumW2Info() const
411 {
412  return m_flexArray.holdsSeparateSumW2Info();
413 }
414 
415 //____________________________________________________________________
416 template <class T>
417 inline double Flex2DHisto<T>::getBinContent(unsigned binx, unsigned biny) const
418 {
419 #ifdef LW_STRICT_ROOT_BEHAVIOUR
420  if (binx>m_nbinsxPlus1) binx = m_nbinsxPlus1;
421  if (biny>m_nbinsyPlus1) biny = m_nbinsyPlus1;
422 #else
423  if (binx>m_nbinsxPlus1||biny>m_nbinsyPlus1)
424  return 0.0;
425 #endif
426  assert(binx<m_nbinsxPlus1+1);
427  assert(biny<m_nbinsyPlus1+1);
428  return m_flexArray.getBinContent(internal_bin(binx,biny));
429 }
430 
431 //____________________________________________________________________
432 template <class T>
433 inline double Flex2DHisto<T>::getBinError(unsigned binx,unsigned biny) const
434 {
435 #ifdef LW_STRICT_ROOT_BEHAVIOUR
436  if (binx>m_nbinsxPlus1) binx = m_nbinsxPlus1;
437  if (biny>m_nbinsyPlus1) biny = m_nbinsyPlus1;
438 #else
439  if (binx>m_nbinsxPlus1||biny>m_nbinsyPlus1)
440  return 0.0;
441 #endif
442  return m_flexArray.getBinError(internal_bin(binx,biny));
443 }
444 
445 //____________________________________________________________________
446 template <class T>
447 inline void Flex2DHisto<T>::getBinContentAndError(unsigned binx, unsigned biny, double& cont, double& err ) const
448 {
449  //float/integer version
450 #ifdef LW_STRICT_ROOT_BEHAVIOUR
451  if (binx>m_nbinsxPlus1) binx = m_nbinsxPlus1;
452  if (biny>m_nbinsyPlus1) biny = m_nbinsyPlus1;
453 #else
454  if (binx>m_nbinsxPlus1||biny>m_nbinsyPlus1) {
455  cont = 0;
456  err = 0;
457  return;
458  }
459 #endif
460  T tmp;
461  m_flexArray.getBinContentAndError(internal_bin(binx,biny),tmp,err);
462  cont = static_cast<double>(tmp);
463 }
464 
465 
466 //____________________________________________________________________
467 template <>
468 inline void Flex2DHisto<double>::getBinContentAndError(unsigned binx, unsigned biny, double& cont, double& err ) const
469 {
470  //double version
471 #ifdef LW_STRICT_ROOT_BEHAVIOUR
472  if (binx>m_nbinsxPlus1) binx = m_nbinsxPlus1;
473  if (biny>m_nbinsyPlus1) biny = m_nbinsyPlus1;
474 #else
475  if (binx>m_nbinsxPlus1||biny>m_nbinsyPlus1) {
476  cont = 0;
477  err = 0;
478  return;
479  }
480 #endif
481 
482  m_flexArray.getBinContentAndError(internal_bin(binx,biny),cont,err);
483 }
484 
485 //____________________________________________________________________
486 template <class T>
487 inline void Flex2DHisto<T>::setBinContent(unsigned binx, unsigned biny, const double& c)
488 {
489  //Weird root behaviour: For 2D version of setBinContent we wrap rather than ignore:
490 #ifdef LW_STRICT_ROOT_BEHAVIOUR
491  if (binx>m_nbinsxPlus1)
492  binx = m_nbinsxPlus1;
493  if (biny>m_nbinsyPlus1)
494  biny = m_nbinsyPlus1;
495 #else
496  if (binx>m_nbinsxPlus1||biny>m_nbinsyPlus1)
497  return;
498 #endif
499  m_flexArray.setBinContent(internal_bin(binx,biny),static_cast<T>(c));
500  //This is how it is done in ROOT, but is it really appropriate??:
501  ++m_nEntries;
502  m_sumW = 0;//As in ROOT
503 }
504 
505 //____________________________________________________________________
506 template <class T>
507 inline void Flex2DHisto<T>::setBinError(unsigned binx, unsigned biny, const double& e)
508 {
509  //Weird root behaviour: For 2D version of setBinError we ignore rather than wrap:
510  if (binx>m_nbinsxPlus1||biny>m_nbinsyPlus1) {//as in root
511  return;
512  }
513  m_flexArray.setBinError(internal_bin(binx,biny),e);
514 }
515 
516 //Weird ROOT behaviour (should file bug report): for get/set
517 //bincontent/binerror 1D/2D we have different choices of what to do
518 //with bins out of range... (constrain to valid value vs. ignore).
519 
520 //____________________________________________________________________
521 template <class T>
522 inline void Flex2DHisto<T>::setBinContentAndError(unsigned binx, unsigned biny, const double& cont, const double& err )
523 {
524  if (binx>m_nbinsxPlus1||biny>m_nbinsyPlus1) {
525 #ifdef LW_STRICT_ROOT_BEHAVIOUR
526  //Annoyingly enough, root has different behaviour for SetBinContent
527  //(wrap) and SetBinError (ignore) in this case:
528  setBinContent(binx,biny,cont);
529 #endif
530  return;
531  }
532  m_flexArray.setBinContentAndError(internal_bin(binx,biny),static_cast<T>(cont),err);
533  //FIXME: It seems that the next two lines are not executed for the last bin in ROOT!!!!!!!
534  ++m_nEntries;
535  m_sumW = 0;//As in ROOT
536 }
537 
538 //____________________________________________________________________
539 template <class T>
540 inline void Flex2DHisto<T>::resetActiveBinLoop()
541 {
542  m_flexArray.resetActiveBinLoop();
543 }
544 
545 //____________________________________________________________________
546 template <class T>
547 inline bool Flex2DHisto<T>::getNextActiveBin(unsigned& binx, unsigned&biny, double& content, double& error)
548 {
549  //float/integer version
550  unsigned internalbin;
551  T tmp;
552  if (m_flexArray.getNextActiveBin(internalbin,tmp,error)) {
553  binx = internalbin % m_nbinsxPlus2;
554  biny = (internalbin-binx) / m_nbinsxPlus2;
555  assert(internal_bin(binx,biny)==internalbin);
556  content = static_cast<double>(tmp);
557  return true;
558  } else {
559  return false;
560  }
561 }
562 
563 //____________________________________________________________________
564 template <>
565 inline bool Flex2DHisto<double>::getNextActiveBin(unsigned& binx, unsigned&biny, double& content, double& error)
566 {
567  //double version
568  unsigned internalbin;
569  if (m_flexArray.getNextActiveBin(internalbin,content,error)) {
570  LWBinUtils::unpack_internal_bin(internalbin,binx,biny,m_nbinsxPlus2);
571  return true;
572  } else {
573  return false;
574  }
575 }
576 
577 //____________________________________________________________________
578 template <class T>
579 inline void Flex2DHisto<T>::scaleContentsAndErrors( const double& fact )
580 {
581  m_sumW *= fact;
582  m_sumW2 *= fact*fact;
583  m_sumWX *= fact;
584  m_sumWX2 *= fact;
585  m_sumWY *= fact;
586  m_sumWY2 *= fact;
587  m_sumWXY *= fact;
588  m_flexArray.scaleContentsAndErrors(fact);
589 }
590 
591 
592 #ifdef LW_DEBUG_HEAVY_USERS
593 
594 //____________________________________________________________________
595 template <class T>
596 inline void Flex2DHisto<T>::countCall(const std::pair<void*,void*>&addresses)
597 {
598  std::map<std::pair<void*,void*>,unsigned long>::iterator it = m_callmap.find(addresses);
599  if (it==m_callmap.end())
600  m_callmap[addresses]=0;
601  else
602  ++(it->second);
603 }
604 
605 //____________________________________________________________________
606 template <class T>
607 inline void Flex2DHisto<T>::produceReport(const char*histname)
608 {
609  std::map<std::pair<void*,void*>,unsigned long>::iterator it,itE(m_callmap.end());
610  for(it=m_callmap.begin();it!=itE;++it) {
611  if (it->second>0) {
612  char * caller = LWHistTraceUtils::getSymbol(it->first.second);
613  char * calledmethod = LWHistTraceUtils::getSymbol(it->first.first);
614  std::cout<<"LWHists WARNING: Method in histogram called "<<it->second<<" times: "<<calledmethod<<" from "<<caller<<" (histogram named \""<<histname<<"\")"<<std::endl;
615  //fixme: do free on caller and calledmethod
616  }
617  }
618 }
619 
620 #endif
Generated on Thu Jul 4 2024 21:15:52 for ATLAS Offline Software by   doxygen 1.8.18