ATLAS Offline Software
Public Member Functions | Static Public Member Functions | Private Attributes | List of all members
Analysis::HistoHelperRoot Class Reference

Helper class for histograming. More...

#include <HistoHelperRoot.h>

Collaboration diagram for Analysis::HistoHelperRoot:

Public Member Functions

 HistoHelperRoot (ITHistSvc *)
 Constructor with histoSvc only. More...
 
 ~HistoHelperRoot ()
 
std::string baseName (const std::string &fullHistoName)
 
void bookHisto (const std::string &histoName, const std::string &histoTitle, unsigned int bins, double minx, double maxx)
 
void bookHisto (const std::string &histoName, const std::string &histoTitle, unsigned int bins, double *edge)
 
void bookHisto (const std::string &histoName, const std::string &histoTitle, unsigned int binsx, double minx, double maxx, unsigned int binsy, double miny, double maxy)
 
void bookHisto (const std::string &histoName, const std::string &histoTitle, unsigned int binsx, double *edgex, unsigned int binsy, double *edgey)
 
void bookHisto (const std::string &histoName, const std::string &histoTitle, unsigned int binsx, double minx, double maxx, unsigned int binsy, double miny, double maxy, unsigned int binsz, double minz, double maxz)
 
void fillHisto (const std::string &histoName, double) const
 
void fillHistoWithWeight (const std::string &histoName, double, double) const
 
void fillHisto (const std::string &histoName, double, double) const
 
void fillHisto (const std::string &histoName, double, double, double) const
 
void setCheckOverflows (bool b)
 
void print ()
 
TH1 *getHisto1D ATLAS_NOT_THREAD_SAFE (const std::string &histoName)
 Histogram accessors (not thread-safe because exposed bare pointer) More...
 
TH2 *getHisto2D ATLAS_NOT_THREAD_SAFE (const std::string &histoName)
 
TH3 *getHisto3D ATLAS_NOT_THREAD_SAFE (const std::string &histoName)
 

Static Public Member Functions

static double Interpol1d (double, TH1 *)
 
static double Interpol2d (double, double, TH2 *)
 
static double Interpol3d (double, double, double, TH3 *)
 
static void smoothASH2D (TH2 *, int m1=3, int m2=3, bool debug=false)
 
static void smoothASH3D (TH3 *, int m1=3, int m2=3, int m3=2, bool debug=false)
 

Private Attributes

ITHistSvc * m_histoSvc
 
std::map< std::string, LockedHandle< TH1 > > m_histoMap1D
 
std::map< std::string, LockedHandle< TH2 > > m_histoMap2D
 
std::map< std::string, LockedHandle< TH3 > > m_histoMap3D
 
std::map< std::string, HistoLimits > m_histoLimitsMap1D
 
std::map< std::string, HistoLimits > m_histoLimitsMap2D
 
std::map< std::string, HistoLimits > m_histoLimitsMap3D
 
bool m_checkOverflows
 

Detailed Description

Helper class for histograming.

The class is made thread-safe through the use of LockedHandle for histogram filling.

Definition at line 35 of file HistoHelperRoot.h.

Constructor & Destructor Documentation

◆ HistoHelperRoot()

Analysis::HistoHelperRoot::HistoHelperRoot ( ITHistSvc *  histoSvc)

Constructor with histoSvc only.

Definition at line 28 of file HistoHelperRoot.cxx.

28  :
29  m_histoSvc(histoSvc),
30  m_checkOverflows(true)
31  {}

◆ ~HistoHelperRoot()

Analysis::HistoHelperRoot::~HistoHelperRoot ( )

Definition at line 34 of file HistoHelperRoot.cxx.

34 {}

Member Function Documentation

◆ ATLAS_NOT_THREAD_SAFE() [1/3]

TH1* getHisto1D Analysis::HistoHelperRoot::ATLAS_NOT_THREAD_SAFE ( const std::string &  histoName)

Histogram accessors (not thread-safe because exposed bare pointer)

◆ ATLAS_NOT_THREAD_SAFE() [2/3]

TH2* getHisto2D Analysis::HistoHelperRoot::ATLAS_NOT_THREAD_SAFE ( const std::string &  histoName)

◆ ATLAS_NOT_THREAD_SAFE() [3/3]

TH3* getHisto3D Analysis::HistoHelperRoot::ATLAS_NOT_THREAD_SAFE ( const std::string &  histoName)

◆ baseName()

std::string Analysis::HistoHelperRoot::baseName ( const std::string &  fullHistoName)

Definition at line 36 of file HistoHelperRoot.cxx.

36  {
37  const std::string delim("/");
38  std::vector<std::string> words;
39  std::string::size_type sPos, sEnd, sLen;
40  sPos = fullHistoName.find_first_not_of(delim);
41  while ( sPos != std::string::npos ) {
42  sEnd = fullHistoName.find_first_of(delim, sPos);
43  if(sEnd==std::string::npos) sEnd = fullHistoName.length();
44  sLen = sEnd - sPos;
45  words.emplace_back(fullHistoName.substr(sPos,sLen));
46  sPos = fullHistoName.find_first_not_of(delim, sEnd);
47  }
48  std::string base = "";
49  if(words.size()>0) base = words[words.size()-1];
50  return base;
51 }

◆ bookHisto() [1/5]

void Analysis::HistoHelperRoot::bookHisto ( const std::string &  histoName,
const std::string &  histoTitle,
unsigned int  bins,
double *  edge 
)

Definition at line 63 of file HistoHelperRoot.cxx.

64 {
65  auto h = std::make_unique<TH1F>(this->baseName(histoName).c_str(),histoTitle.c_str(),bins,edge);
66  if (!m_histoSvc->regShared(histoName, std::move(h), m_histoMap1D[histoName]).isSuccess()) {
67  std::cout <<"Cannot book histo " << histoName << " with Title " << histoTitle <<std::endl;
68  }
69  m_histoLimitsMap1D[histoName] = HistoLimits(bins, edge[0], edge[bins]);
70 }

