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  std::string word = fullHistoName.substr(sPos,sLen);
46  words.push_back(word);
47  sPos = fullHistoName.find_first_not_of(delim, sEnd);
48  }
49  std::string base = "";
50  if(words.size()>0) base = words[words.size()-1];
51  //std::cout<<"baseName: decoding "<<fullHistoName<<" --> "<<base<<std::endl;
52  return base;
53 }

◆ bookHisto() [1/5]

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

Definition at line 65 of file HistoHelperRoot.cxx.

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

◆ 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 55 of file HistoHelperRoot.cxx.

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

◆ 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 83 of file HistoHelperRoot.cxx.

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

◆ 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 74 of file HistoHelperRoot.cxx.

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

◆ 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 92 of file HistoHelperRoot.cxx.

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

◆ 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 103 of file HistoHelperRoot.cxx.

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

◆ 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 131 of file HistoHelperRoot.cxx.

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

◆ 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 147 of file HistoHelperRoot.cxx.

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

◆ 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 117 of file HistoHelperRoot.cxx.

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

◆ Interpol1d()

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

Definition at line 545 of file HistoHelperRoot.cxx.

546 {
547  int bin((h->GetXaxis())->FindBin(x));
548  double bincenter(h->GetBinCenter(bin));
549  double nextbincenter(bincenter);
550  double nextbincontent(0);
551  if(x>bincenter && bin < h->GetNbinsX())
552  {
553  nextbincenter = h->GetBinCenter(bin+1);
554  nextbincontent = h->GetBinContent(bin+1);
555  }
556  else if(x<bincenter && bin > 1)
557  {
558  nextbincenter = h->GetBinCenter(bin-1);
559  nextbincontent = h->GetBinContent(bin-1);
560  }
561  double tmp(h->GetBinContent(bin));
562  if(nextbincenter != bincenter) tmp = ( nextbincontent*(x-bincenter) + tmp*(nextbincenter-x) ) / (nextbincenter - bincenter);
563  return tmp;
564 }

◆ Interpol2d()

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

Definition at line 569 of file HistoHelperRoot.cxx.

570 {
571  int nx = h->GetNbinsX(), ny = h->GetNbinsY();
572  int ibx = h->GetXaxis()->FindBin(x), iby = h->GetYaxis()->FindBin(y);
573  int ibx2,iby2;
574  double z00,z11,z01,z10,xc,yc,xc2,yc2,u,t,r;
575  if (ibx > nx) ibx = nx;
576  if (iby > ny) iby = ny;
577  xc = h->GetXaxis()->GetBinCenter(ibx);
578  yc = h->GetYaxis()->GetBinCenter(iby);
579  z11 = h->GetBinContent(ibx,iby);
580  if ( (ibx > 1 || (ibx == 1 && x > xc) ) &&
581  (ibx < nx || (ibx == nx && x < xc) ) ) {
582  if (x > xc) {ibx2 = ibx + 1;} else {ibx2 = ibx - 1;}
583  xc2 = h->GetXaxis()->GetBinCenter(ibx2);
584  if ( (iby > 1 || (iby == 1 && y > yc) ) &&
585  (iby < ny || (iby == ny && y < yc) ) ) {
586  if (y > yc) {iby2 = iby + 1;} else {iby2 = iby - 1;}
587  yc2 = h->GetYaxis()->GetBinCenter(iby2);
588  z00 = h->GetBinContent(ibx2,iby2);
589  z10 = h->GetBinContent(ibx,iby2);
590  z01 = h->GetBinContent(ibx2,iby);
591  t = (x - xc2)/(xc-xc2);
592  u = (y - yc2)/(yc-yc2);
593  r = z11*t*u+z00*(1.-t)*(1.-u)+z01*(1.-t)*u+z10*t*(1.-u);
594  } else {
595  z01 = h->GetBinContent(ibx2,iby);
596  t = (x - xc2)/(xc-xc2);
597  r = z11*t + z01*(1.-t);
598  }
599  } else if ((iby > 1 || (iby == 1 && y > yc) ) &&
600  (iby < ny || (iby == ny && y < yc) ) ) {
601  if (y > yc) {iby2 = iby + 1;} else {iby2 = iby - 1;}
602  z10 = h->GetBinContent(ibx,iby2);
603  yc2 = h->GetYaxis()->GetBinCenter(iby2);
604  u = (y - yc2)/(yc-yc2);
605  r = z11*u + z10*(1.-u);
606  } else {
607  r = z11;
608  }
609  return r;
610 }

