The Gaudi Framework  master (d98a2936)
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
H2D.cpp
Go to the documentation of this file.
1 /***********************************************************************************\
2 * (c) Copyright 1998-2025 CERN for the benefit of the LHCb and ATLAS collaborations *
3 * *
4 * This software is distributed under the terms of the Apache version 2 licence, *
5 * copied verbatim in the file "LICENSE". *
6 * *
7 * In applying this licence, CERN does not waive the privileges and immunities *
8 * granted to it by virtue of its status as an Intergovernmental Organization *
9 * or submit itself to any jurisdiction. *
10 \***********************************************************************************/
11 #ifdef __ICC
12 // disable icc remark #2259: non-pointer conversion from "X" to "Y" may lose significant bits
13 // TODO: To be removed, since it comes from ROOT TMathBase.h
14 # pragma warning( disable : 2259 )
15 #endif
16 
17 #include "GaudiPI.h"
18 #include <Gaudi/MonitoringHub.h>
19 #include <GaudiCommonSvc/H1D.h>
20 #include <GaudiCommonSvc/H2D.h>
21 #include <GaudiCommonSvc/HistogramUtility.h>
22 #include <GaudiCommonSvc/P1D.h>
23 #include <GaudiKernel/ObjectFactory.h>
24 #include <TH2D.h>
25 #include <TProfile.h>
26 #include <array>
27 #include <cmath>
28 
29 namespace {
30  using AIDA::IHistogram1D;
31  using AIDA::IHistogram2D;
32  using AIDA::IProfile1D;
33 } // namespace
34 
35 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( ISvcLocator* svcLocator, const std::string& path,
36  const std::string& title, int binsX, double iminX, double imaxX,
37  int binsY, double iminY, double imaxY ) {
38  auto p = new Histogram2D( new TH2D( title.c_str(), title.c_str(), binsX, iminX, imaxX, binsY, iminY, imaxY ) );
39  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
40  return { p, p };
41 }
42 
43 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( ISvcLocator* svcLocator, const std::string& path,
44  const std::string& title, const Edges& eX, const Edges& eY ) {
45  auto p = new Histogram2D(
46  new TH2D( title.c_str(), title.c_str(), eX.size() - 1, &eX.front(), eY.size() - 1, &eY.front() ) );
47  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
48  return { p, p };
49 }
50 
51 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( ISvcLocator* svcLocator, const std::string& path, TH2D* rep ) {
52  auto p = new Histogram2D( rep );
53  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
54  return { p, p };
55 }
56 
57 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( ISvcLocator* svcLocator, const std::string& path,
58  const IHistogram2D& hist ) {
59  TH2D* h = getRepresentation<AIDA::IHistogram2D, TH2D>( hist );
60  Histogram2D* n = h ? new Histogram2D( new TH2D( *h ) ) : nullptr;
61  if ( n ) { svcLocator->monitoringHub().registerEntity( path, h->GetName(), "histogram:Histogram:double", *n ); }
62  return { n, n };
63 }
64 
65 std::pair<DataObject*, IHistogram1D*> Gaudi::slice1DX( const std::string& nam, const IHistogram2D& hist, int first,
66  int last ) {
67  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
68  TH1D* t = r ? r->ProjectionX( "_px", first, last, "e" ) : nullptr;
69  if ( t ) t->SetName( nam.c_str() );
70  Histogram1D* p = ( t ? new Histogram1D( t ) : nullptr );
71  return { p, p };
72 }
73 
74 std::pair<DataObject*, IHistogram1D*> Gaudi::slice1DY( const std::string& nam, const IHistogram2D& hist, int first,
75  int last ) {
76  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
77  TH1D* t = r ? r->ProjectionY( "_py", first, last, "e" ) : nullptr;
78  if ( t ) t->SetName( nam.c_str() );
79  Histogram1D* p = ( t ? new Histogram1D( t ) : nullptr );
80  return { p, p };
81 }
82 
83 std::pair<DataObject*, IHistogram1D*> Gaudi::project1DY( const std::string& nam, const IHistogram2D& hist, int first,
84  int last ) {
85  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
86  TH1D* t = r ? r->ProjectionY( "_px", first, last, "e" ) : nullptr;
87  if ( t ) t->SetName( nam.c_str() );
88  Histogram1D* p = ( t ? new Histogram1D( t ) : nullptr );
89  return { p, p };
90 }
91 
92 std::pair<DataObject*, IProfile1D*> Gaudi::profile1DX( const std::string& nam, const IHistogram2D& hist, int first,
93  int last ) {
94  TH2* r = Gaudi::getRepresentation<IHistogram2D, TH2>( hist );
95  TProfile* t = r ? r->ProfileX( "_pfx", first, last, "e" ) : nullptr;
96  if ( t ) t->SetName( nam.c_str() );
97  Profile1D* p = ( t ? new Profile1D( t ) : nullptr );
98  return { p, p };
99 }
100 
101 std::pair<DataObject*, IProfile1D*> Gaudi::profile1DY( const std::string& nam, const IHistogram2D& hist, int first,
102  int last ) {
103  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
104  TProfile* t = r ? r->ProfileY( "_pfx", first, last, "e" ) : nullptr;
105  if ( t ) t->SetName( nam.c_str() );
106  Profile1D* p = ( t ? new Profile1D( t ) : nullptr );
107  return { p, p };
108 }
109 
110 namespace Gaudi {
111  template <>
112  void* Generic2D<IHistogram2D, TH2D>::cast( const std::string& className ) const {
113  if ( className == "AIDA::IHistogram2D" )
114  return (IHistogram2D*)this;
115  else if ( className == "AIDA::IHistogram" )
116  return (IHistogram*)this;
117  return nullptr;
118  }
119 
120  template <>
121  int Generic2D<IHistogram2D, TH2D>::binEntries( int indexX, int indexY ) const {
122  if ( binHeight( indexX, indexY ) <= 0 ) return 0;
123  double xx = binHeight( indexX, indexY ) / binError( indexX, indexY );
124  return int( xx * xx + 0.5 );
125  }
126 
127  template <>
128  void Generic2D<AIDA::IHistogram2D, TH2D>::adoptRepresentation( TObject* rep ) {
129  TH2D* imp = dynamic_cast<TH2D*>( rep );
130  if ( !imp ) throw std::runtime_error( "Cannot adopt native histogram representation." );
131  m_rep.reset( imp );
132  m_xAxis.initialize( m_rep->GetXaxis(), true );
133  m_yAxis.initialize( m_rep->GetYaxis(), true );
134  const TArrayD* a = m_rep->GetSumw2();
135  if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
136  setTitle( m_rep->GetTitle() );
137  }
138 } // namespace Gaudi
139 
140 Gaudi::Histogram2D::Histogram2D() : Base( new TH2D() ) {
141  m_rep->Sumw2();
142  m_sumwx = m_sumwy = 0;
143  setTitle( "" );
144  m_rep->SetDirectory( nullptr );
145 }
146 
147 Gaudi::Histogram2D::Histogram2D( TH2D* rep ) {
148  adoptRepresentation( rep );
149  m_sumwx = m_sumwy = 0;
150  m_rep->SetDirectory( nullptr );
151 }
152 
153 bool Gaudi::Histogram2D::setBinContents( int i, int j, int entries, double height, double error, double centreX,
154  double centreY ) {
155  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), height );
156  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), error );
157  // accumulate sumwx for in range bins
158  if ( i >= 0 && j >= 0 ) {
159  m_sumwx += centreX * height;
160  m_sumwy += centreY * height;
161  }
162  m_sumEntries += entries;
163  return true;
164 }
165 
166 bool Gaudi::Histogram2D::reset() {
167  m_sumwx = 0;
168  m_sumwy = 0;
169  return Base::reset();
170 }
171 
172 bool Gaudi::Histogram2D::fill( double x, double y, double weight ) {
173  // avoid race conditiosn when filling the histogram
174  auto guard = std::scoped_lock{ m_fillSerialization };
175  m_rep->Fill( x, y, weight );
176  return true;
177 }
178 
179 bool Gaudi::Histogram2D::setRms( double rmsX, double rmsY ) {
180  m_rep->SetEntries( m_sumEntries );
181  std::vector<double> stat( 11 );
182  stat[0] = sumBinHeights();
183  stat[1] = 0;
184  if ( std::abs( equivalentBinEntries() ) > std::numeric_limits<double>::epsilon() )
185  stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
186  stat[2] = m_sumwx;
187  stat[4] = m_sumwy;
188  double meanX = 0;
189  double meanY = 0;
190  if ( std::abs( sumBinHeights() ) > std::numeric_limits<double>::epsilon() ) {
191  meanX = m_sumwx / sumBinHeights();
192  meanY = m_sumwy / sumBinHeights();
193  }
194  stat[3] = ( meanX * meanX + rmsX * rmsX ) * sumBinHeights();
195  stat[5] = ( meanY * meanY + rmsY * rmsY ) * sumBinHeights();
196  stat[6] = 0;
197  m_rep->PutStats( &stat.front() );
198  return true;
199 }
200 
201 void Gaudi::Histogram2D::copyFromAida( const IHistogram2D& h ) {
202  // implement here the copy
203  std::string titlestr = h.title();
204  const char* title = titlestr.c_str();
205  if ( h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning() )
206  m_rep.reset( new TH2D( title, title, h.xAxis().bins(), h.xAxis().lowerEdge(), h.xAxis().upperEdge(),
207  h.yAxis().bins(), h.yAxis().lowerEdge(), h.yAxis().upperEdge() ) );
208  else {
209  Edges eX, eY;
210  for ( int i = 0; i < h.xAxis().bins(); ++i ) eX.push_back( h.xAxis().binLowerEdge( i ) );
211  // add also upperedges at the end
212  eX.push_back( h.xAxis().upperEdge() );
213  for ( int i = 0; i < h.yAxis().bins(); ++i ) eY.push_back( h.yAxis().binLowerEdge( i ) );
214  // add also upperedges at the end
215  eY.push_back( h.yAxis().upperEdge() );
216  m_rep.reset( new TH2D( title, title, eX.size() - 1, &eX.front(), eY.size() - 1, &eY.front() ) );
217  }
218  m_xAxis.initialize( m_rep->GetXaxis(), true );
219  m_yAxis.initialize( m_rep->GetYaxis(), true );
220  m_rep->Sumw2();
221  m_sumEntries = 0;
222  m_sumwx = 0;
223  m_sumwy = 0;
224  // statistics
225  double sumw = h.sumBinHeights();
226  double sumw2 = 0;
227  if ( std::abs( h.equivalentBinEntries() ) > std::numeric_limits<double>::epsilon() )
228  sumw2 = ( sumw * sumw ) / h.equivalentBinEntries();
229  double sumwx = h.meanX() * h.sumBinHeights();
230  double sumwx2 = ( h.meanX() * h.meanX() + h.rmsX() * h.rmsX() ) * h.sumBinHeights();
231  double sumwy = h.meanY() * h.sumBinHeights();
232  double sumwy2 = ( h.meanY() * h.meanY() + h.rmsY() * h.rmsY() ) * h.sumBinHeights();
233  double sumwxy = 0;
234 
235  // copy the contents in (AIDA underflow/overflow are -2,-1)
236  for ( int i = -2; i < xAxis().bins(); ++i ) {
237  for ( int j = -2; j < yAxis().bins(); ++j ) {
238  // root binning starts from one !
239  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), h.binHeight( i, j ) );
240  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), h.binError( i, j ) );
241  // calculate statistics
242  if ( i >= 0 && j >= 0 ) { sumwxy += h.binHeight( i, j ) * h.binMeanX( i, j ) * h.binMeanY( i, j ); }
243  }
244  }
245  // need to do set entries after setting contents otherwise root will recalculate them
246  // taking into account how many time SetBinContents() has been called
247  m_rep->SetEntries( h.allEntries() );
248  // fill stat vector
249  std::array<double, 11> stat = { { sumw, sumw2, sumwx, sumwx2, sumwy, sumwy2, sumwxy } };
250  m_rep->PutStats( stat.data() );
251 }
Gaudi::Histogram2D::reset
bool reset() override
Definition: H2D.cpp:166
Gaudi::profile1DX
std::pair< DataObject *, AIDA::IProfile1D * > profile1DX(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D profile in X from 2D histogram.
Gaudi::profile1DY
std::pair< DataObject *, AIDA::IProfile1D * > profile1DY(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D profile in Y from 2D histogram.
Gaudi::Histogram2D::m_sumwx
double m_sumwx
Definition: H2D.h:58
ISvcLocator::monitoringHub
Gaudi::Monitoring::Hub & monitoringHub()
Definition: ISvcLocator.cpp:24
AtlasMCRecoFullPrecedenceDump.path
path
Definition: AtlasMCRecoFullPrecedenceDump.py:49
MonitoringHub.h
Gaudi::Histogram2D::m_sumwy
double m_sumwy
Definition: H2D.h:59
ISvcLocator
Definition: ISvcLocator.h:42
P1D.h
Gaudi::slice1DY
std::pair< DataObject *, AIDA::IHistogram1D * > slice1DY(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D slice from 2D histogram.
ObjectFactory.h
Gaudi::Generic2D< AIDA::IHistogram2D, TH2D >::setTitle
bool setTitle(const std::string &title) override
Set the title of the object.
Definition: Generic2D.h:169
Gaudi::Generic2D::cast
void * cast(const std::string &className) const override
Introspection method.
HistogramUtility.h
Gaudi::Generic2D::binEntries
int binEntries(int indexX, int indexY) const override
The number of entries (ie the number of times fill was called for this bin).
bug_34121.t
t
Definition: bug_34121.py:31
Gaudi::svcLocator
GAUDI_API ISvcLocator * svcLocator()
Gaudi::Generic2D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
Gaudi::Histogram2D::setBinContents
virtual bool setBinContents(int binIndexX, int binIndexY, int entries, double height, double error, double centreX, double centreY)
Fast filling method for a given bin. It can be also the over/underflow bin.
Definition: H2D.cpp:153
ProduceConsume.j
j
Definition: ProduceConsume.py:104
AlgSequencer.h
h
Definition: AlgSequencer.py:31
GaudiPython.Bindings.nullptr
nullptr
Definition: Bindings.py:87
Gaudi::Monitoring::Hub::registerEntity
void registerEntity(std::string c, std::string n, std::string t, T &ent)
Definition: MonitoringHub.h:137
H1D.h
Gaudi::Histogram2D::copyFromAida
void copyFromAida(const AIDA::IHistogram2D &h)
Create new histogram from any AIDA based histogram.
Definition: H2D.cpp:201
Gaudi
This file provides a Grammar for the type Gaudi::Accumulators::Axis It allows to use that type from p...
Definition: __init__.py:1
cpluginsvc.n
n
Definition: cpluginsvc.py:234
H2D.h
Gaudi::Histogram2D::Histogram2D
Histogram2D()
Standard Constructor.
Definition: H2D.cpp:140
Gaudi::Generic2D< AIDA::IHistogram2D, TH2D >::m_rep
std::unique_ptr< TH2D > m_rep
Reference to underlying implementation.
Definition: Generic2D.h:161
GaudiPI.h
Gaudi::Edges
std::vector< double > Edges
Definition: GaudiPI.h:27
Gaudi::slice1DX
std::pair< DataObject *, AIDA::IHistogram1D * > slice1DX(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D slice from 2D histogram.
Gaudi::project1DY
std::pair< DataObject *, AIDA::IHistogram1D * > project1DY(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D projection in Y from 2D histogram.
Gaudi::Generic2D< AIDA::IHistogram2D, TH2D >
Gaudi::Histogram2D::setRms
bool setRms(double rmsX, double rmsY)
Sets the rms of the histogram.
Definition: H2D.cpp:179
Gaudi::createH2D
std::pair< DataObject *, AIDA::IHistogram2D * > createH2D(ISvcLocator *svcLocator, const std::string &path, const AIDA::IHistogram2D &hist)
Copy constructor.
Gaudi::Histogram2D::fill
bool fill(double x, double y, double weight=1.) override
Fill the Histogram2D with a value and the.
Definition: H2D.cpp:172