◆ bookHisto() [2/5]

void Analysis::HistoHelperRoot::bookHisto ( const std::string &  histoName,
const std::string &  histoTitle,
unsigned int  bins,
double  minx,
double  maxx 
)

Definition at line 53 of file HistoHelperRoot.cxx.

54 {
55  auto h = std::make_unique<TH1F>(this->baseName(histoName).c_str(),histoTitle.c_str(),bins,minx,maxx);
56  //std::cout<<"HistoHelperRoot: booking 1D "<<histoName<<std::endl;
57  if (!m_histoSvc->regShared(histoName, std::move(h), m_histoMap1D[histoName]).isSuccess()) {
58  std::cout <<"Cannot book histo " << histoName << " with Title " << histoTitle <<std::endl;
59  }
60  m_histoLimitsMap1D[histoName] = HistoLimits(bins, minx, maxx);
61 }

◆ bookHisto() [3/5]

void Analysis::HistoHelperRoot::bookHisto ( const std::string &  histoName,
const std::string &  histoTitle,
unsigned int  binsx,
double *  edgex,
unsigned int  binsy,
double *  edgey 
)

Definition at line 81 of file HistoHelperRoot.cxx.

82 {
83  auto h = std::make_unique<TH2F>(this->baseName(histoName).c_str(),histoTitle.c_str(),binsx,edgex,binsy,edgey);
84  if (!m_histoSvc->regShared(histoName, std::move(h), m_histoMap2D[histoName]).isSuccess()) {
85  std::cout <<"Cannot book histo " << histoName << " with Title " << histoTitle <<std::endl;
86  }
87  m_histoLimitsMap2D[histoName] = HistoLimits(binsx, edgex[0], edgex[binsx], binsy, edgey[0], edgey[binsy]);
88 }

◆ bookHisto() [4/5]

void Analysis::HistoHelperRoot::bookHisto ( const std::string &  histoName,
const std::string &  histoTitle,
unsigned int  binsx,
double  minx,
double  maxx,
unsigned int  binsy,
double  miny,
double  maxy 
)

Definition at line 72 of file HistoHelperRoot.cxx.

73 {
74  auto h = std::make_unique<TH2F>(this->baseName(histoName).c_str(),histoTitle.c_str(),binsx,minx,maxx,binsy,miny,maxy);
75  if (!m_histoSvc->regShared(histoName, std::move(h), m_histoMap2D[histoName]).isSuccess()) {
76  std::cout <<"Cannot book histo " << histoName << " with Title " << histoTitle <<std::endl;
77  }
78  m_histoLimitsMap2D[histoName] = HistoLimits(binsx, minx, maxx, binsy, miny, maxy);
79 }

◆ bookHisto() [5/5]

void Analysis::HistoHelperRoot::bookHisto ( const std::string &  histoName,
const std::string &  histoTitle,
unsigned int  binsx,
double  minx,
double  maxx,
unsigned int  binsy,
double  miny,
double  maxy,
unsigned int  binsz,
double  minz,
double  maxz 
)

Definition at line 90 of file HistoHelperRoot.cxx.

91 {
92  auto h = std::make_unique<TH3F>(this->baseName(histoName).c_str(),histoTitle.c_str(),binsx,minx,maxx,binsy,miny,maxy,binsz,minz,maxz);
93  if (!m_histoSvc->regShared(histoName, std::move(h), m_histoMap3D[histoName]).isSuccess()) {
94  std::cout <<"Cannot book histo " << histoName << " with Title " << histoTitle <<std::endl;
95  }
96  m_histoLimitsMap3D[histoName] = HistoLimits(binsx, minx, maxx, binsy, miny, maxy, binsz, minz, maxz);
97 }

◆ fillHisto() [1/3]

void Analysis::HistoHelperRoot::fillHisto ( const std::string &  histoName,
double  value 
) const

make sure that underflow/overflow is filled in 1st/last bin and not in 1st-1 and last+1 bin

Definition at line 101 of file HistoHelperRoot.cxx.

102 {
105  if (m_checkOverflows) {
106  const HistoLimits& l = (*(m_histoLimitsMap1D.find(histoName))).second;
107  if ( value < l.xmin ) value = l.xmin+0.0001;
108  else if ( value > l.xmax ) value = l.xmax-0.0001;
109  }
110  // Thread-safe because handle uses internal lock
111  auto locked_handle ATLAS_THREAD_SAFE = const_cast<LockedHandle<TH1>&>(m_histoMap1D.at(histoName));
112  locked_handle->Fill(value);
113 }

◆ fillHisto() [2/3]

void Analysis::HistoHelperRoot::fillHisto ( const std::string &  histoName,
double  valuex,
double  valuey 
) const

make sure that underflow/overflow is filled in 1st/last bin and not in 1st-1 and last+1 bin

Definition at line 129 of file HistoHelperRoot.cxx.

130 {
133  if (m_checkOverflows) {
134  const HistoLimits& l = (*(m_histoLimitsMap2D.find(histoName))).second;
135  if ( valuex < l.xmin ) valuex = l.xmin+0.0001;
136  else if ( valuex > l.xmax ) valuex = l.xmax-0.0001;
137  if ( valuey < l.ymin ) valuey = l.ymin+0.0001;
138  else if ( valuey > l.ymax ) valuey = l.ymax-0.0001;
139  }
140  // Thread-safe because handle uses internal lock
141  auto locked_handle ATLAS_THREAD_SAFE = const_cast<LockedHandle<TH2>&>(m_histoMap2D.at(histoName));
142  locked_handle->Fill(valuex,valuey);
143 }

◆ fillHisto() [3/3]

void Analysis::HistoHelperRoot::fillHisto ( const std::string &  histoName,
double  valuex,
double  valuey,
double  valuez 
) const

make sure that underflow/overflow is filled in 1st/last bin and not in 1st-1 and last+1 bin