◆ Interpol3d()

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

Definition at line 615 of file HistoHelperRoot.cxx.

616 {
617  int nx = h->GetNbinsX(), ny = h->GetNbinsY(), nz = h->GetNbinsZ();
618  int ibx = h->GetXaxis()->FindBin(x), iby = h->GetYaxis()->FindBin(y), ibz = h->GetZaxis()->FindBin(z);
619  int ibx2,iby2,ibz2;
620  double z000,z010,z110,z100,z001,z011,z111,z101,xc,yc,zc,xc2,yc2,zc2,u,t,v,r;
621  if (ibx > nx) ibx = nx;
622  if (iby > ny) iby = ny;
623  if (ibz > nz) ibz = nz;
624  xc = h->GetXaxis()->GetBinCenter(ibx);
625  yc = h->GetYaxis()->GetBinCenter(iby);
626  zc = h->GetZaxis()->GetBinCenter(ibz);
627  //
628  z111 = h->GetBinContent(ibx,iby,ibz);
629  if ( (ibx > 1 || (ibx == 1 && x > xc) ) &&
630  (ibx < nx || (ibx == nx && x < xc) ) ) {
631  if (x > xc) {ibx2 = ibx + 1;} else {ibx2 = ibx - 1;}
632  xc2 = h->GetXaxis()->GetBinCenter(ibx2);
633  if ( (iby > 1 || (iby == 1 && y > yc) ) &&
634  (iby < ny || (iby == ny && y < yc) ) ) {
635  if (y > yc) {iby2 = iby + 1;} else {iby2 = iby - 1;}
636  yc2 = h->GetYaxis()->GetBinCenter(iby2);
637  if ( (ibz > 1 || (ibz == 1 && z > zc) ) &&
638  (ibz < nz || (ibz == nz && z < zc) ) ) {
639  if (z > zc) {ibz2 = ibz + 1;} else {ibz2 = ibz - 1;}
640  zc2 = h->GetZaxis()->GetBinCenter(ibz2);
641  //
642  z000 = h->GetBinContent(ibx2,iby2,ibz2);
643  z100 = h->GetBinContent(ibx,iby2,ibz2);
644  z010 = h->GetBinContent(ibx2,iby,ibz2);
645  z110 = h->GetBinContent(ibx,iby,ibz2);
646  z001 = h->GetBinContent(ibx2,iby2,ibz);
647  z101 = h->GetBinContent(ibx,iby2,ibz);
648  z011 = h->GetBinContent(ibx2,iby,ibz);
649  t = (x - xc2)/(xc-xc2);
650  u = (y - yc2)/(yc-yc2);
651  v = (z - zc2)/(zc-zc2);
652  r = z111*t*u*v + z001*(1.-t)*(1.-u)*v
653  + z011*(1.-t)*u*v + z101*t*(1.-u)*v
654  + z110*t*u*(1.-v) + z000*(1.-t)*(1.-u)*(1.-v)
655  + z010*(1.-t)*u*(1.-v) + z100*t*(1.-u)*(1.-v);
656  } else {
657  z011 = h->GetBinContent(ibx2,iby,ibz);
658  z001 = h->GetBinContent(ibx2,iby2,ibz);
659  z101 = h->GetBinContent(ibx,iby2,ibz);
660  t = (x - xc2)/(xc-xc2);
661  u = (y - yc2)/(yc-yc2);
662  r = z111*t*u + z011*(1.-t)*u + z101*t*(1.-u) + z001*(1.-t)*(1.-u) ;
663  }
664  } else if ((ibz > 1 || (ibz == 1 && z > zc) ) &&
665  (ibz < nz || (ibz == nz && z < zc) ) ) {
666  if (z > zc) {ibz2 = ibz + 1;} else {ibz2 = ibz - 1;}
667  z110 = h->GetBinContent(ibx,iby,ibz2);
668  z010 = h->GetBinContent(ibx2,iby,ibz2);
669  z011 = h->GetBinContent(ibx2,iby,ibz);
670  zc2 = h->GetYaxis()->GetBinCenter(ibz2);
671  t = (x - xc2)/(xc-xc2);
672  v = (z - zc2)/(zc-zc2);
673  r = z111*t*v + z011*(1.-t)*v + z110*t*(1.-v) +z010*(1.-t)*(1.-v);
674  } else {
675  z011 = h->GetBinContent(ibx2,iby,ibz);
676  t = (x - xc2)/(xc-xc2);
677  r = z111*t + z011*(1.-t);
678  }
679  } else {
680  if ( (iby > 1 || (iby == 1 && y > yc) ) &&
681  (iby < ny || (iby == ny && y < yc) ) ) {
682  if (y > yc) {iby2 = iby + 1;} else {iby2 = iby - 1;}
683  yc2 = h->GetYaxis()->GetBinCenter(iby2);
684  if ( (ibz > 1 || (ibz == 1 && z > zc) ) &&
685  (ibz < nz || (ibz == nz && z < zc) ) ) {
686  if (z > zc) {ibz2 = ibz + 1;} else {ibz2 = ibz - 1;}
687  zc2 = h->GetZaxis()->GetBinCenter(ibz2);
688  //
689  z100 = h->GetBinContent(ibx,iby2,ibz2);
690  z110 = h->GetBinContent(ibx,iby,ibz2);
691  z101 = h->GetBinContent(ibx,iby2,ibz);
692  u = (y - yc2)/(yc-yc2);
693  v = (z - zc2)/(zc-zc2);
694  r = z111*u*v + z101*(1.-u)*v
695  + z110*u*(1.-v) + z100*(1.-u)*(1.-v);
696  } else {
697  z101 = h->GetBinContent(ibx,iby2,ibz);
698  u = (y - yc2)/(yc-yc2);
699  r = z111*u + z101*(1.-u);
700  }
701  } else if ((ibz > 1 || (ibz == 1 && z > zc) ) &&
702  (ibz < nz || (ibz == nz && z < zc) ) ) {
703  if (z > zc) {ibz2 = ibz + 1;} else {ibz2 = ibz - 1;}
704  z110 = h->GetBinContent(ibx,iby,ibz2);
705  zc2 = h->GetYaxis()->GetBinCenter(ibz2);
706  v = (z - zc2)/(zc-zc2);
707  r = z111*v + z110*(1.-v);
708  } else {
709  r = z111;
710  }
711  }
712  return r;
713 }