Definition at line 145 of file HistoHelperRoot.cxx.

146 {
149  if (m_checkOverflows) {
150  const HistoLimits& l = (*(m_histoLimitsMap3D.find(histoName))).second;
151  if ( valuex < l.xmin ) valuex = l.xmin+0.0001;
152  else if ( valuex > l.xmax ) valuex = l.xmax-0.0001;
153  if ( valuey < l.ymin ) valuey = l.ymin+0.0001;
154  else if ( valuey > l.ymax ) valuey = l.ymax-0.0001;
155  if ( valuez < l.zmin ) valuez = l.zmin+0.0001;
156  else if ( valuez > l.zmax ) valuez = l.zmax-0.0001;
157  }
158  // Thread-safe because handle uses internal lock
159  auto locked_handle ATLAS_THREAD_SAFE = const_cast<LockedHandle<TH3>&>(m_histoMap3D.at(histoName));
160  locked_handle->Fill(valuex,valuey,valuez);
161 }

◆ fillHistoWithWeight()

void Analysis::HistoHelperRoot::fillHistoWithWeight ( const std::string &  histoName,
double  value,
double  weight 
) const

make sure that underflow/overflow is filled in 1st/last bin and not in 1st-1 and last+1 bin

Definition at line 115 of file HistoHelperRoot.cxx.

116 {
119  if (m_checkOverflows) {
120  const HistoLimits& l = (*(m_histoLimitsMap1D.find(histoName))).second;
121  if ( value < l.xmin ) value = l.xmin+0.0001;
122  else if ( value > l.xmax ) value = l.xmax-0.0001;
123  }
124  // Thread-safe because handle uses internal lock
125  auto locked_handle ATLAS_THREAD_SAFE = const_cast<LockedHandle<TH1>&>(m_histoMap1D.at(histoName));
126  locked_handle->Fill(value, weight);
127 }

◆ Interpol1d()

double Analysis::HistoHelperRoot::Interpol1d ( double  x,
TH1 *  h 
)
static

Definition at line 578 of file HistoHelperRoot.cxx.

579 {
580  int bin((h->GetXaxis())->FindBin(x));
581  double bincenter(h->GetBinCenter(bin));
582  double nextbincenter(bincenter);
583  double nextbincontent(0);
584  if(x>bincenter && bin < h->GetNbinsX())
585  {
586  nextbincenter = h->GetBinCenter(bin+1);
587  nextbincontent = h->GetBinContent(bin+1);
588  }
589  else if(x<bincenter && bin > 1)
590  {
591  nextbincenter = h->GetBinCenter(bin-1);
592  nextbincontent = h->GetBinContent(bin-1);
593  }
594  double tmp(h->GetBinContent(bin));
595  if(nextbincenter != bincenter) tmp = ( nextbincontent*(x-bincenter) + tmp*(nextbincenter-x) ) / (nextbincenter - bincenter);
596  return tmp;
597 }

◆ Interpol2d()

double Analysis::HistoHelperRoot::Interpol2d ( double  x,
double  y,
TH2 *  h 
)
static

Definition at line 602 of file HistoHelperRoot.cxx.

603 {
604  int nx = h->GetNbinsX(), ny = h->GetNbinsY();
605  int ibx = h->GetXaxis()->FindBin(x), iby = h->GetYaxis()->FindBin(y);
606  int ibx2,iby2;
607  double z00,z11,z01,z10,xc,yc,xc2,yc2,u,t,r;
608  if (ibx > nx) ibx = nx;
609  if (iby > ny) iby = ny;
610  xc = h->GetXaxis()->GetBinCenter(ibx);
611  yc = h->GetYaxis()->GetBinCenter(iby);
612  z11 = h->GetBinContent(ibx,iby);
613  if ( (ibx > 1 || (ibx == 1 && x > xc) ) &&
614  (ibx < nx || (ibx == nx && x < xc) ) ) {
615  if (x > xc) {ibx2 = ibx + 1;} else {ibx2 = ibx - 1;}
616  xc2 = h->GetXaxis()->GetBinCenter(ibx2);
617  if ( (iby > 1 || (iby == 1 && y > yc) ) &&
618  (iby < ny || (iby == ny && y < yc) ) ) {
619  if (y > yc) {iby2 = iby + 1;} else {iby2 = iby - 1;}
620  yc2 = h->GetYaxis()->GetBinCenter(iby2);
621  z00 = h->GetBinContent(ibx2,iby2);
622  z10 = h->GetBinContent(ibx,iby2);
623  z01 = h->GetBinContent(ibx2,iby);
624  t = (x - xc2)/(xc-xc2);
625  u = (y - yc2)/(yc-yc2);
626  r = z11*t*u+z00*(1.-t)*(1.-u)+z01*(1.-t)*u+z10*t*(1.-u);
627  } else {
628  z01 = h->GetBinContent(ibx2,iby);
629  t = (x - xc2)/(xc-xc2);
630  r = z11*t + z01*(1.-t);
631  }
632  } else if ((iby > 1 || (iby == 1 && y > yc) ) &&
633  (iby < ny || (iby == ny && y < yc) ) ) {
634  if (y > yc) {iby2 = iby + 1;} else {iby2 = iby - 1;}
635  z10 = h->GetBinContent(ibx,iby2);
636  yc2 = h->GetYaxis()->GetBinCenter(iby2);
637  u = (y - yc2)/(yc-yc2);
638  r = z11*u + z10*(1.-u);
639  } else {
640  r = z11;
641  }
642  return r;
643 }

◆ Interpol3d()

double Analysis::HistoHelperRoot::Interpol3d ( double  x,
double  y,
double  z,
TH3 *  h 
)
static

Definition at line 648 of file HistoHelperRoot.cxx.

649 {
650  int nx = h->GetNbinsX(), ny = h->GetNbinsY(), nz = h->GetNbinsZ();
651  int ibx = h->GetXaxis()->FindBin(x), iby = h->GetYaxis()->FindBin(y), ibz = h->GetZaxis()->FindBin(z);
652  int ibx2,iby2,ibz2;
653  double z000,z010,z110,z100,z001,z011,z111,z101,xc,yc,zc,xc2,yc2,zc2,u,t,v,r;
654  if (ibx > nx) ibx = nx;
655  if (iby > ny) iby = ny;
656  if (ibz > nz) ibz = nz;
657  xc = h->GetXaxis()->GetBinCenter(ibx);
658  yc = h->GetYaxis()->GetBinCenter(iby);
659  zc = h->GetZaxis()->GetBinCenter(ibz);
660  //
661  z111 = h->GetBinContent(ibx,iby,ibz);
662  if ( (ibx > 1 || (ibx == 1 && x > xc) ) &&
663  (ibx < nx || (ibx == nx && x < xc) ) ) {
664  if (x > xc) {ibx2 = ibx + 1;} else {ibx2 = ibx - 1;}
665  xc2 = h->GetXaxis()->GetBinCenter(ibx2);
666  if ( (iby > 1 || (iby == 1 && y > yc) ) &&
667  (iby < ny || (iby == ny && y < yc) ) ) {
668  if (y > yc) {iby2 = iby + 1;} else {iby2 = iby - 1;}
669  yc2 = h->GetYaxis()->GetBinCenter(iby2);
670  if ( (ibz > 1 || (ibz == 1 && z > zc) ) &&
671  (ibz < nz || (ibz == nz && z < zc) ) ) {
672  if (z > zc) {ibz2 = ibz + 1;} else {ibz2 = ibz - 1;}
673  zc2 = h->GetZaxis()->GetBinCenter(ibz2);
674  //
675  z000 = h->GetBinContent(ibx2,iby2,ibz2);
676  z100 = h->GetBinContent(ibx,iby2,ibz2);
677  z010 = h->GetBinContent(ibx2,iby,ibz2);
678  z110 = h->GetBinContent(ibx,iby,ibz2);
679  z001 = h->GetBinContent(ibx2,iby2,ibz);
680  z101 = h->GetBinContent(ibx,iby2,ibz);
681  z011 = h->GetBinContent(ibx2,iby,ibz);
682  t = (x - xc2)/(xc-xc2);
683  u = (y - yc2)/(yc-yc2);
684  v = (z - zc2)/(zc-zc2);
685  r = z111*t*u*v + z001*(1.-t)*(1.-u)*v
686  + z011*(1.-t)*u*v + z101*t*(1.-u)*v
687  + z110*t*u*(1.-v) + z000*(1.-t)*(1.-u)*(1.-v)
688  + z010*(1.-t)*u*(1.-v) + z100*t*(1.-u)*(1.-v);
689  } else {
690  z011 = h->GetBinContent(ibx2,iby,ibz);
691  z001 = h->GetBinContent(ibx2,iby2,ibz);
692  z101 = h->GetBinContent(ibx,iby2,ibz);
693  t = (x - xc2)/(xc-xc2);
694  u = (y - yc2)/(yc-yc2);
695  r = z111*t*u + z011*(1.-t)*u + z101*t*(1.-u) + z001*(1.-t)*(1.-u) ;
696  }
697  } else if ((ibz > 1 || (ibz == 1 && z > zc) ) &&
698  (ibz < nz || (ibz == nz && z < zc) ) ) {
699  if (z > zc) {ibz2 = ibz + 1;} else {ibz2 = ibz - 1;}
700  z110 = h->GetBinContent(ibx,iby,ibz2);
701  z010 = h->GetBinContent(ibx2,iby,ibz2);
702  z011 = h->GetBinContent(ibx2,iby,ibz);
703  zc2 = h->GetYaxis()->GetBinCenter(ibz2);
704  t = (x - xc2)/(xc-xc2);
705  v = (z - zc2)/(zc-zc2);
706  r = z111*t*v + z011*(1.-t)*v + z110*t*(1.-v) +z010*(1.-t)*(1.-v);
707  } else {
708  z011 = h->GetBinContent(ibx2,iby,ibz);
709  t = (x - xc2)/(xc-xc2);
710  r = z111*t + z011*(1.-t);
711  }
712  } else {
713  if ( (iby > 1 || (iby == 1 && y > yc) ) &&
714  (iby < ny || (iby == ny && y < yc) ) ) {
715  if (y > yc) {iby2 = iby + 1;} else {iby2 = iby - 1;}
716  yc2 = h->GetYaxis()->GetBinCenter(iby2);
717  if ( (ibz > 1 || (ibz == 1 && z > zc) ) &&
718  (ibz < nz || (ibz == nz && z < zc) ) ) {
719  if (z > zc) {ibz2 = ibz + 1;} else {ibz2 = ibz - 1;}
720  zc2 = h->GetZaxis()->GetBinCenter(ibz2);
721  //
722  z100 = h->GetBinContent(ibx,iby2,ibz2);
723  z110 = h->GetBinContent(ibx,iby,ibz2);
724  z101 = h->GetBinContent(ibx,iby2,ibz);
725  u = (y - yc2)/(yc-yc2);
726  v = (z - zc2)/(zc-zc2);
727  r = z111*u*v + z101*(1.-u)*v
728  + z110*u*(1.-v) + z100*(1.-u)*(1.-v);
729  } else {
730  z101 = h->GetBinContent(ibx,iby2,ibz);
731  u = (y - yc2)/(yc-yc2);
732  r = z111*u + z101*(1.-u);
733  }
734  } else if ((ibz > 1 || (ibz == 1 && z > zc) ) &&
735  (ibz < nz || (ibz == nz && z < zc) ) ) {
736  if (z > zc) {ibz2 = ibz + 1;} else {ibz2 = ibz - 1;}
737  z110 = h->GetBinContent(ibx,iby,ibz2);
738  zc2 = h->GetYaxis()->GetBinCenter(ibz2);
739  v = (z - zc2)/(zc-zc2);
740  r = z111*v + z110*(1.-v);
741  } else {
742  r = z111;
743  }
744  }
745  return r;
746 }