◆ print()

void Analysis::HistoHelperRoot::print ( )

Definition at line 178 of file HistoHelperRoot.cxx.

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

◆ 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 190 of file HistoHelperRoot.cxx.

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

◆ smoothASH3D()

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

Definition at line 324 of file HistoHelperRoot.cxx.

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

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:78
beamspotman.r
def r
Definition: beamspotman.py:676
python.SystemOfUnits.m2
int m2
Definition: SystemOfUnits.py:92
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
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:158
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:189
Analysis::HistoHelperRoot::m_histoLimitsMap1D
std::map< std::string, HistoLimits > m_histoLimitsMap1D
Definition: HistoHelperRoot.h:80
python.changerun.m1
m1
Definition: changerun.py:32
lumiFormat.i
int i
Definition: lumiFormat.py:85
z
#define z
beamspotman.n
n
Definition: beamspotman.py:731
extractSporadic.h
list h
Definition: extractSporadic.py:97
res
std::pair< std::vector< unsigned int >, bool > res
Definition: JetGroupProductTest.cxx:14
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:228
debug
const bool debug
Definition: MakeUncertaintyPlots.cxx:53
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:83
Analysis::HistoHelperRoot::m_histoSvc
ITHistSvc * m_histoSvc
Definition: HistoHelperRoot.h:76
Analysis::HistoHelperRoot::m_checkOverflows
bool m_checkOverflows
Definition: HistoHelperRoot.h:83
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
ATLAS_THREAD_SAFE
#define ATLAS_THREAD_SAFE
Definition: checker_macros.h:211
checkFileSG.words
words
Definition: checkFileSG.py:76
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
readCCLHist.float
float
Definition: readCCLHist.py:83
python.SystemOfUnits.m3
int m3
Definition: SystemOfUnits.py:93
fitman.k
k
Definition: fitman.py:528
Generated on Sun Dec 22 2024 21:55:55 for ATLAS Offline Software by   doxygen 1.8.18