◆ print()

void Analysis::HistoHelperRoot::print ( )

Definition at line 176 of file HistoHelperRoot.cxx.

176  {
177  for (const auto& [name, h] : m_histoMap1D) {
178  std::cout <<"The name of the 1D Histo " << name << std::endl;
179  }
180  for (const auto& [name, h] : m_histoMap2D) {
181  std::cout <<"The name of the 2D Histo " << name << std::endl;
182  }
183  for (const auto& [name, h] : m_histoMap3D) {
184  std::cout <<"The name of the 3D Histo " << name << std::endl;
185  }
186 }

◆ setCheckOverflows()

void Analysis::HistoHelperRoot::setCheckOverflows ( bool  b)
inline

Definition at line 71 of file HistoHelperRoot.h.

71 {m_checkOverflows = b;};

◆ smoothASH2D()

void Analysis::HistoHelperRoot::smoothASH2D ( TH2 *  input2D,
int  m1 = 3,
int  m2 = 3,
bool  debug = false 
)
static

Definition at line 188 of file HistoHelperRoot.cxx.

188  {
189  if(debug) {
190  std::cout << "Smoothing a two dimensional histogram "<< input2D->GetName()
191  << " " << m1 << " " << m2 << std::endl;
192  std::cout << "First (1-3, 1-3) 3x3 bins before smoothing: " << std::endl;
193  for(int i=1;i<4;i++) {
194  for(int j=1;j<4;j++) {
195  std::cout<<i<<" "<<j<<" : "<<input2D->GetBinContent(i,j)<< " / " ;
196  }
197  }
198  std::cout << std::endl;
199  int ioffset = input2D->GetNbinsX() / 2 ;
200  int joffset = input2D->GetNbinsY() / 2 ;
201  std::cout << "Middle (" << ioffset+1 << "-" << ioffset+4 << ", ("
202  << joffset+1 << "-" << joffset+4 << ") 3x3 bins before smoothing: " << std::endl;
203  for(int i=1;i<4;i++) {
204  for(int j=1;j<4;j++) {
205  std::cout<<i<<" "<<j<<" : "<<input2D->GetBinContent(i+ioffset,j+joffset)<< " / " ;
206  }
207  }
208  std::cout << std::endl;
209  }
210  //
211  const int lsup = 20;
212  if (m1 > lsup || m2 > lsup) {
213  std::cout <<"HistoHelperRoot::smoothASH2D: m1 or m2 too big !"<<std::endl;
214  return;
215  } else {
216  int nx = input2D->GetNbinsX()+1;
217  int ny = input2D->GetNbinsY()+1;
218  float **h, **res;
219  h = new float*[nx-1];
220  res = new float*[nx-1];
221  for (int i = 0;i < nx-1;i++) {
222  h[i] = new float[ny-1];
223  res[i] = new float[ny-1];
224  }
225  for (int iy = 1;iy<ny;iy++) {
226  for (int ix = 1;ix<nx;ix++) {
227  h[ix-1][iy-1] = (float) input2D->GetBinContent(ix,iy);
228  }
229  }
230  //
231  int i,j,k,l;
232  float wk1[41],wk2[41];
233  //avoid large stack use
234  std::vector<std::vector<float>> wgt(100, std::vector<float>(100));
235  std::vector<std::vector<float>> wk(41, std::vector<float>(41));
236  double wks = 0.;
237  float ai,am1 = float(m1), am2 = float(m2);
238  const float am12 = am1*am1, am22 = am2*am2;
239  const float inv_am1_am2 = 1. / (am1 * am2);
240  const float inv_am12 = 1. / am12;
241  const float inv_am22 = 1. / am22;
242  // Initialisation
243  for (k = 0;k<nx-1;k++) {
244  for (l = 0;l<ny-1;l++) {
245  res[k][l] = 0.; wgt[k][l] = 0.;
246  }
247  }
248  // Weights
249  for (i = lsup+1-m1;i<lsup+m1;i++) {
250  ai = float(i-lsup)*float(i-lsup);
251  wk1[i] = 15./16.*(1.-ai*inv_am12)*(1.-ai*inv_am12);
252  wks = wks + wk1[i];
253  }
254  if (wks == 0)[[unlikely]]{
255  std::cout <<"HistoHelperRoot::smoothASH2D: wks is zero! "<<std::endl;
256  for (int i = 0; i < nx-1; ++i) {
257  delete[] h[i];
258  delete[] res[i];
259  }
260 
261  delete[] h;
262  delete[] res;
263  return;
264  }
265  const double fac1 = am1 / wks;
266  for (i = lsup+1-m1;i<lsup+m1;i++) {
267  wk1[i] = wk1[i]*fac1;
268  }
269  wks = 0.;
270  for (i = lsup+1-m2;i<lsup+m2;i++) {
271  ai = float(i-lsup)*float(i-lsup);
272  wk2[i] = 15./16.*(1.-ai*inv_am22)*(1.-ai*inv_am22);
273  wks = wks + wk2[i];
274  }
275  if (wks == 0)[[unlikely]]{
276  std::cout <<"HistoHelperRoot::smoothASH2D: wks is zero! "<<std::endl;
277  for (int i = 0; i < nx-1; ++i) {
278  delete[] h[i];
279  delete[] res[i];
280  }
281  delete[] h;
282  delete[] res;
283 
284  return;
285  }
286  const double fac2 = am2 / wks;
287  for (i = lsup+1-m2;i<lsup+m2;i++) {
288  wk2[i] = wk2[i]*fac2;
289  }
290  for (i = lsup+1-m1;i<lsup+m1;i++) {
291  for (j = lsup+1-m2;j<lsup+m2;j++) {
292  wk[i][j] = wk1[i]*wk2[j];
293  }
294  }
295  //
296  for (k = 0;k<nx-1;k++) {
297  for (l = 0;l<ny-1;l++) {
298  for (i = std::max(0,k-m1+1);i<std::min(nx-1,k+m1);i++) {
299  for (j = std::max(0,l-m2+1);j<std::min(ny-1,l+m2);j++) {
300  res[i][j] = res[i][j] + wk[lsup+i-k][lsup+j-l]*h[k][l];
301  wgt[i][j] = wgt[i][j] + wk[lsup+i-k][lsup+j-l];
302  }
303  }
304  }
305  }
306  for (k = 0;k<nx-1;k++) {
307  for (l = 0;l<ny-1;l++) {
308  res[k][l] = res[k][l]*inv_am1_am2;
309  if (wgt[k][l] != 0.) {res[k][l] = res[k][l]/wgt[k][l];}
310  }
311  }
312  // replace the histo content with the result of smoothing:
313  input2D->Reset();
314  for (int iy = 1;iy<ny;iy++) {
315  for (int ix = 1;ix<nx;ix++) {
316  input2D->SetBinContent(ix,iy,res[ix-1][iy-1]);
317  }
318  }
319  for (i = 0; i < nx-1; i++){
320  delete[] h[i];
321  delete[] res[i];
322  }
323  delete[] h;
324  delete[] res;
325  }
326  if(debug) {
327  std::cout << "First (1-3, 1-3) 3x3 bins after smoothing: " << std::endl;
328  for(int i=1;i<4;i++) {
329  for(int j=1;j<4;j++) {
330  std::cout<<i<<" "<<j<<" : "<<input2D->GetBinContent(i,j)<< " / " ;
331  }
332  }
333  std::cout << std::endl;
334  int ioffset = input2D->GetNbinsX() / 2 ;
335  int joffset = input2D->GetNbinsY() / 2 ;
336  std::cout << "Middle (" << ioffset+1 << "-" << ioffset+4 << ", ("
337  << joffset+1 << "-" << joffset+4 << ") 3x3 bins after smoothing: " << std::endl;
338  for(int i=1;i<4;i++) {
339  for(int j=1;j<4;j++) {
340  std::cout<<i<<" "<<j<<" : "<<input2D->GetBinContent(i+ioffset,j+joffset)<< " / " ;
341  }
342  }
343  std::cout << std::endl;
344  }
345 }

◆ smoothASH3D()

void Analysis::HistoHelperRoot::smoothASH3D ( TH3 *  input3D,
int  m1 = 3,
int  m2 = 3,
int  m3 = 2,
bool  debug = false 
)
static

Definition at line 347 of file HistoHelperRoot.cxx.

347  {
348  if(debug) {
349  std::cout << "Smoothing a three dimensional histogram "<< input3D->GetName()
350  << " " << m1 << " " << m2 << " " << m3 << std::endl;
351  std::cout << "First 2x2x2 bins before smoothing: " << std::endl;
352  for(int i=1;i<3;i++) {
353  for(int j=1;j<3;j++) {
354  for(int k=1;k<3;k++) {
355  std::cout<<i<<" "<<j<<" "<<k<<" : "<<input3D->GetBinContent(i,j,k)<< " / " ;
356  }
357  }
358  }
359  std::cout << std::endl;
360  int ioffset = input3D->GetNbinsX() / 2 ;
361  int joffset = input3D->GetNbinsY() / 2 ;
362  int koffset = input3D->GetNbinsZ() / 2 ;
363  std::cout << "Middle (" << ioffset+1 << "-" << ioffset+3 << ", "
364  << joffset+1 << "-" << joffset+3 << ", "
365  << koffset+1 << "-" << koffset+3 << ") 2x2 bins before smoothing: " << std::endl;
366  for(int i=1;i<3;i++) {
367  for(int j=1;j<3;j++) {
368  for(int k=1;k<3;k++) {
369  std::cout<<i<<" "<<j<<" "<<k<<" : "<<input3D->GetBinContent(i+ioffset,j+joffset,k+koffset)<< " / " ;
370  }
371  }
372  }
373  std::cout << std::endl;
374  }
375  //
376  const int lsup = 20;
377  if (m1 > lsup || m2 > lsup || m3 > lsup) {
378  std::cout <<"HistoHelperRoot::smoothASH3D: m1, m2 or m3 too big !"<<std::endl;
379  return;
380  } else {
381  int nx = input3D->GetNbinsX()+1;
382  int ny = input3D->GetNbinsY()+1;
383  int nz = input3D->GetNbinsZ()+1;
384  float ***h, ***res;
385  h = new float**[nx-1];
386  res = new float**[nx-1];
387  for (int i = 0;i < nx-1; i++) {
388  h[i] = new float*[ny-1];
389  res[i] = new float*[ny-1];
390  for (int j = 0; j < ny-1; j++) {
391  h[i][j] = new float[nz-1];
392  res[i][j] = new float[nz-1];
393  }
394  }
395  for (int iz = 1;iz<nz;iz++) {
396  for (int iy = 1;iy<ny;iy++) {
397  for (int ix = 1;ix<nx;ix++) {
398  h[ix-1][iy-1][iz-1] = (float) input3D->GetBinContent(ix,iy,iz);
399  }
400  }
401  }
402  //
403  int i,j,k,l,m,n;
404  float wk1[41],wk2[41],wk3[41];
405  //float wgt[100][100][100]; // Trop gros pour certaines machines !!??
406  //avoid large stack use
407  auto wk = std::vector(41,std::vector(41,std::vector<float>(41)));
408  double wks = 0.;
409  float ai,am1 = float(m1), am2 = float(m2), am3 = float(m3);
410  const float am12 = am1*am1, am22 = am2*am2, am32 = am3*am3;
411  const float inv_am1_am2 = 1. / (am1*am2);
412  const float inv_am12 = 1. / am12;
413  const float inv_am22 = 1. / am22;
414  const float inv_am32 = 1. / am32;
415 
416  // Initialisation
417  for (k = 0;k<nx-1;k++) {
418  for (l = 0;l<ny-1;l++) {
419  for (m = 0;m<nz-1;m++) {
420  res[k][l][m] = 0.;
421  //wgt[k][l][m] = 0.;
422  }
423  }
424  }
425  //
426  // Bidouille car selon la machine un tableau wgt[100][100][100] peut poser probleme !
427  //
428  float ***wgt;
429  int dimension = 100; int x, y, z = 0;
430  wgt = new float**[dimension];
431  if (!wgt) {
432  std::cout <<"HistoHelperRoot::smoothASH3D: problem to allocate memory, exiting"<<std::endl;
433  std::abort();
434  }
435  for(x = 0; x < dimension; x++) {
436  wgt[x] = new float*[dimension];
437  if (!wgt) {
438  std::cout <<"HistoHelperRoot::smoothASH3D: problem to allocate memory, exiting"<<std::endl;
439  std::abort();
440  }
441  for(y = 0; y < dimension; y++) {
442  wgt[x][y] = new float[dimension];
443  if (!wgt) {
444  std::cout <<"HistoHelperRoot::smoothASH3D: problem to allocate memory, exiting"<<std::endl;
445  std::abort();
446  }
447  for(z = 0; z < dimension; z++)
448  wgt[x][y][z] = 0;
449  }
450  }
451  // Weights
452  for (i = lsup+1-m1;i<lsup+m1;i++) {
453  ai = float(i-lsup)*float(i-lsup);
454  wk1[i] = 15./16.*(1.-ai*inv_am12)*(1.-ai*inv_am12);
455  wks = wks + wk1[i];
456  }
457  const double fac1 = am1 / wks;
458  for (i = lsup+1-m1;i<lsup+m1;i++) {
459  wk1[i] = wk1[i]*fac1;
460  }
461  wks = 0.;
462  for (i = lsup+1-m2;i<lsup+m2;i++) {
463  ai = float(i-lsup)*float(i-lsup);
464  wk2[i] = 15./16.*(1.-ai*inv_am22)*(1.-ai*inv_am22);
465  wks = wks + wk2[i];
466  }
467  if (wks == 0){
468  std::cout <<"HistoHelperRoot::smoothASH3D: wks is zero! "<<std::endl;
469  return;
470  }
471  const double fac2 = am2 / wks;
472  for (i = lsup+1-m2;i<lsup+m2;i++) {
473  wk2[i] = wk2[i]*fac2;
474  }
475  wks = 0.;
476  for (i = lsup+1-m3;i<lsup+m3;i++) {
477  ai = float(i-lsup)*float(i-lsup);
478  wk3[i] = 15./16.*(1.-ai*inv_am32)*(1.-ai*inv_am32);
479  wks = wks + wk3[i];
480  }
481  if (wks == 0){
482  std::cout <<"HistoHelperRoot::smoothASH3D: wks is zero! "<<std::endl;
483  return;
484  }
485  const double fac3 = am3 / wks;
486  for (i = lsup+1-m3;i<lsup+m3;i++) {
487  wk3[i] = wk3[i]*fac3;
488  }
489 
490  for (i = lsup+1-m1;i<lsup+m1;i++) {
491  for (j = lsup+1-m2;j<lsup+m2;j++) {
492  for (k = lsup+1-m3;k<lsup+m3;k++) {
493  wk[i][j][k] = wk1[i]*wk2[j]*wk3[k];
494  }
495  }
496  }
497  //
498  for (k = 0;k<nx-1;k++) {
499  for (l = 0;l<ny-1;l++) {
500  for (m = 0;m<nz-1;m++) {
501  for (i = std::max(0,k-m1+1);i<std::min(nx-1,k+m1);i++) {
502  for (j = std::max(0,l-m2+1);j<std::min(ny-1,l+m2);j++) {
503  for (n = std::max(0,m-m3+1);n<std::min(nz-1,m+m3);n++) {
504  res[i][j][n] = res[i][j][n] + wk[lsup+i-k][lsup+j-l][lsup+n-m]*h[k][l][m];
505  wgt[i][j][n] = wgt[i][j][n] + wk[lsup+i-k][lsup+j-l][lsup+n-m];
506  }
507  }
508  }
509  }
510  }
511  }
512  for (k = 0;k<nx-1;k++) {
513  for (l = 0;l<ny-1;l++) {
514  for (m = 0;m<nz-1;m++) {
515  res[k][l][m] = res[k][l][m]*inv_am1_am2;
516  if (wgt[k][l][m] != 0.) {res[k][l][m] = res[k][l][m]/wgt[k][l][m];}
517  }
518  }
519  }
520  // replace the content of histo with results of smoothing:
521  input3D->Reset();
522  for (int iz = 1;iz<nz;iz++) {
523  for (int iy = 1;iy<ny;iy++) {
524  for (int ix = 1;ix<nx;ix++) {
525  input3D->SetBinContent(ix,iy,iz,res[ix-1][iy-1][iz-1]);
526  }
527  }
528  }
529  for (i = 0; i < nx-1; i++){
530  for (j = 0; j < ny-1; j++) {
531  delete[] h[i][j];
532  delete[] res[i][j];
533  }
534  delete[] h[i];
535  delete[] res[i];
536  }
537  delete[] h;
538  delete[] res;
539  //
540  for (i = 0; i < dimension; i++){
541  for (j = 0; j < dimension; j++)
542  delete[] wgt[i][j];
543  delete[] wgt[i];
544  }
545  delete[] wgt;
546  }
547  if(debug) {
548  std::cout << "First 2x2x2 bins after smoothing: " << std::endl;
549  for(int i=1;i<3;i++) {
550  for(int j=1;j<3;j++) {
551  for(int k=1;k<3;k++) {
552  std::cout<<i<<" "<<j<<" "<<k<<" : "<<input3D->GetBinContent(i,j,k)<< " / " ;
553  }
554  }
555  }
556  std::cout << std::endl;
557  int ioffset = input3D->GetNbinsX() / 2 ;
558  int joffset = input3D->GetNbinsY() / 2 ;
559  int koffset = input3D->GetNbinsZ() / 2 ;
560  std::cout << "Middle (" << ioffset+1 << "-" << ioffset+3 << ", "
561  << joffset+1 << "-" << joffset+3 << ", "
562  << koffset+1 << "-" << koffset+3 << ") 2x2 bins after smoothing: " << std::endl;
563  for(int i=1;i<3;i++) {
564  for(int j=1;j<3;j++) {
565  for(int k=1;k<3;k++) {
566  std::cout<<i<<" "<<j<<" "<<k<<" : "<<input3D->GetBinContent(i+ioffset,j+joffset,k+koffset)<< " / " ;
567  }
568  }
569  }
570  std::cout << std::endl;
571  }
572 }

Member Data Documentation

◆ m_checkOverflows

bool Analysis::HistoHelperRoot::m_checkOverflows
private

Definition at line 83 of file HistoHelperRoot.h.

◆ m_histoLimitsMap1D

std::map<std::string, HistoLimits> Analysis::HistoHelperRoot::m_histoLimitsMap1D
private

Definition at line 80 of file HistoHelperRoot.h.

◆ m_histoLimitsMap2D

std::map<std::string, HistoLimits> Analysis::HistoHelperRoot::m_histoLimitsMap2D
private

Definition at line 81 of file HistoHelperRoot.h.

◆ m_histoLimitsMap3D

std::map<std::string, HistoLimits> Analysis::HistoHelperRoot::m_histoLimitsMap3D
private

Definition at line 82 of file HistoHelperRoot.h.

◆ m_histoMap1D

std::map<std::string, LockedHandle<TH1> > Analysis::HistoHelperRoot::m_histoMap1D
private

Definition at line 77 of file HistoHelperRoot.h.

◆ m_histoMap2D

std::map<std::string, LockedHandle<TH2> > Analysis::HistoHelperRoot::m_histoMap2D
private

Definition at line 78 of file HistoHelperRoot.h.

◆ m_histoMap3D

std::map<std::string, LockedHandle<TH3> > Analysis::HistoHelperRoot::m_histoMap3D
private

Definition at line 79 of file HistoHelperRoot.h.

◆ m_histoSvc

ITHistSvc* Analysis::HistoHelperRoot::m_histoSvc
private

Definition at line 76 of file HistoHelperRoot.h.

The documentation for this class was generated from the following files:
base
std::string base
Definition: hcg.cxx:81
beamspotman.r
def r
Definition: beamspotman.py:672
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
bin
Definition: BinsDiffFromStripMedian.h:43
athena.value
value
Definition: athena.py:124
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:157
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
x
#define x
Trk::u
@ u
Enums for curvilinear frames.
Definition: ParamDefs.h:77
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:190
Analysis::HistoHelperRoot::m_histoLimitsMap1D
std::map< std::string, HistoLimits > m_histoLimitsMap1D
Definition: HistoHelperRoot.h:80
python.changerun.m1
m1
Definition: changerun.py:30
lumiFormat.i
int i
Definition: lumiFormat.py:85
z
#define z
beamspotman.n
n
Definition: beamspotman.py:727
extractSporadic.h
list h
Definition: extractSporadic.py:96
res
std::pair< std::vector< unsigned int >, bool > res
Definition: JetGroupProductTest.cxx:11
plotting.yearwise_luminosity_vs_mu.bins
bins
Definition: yearwise_luminosity_vs_mu.py:30
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
Analysis::HistoHelperRoot::baseName
std::string baseName(const std::string &fullHistoName)
Definition: HistoHelperRoot.cxx:36
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
debug
const bool debug
Definition: MakeUncertaintyPlots.cxx:53
PlotCalibFromCool.words
words
Definition: PlotCalibFromCool.py:51
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:82
Analysis::HistoHelperRoot::m_histoSvc
ITHistSvc * m_histoSvc
Definition: HistoHelperRoot.h:76
Analysis::HistoHelperRoot::m_checkOverflows
bool m_checkOverflows
Definition: HistoHelperRoot.h:83
unlikely
#define unlikely(x)
Definition: dictionary.h:41
python.PyAthena.v
v
Definition: PyAthena.py:154
Analysis::HistoHelperRoot::m_histoLimitsMap3D
std::map< std::string, HistoLimits > m_histoLimitsMap3D
Definition: HistoHelperRoot.h:82
y
#define y
h
Analysis::HistoHelperRoot::m_histoLimitsMap2D
std::map< std::string, HistoLimits > m_histoLimitsMap2D
Definition: HistoHelperRoot.h:81
Analysis::HistoHelperRoot::m_histoMap1D
std::map< std::string, LockedHandle< TH1 > > m_histoMap1D
Definition: HistoHelperRoot.h:77
python.SystemOfUnits.m3
float m3
Definition: SystemOfUnits.py:108
ATLAS_THREAD_SAFE
#define ATLAS_THREAD_SAFE
Definition: checker_macros.h:211
python.SystemOfUnits.m2
float m2
Definition: SystemOfUnits.py:107
Analysis::HistoHelperRoot::m_histoMap2D
std::map< std::string, LockedHandle< TH2 > > m_histoMap2D
Definition: HistoHelperRoot.h:78
Analysis::HistoHelperRoot::m_histoMap3D
std::map< std::string, LockedHandle< TH3 > > m_histoMap3D
Definition: HistoHelperRoot.h:79
fitman.k
k
Definition: fitman.py:528
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
python.LArMinBiasAlgConfig.float
float
Definition: LArMinBiasAlgConfig.py:65
Generated on Fri Nov 21 2025 21:53:34 for ATLAS Offline Software by   doxygen 1.8.18