The Gaudi Framework  master (37c0b60a)
H3D.cpp
Go to the documentation of this file.
1 /***********************************************************************************\
2 * (c) Copyright 1998-2024 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 #ifdef WIN32
17 // Disable warning
18 // warning C4996: 'sprintf': This function or variable may be unsafe.
19 // coming from TString.h
20 # pragma warning( disable : 4996 )
21 #endif
22 
23 #include <Gaudi/MonitoringHub.h>
24 #include <GaudiCommonSvc/Generic3D.h>
25 #include <GaudiCommonSvc/H3D.h>
26 #include <GaudiCommonSvc/HistogramUtility.h>
27 #include <GaudiKernel/DataObject.h>
28 #include <GaudiKernel/ObjectFactory.h>
29 
30 #include "GaudiPI.h"
31 #include <TH3D.h>
32 
33 #include <cmath>
34 
35 namespace Gaudi {
36  template <>
37  void Generic3D<AIDA::IHistogram3D, TH3D>::adoptRepresentation( TObject* rep ) {
38  TH3D* imp = dynamic_cast<TH3D*>( rep );
39  if ( !imp ) throw std::runtime_error( "Cannot adopt native histogram representation." );
40  m_rep.reset( imp );
41  m_xAxis.initialize( m_rep->GetXaxis(), true );
42  m_yAxis.initialize( m_rep->GetYaxis(), true );
43  m_zAxis.initialize( m_rep->GetZaxis(), true );
44  const TArrayD* a = m_rep->GetSumw2();
45  if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
46  setTitle( m_rep->GetTitle() );
47  }
48 } // namespace Gaudi
49 
51 std::pair<DataObject*, AIDA::IHistogram3D*> Gaudi::createH3D( ISvcLocator* svcLocator, const std::string& path,
52  const std::string& title, int nBinsX, double xlow,
53  double xup, int nBinsY, double ylow, double yup,
54  int nBinsZ, double zlow, double zup ) {
55  auto p = new Histogram3D(
56  new TH3D( title.c_str(), title.c_str(), nBinsX, xlow, xup, nBinsY, ylow, yup, nBinsZ, zlow, zup ) );
57  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
58  return { p, p };
59 }
60 
62 std::pair<DataObject*, AIDA::IHistogram3D*> Gaudi::createH3D( ISvcLocator* svcLocator, const std::string& path,
63  const std::string& title, const Edges& eX,
64  const Edges& eY, const Edges& eZ ) {
65  auto p = new Histogram3D( new TH3D( title.c_str(), title.c_str(), eX.size() - 1, &eX.front(), eY.size() - 1,
66  &eY.front(), eZ.size() - 1, &eZ.front() ) );
67  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
68  return { p, p };
69 }
70 
71 std::pair<DataObject*, AIDA::IHistogram3D*> Gaudi::createH3D( ISvcLocator* svcLocator, const std::string& path,
72  const AIDA::IHistogram3D& hist ) {
73  TH3D* h = getRepresentation<AIDA::IHistogram3D, TH3D>( hist );
74  Histogram3D* n = h ? new Histogram3D( new TH3D( *h ) ) : nullptr;
75  if ( n ) { svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *n ); }
76  return { n, n };
77 }
78 
79 Gaudi::Histogram3D::Histogram3D() : Base( new TH3D() ) {
80  setTitle( "" );
81  m_rep->Sumw2();
82  m_sumwx = 0;
83  m_sumwy = 0;
84  m_sumwz = 0;
85  m_rep->SetDirectory( nullptr );
86 }
87 
88 Gaudi::Histogram3D::Histogram3D( TH3D* rep ) {
89  adoptRepresentation( rep );
90  m_sumwx = 0;
91  m_sumwy = 0;
92  m_sumwz = 0;
93  m_rep->SetDirectory( nullptr );
94 }
95 
96 // set bin content (entries and centre are not used )
97 bool Gaudi::Histogram3D::setBinContents( int i, int j, int k, int entries, double height, double error, double centreX,
98  double centreY, double centreZ ) {
99  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), rIndexZ( k ), height );
100  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), rIndexZ( k ), error );
101  // accumulate sum bin centers
102  if ( i >= 0 && j >= 0 && k >= 0 ) {
103  m_sumwx += centreX * height;
104  m_sumwy += centreY * height;
105  m_sumwz += centreZ * height;
106  }
107  m_sumEntries += entries;
108  return true;
109 }
110 
111 bool Gaudi::Histogram3D::reset() {
112  m_sumwx = 0;
113  m_sumwy = 0;
114  m_sumwz = 0;
115  m_sumEntries = 0;
116  m_rep->Reset();
117  return true;
118 }
119 
120 bool Gaudi::Histogram3D::fill( double x, double y, double z, double weight ) {
121  // avoid race conditiosn when filling the histogram
122  auto guard = std::scoped_lock{ m_fillSerialization };
123  m_rep->Fill( x, y, z, weight );
124  return true;
125 }
126 
127 void* Gaudi::Histogram3D::cast( const std::string& className ) const {
128  if ( className == "AIDA::IHistogram3D" ) {
129  return (AIDA::IHistogram3D*)this;
130  } else if ( className == "AIDA::IHistogram" ) {
131  return (AIDA::IHistogram*)this;
132  }
133  return nullptr;
134 }
135 
136 bool Gaudi::Histogram3D::setRms( double rmsX, double rmsY, double rmsZ ) {
137  m_rep->SetEntries( m_sumEntries );
138  std::vector<double> stat( 11 );
139  // sum weights
140  stat[0] = sumBinHeights();
141  stat[1] = 0;
142  if ( std::abs( equivalentBinEntries() ) > std::numeric_limits<double>::epsilon() )
143  stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
144  stat[2] = m_sumwx;
145  stat[4] = m_sumwy;
146  stat[6] = 0;
147  stat[7] = m_sumwz;
148  double meanX = 0;
149  double meanY = 0;
150  double meanZ = 0;
151  if ( std::abs( sumBinHeights() ) > std::numeric_limits<double>::epsilon() ) {
152  meanX = m_sumwx / sumBinHeights();
153  meanY = m_sumwy / sumBinHeights();
154  meanZ = m_sumwz / sumBinHeights();
155  }
156  stat[3] = ( meanX * meanX + rmsX * rmsX ) * sumBinHeights();
157  stat[5] = ( meanY * meanY + rmsY * rmsY ) * sumBinHeights();
158  stat[8] = ( meanZ * meanZ + rmsZ * rmsZ ) * sumBinHeights();
159  // do not need to use sumwxy sumwxz and sumwyz
160  m_rep->PutStats( &stat.front() );
161  return true;
162 }
163 
164 void Gaudi::Histogram3D::copyFromAida( const AIDA::IHistogram3D& h ) {
165  // implement here the copy
166  std::string titlestr = h.title();
167  const char* title = titlestr.c_str();
168  if ( h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning() && h.zAxis().isFixedBinning() ) {
169  m_rep = std::make_unique<TH3D>( title, title, h.xAxis().bins(), h.xAxis().lowerEdge(), h.xAxis().upperEdge(),
170  h.yAxis().bins(), h.yAxis().lowerEdge(), h.yAxis().upperEdge(), h.zAxis().bins(),
171  h.zAxis().lowerEdge(), h.zAxis().upperEdge() );
172  } else {
173  Edges eX, eY, eZ;
174  for ( int i = 0; i < h.xAxis().bins(); ++i ) eX.push_back( h.xAxis().binLowerEdge( i ) );
175  // add also upperedges at the end
176  eX.push_back( h.xAxis().upperEdge() );
177  for ( int i = 0; i < h.yAxis().bins(); ++i ) eY.push_back( h.yAxis().binLowerEdge( i ) );
178  // add also upperedges at the end
179  eY.push_back( h.yAxis().upperEdge() );
180  for ( int i = 0; i < h.zAxis().bins(); ++i ) eZ.push_back( h.zAxis().binLowerEdge( i ) );
181  // add also upperedges at the end
182  eZ.push_back( h.zAxis().upperEdge() );
183  m_rep.reset(
184  new TH3D( title, title, eX.size() - 1, &eX.front(), eY.size() - 1, &eY.front(), eZ.size() - 1, &eZ.front() ) );
185  }
186  m_xAxis.initialize( m_rep->GetXaxis(), true );
187  m_yAxis.initialize( m_rep->GetYaxis(), true );
188  m_zAxis.initialize( m_rep->GetZaxis(), true );
189  const TArrayD* a = m_rep->GetSumw2();
190  if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
191  m_sumEntries = 0;
192  m_sumwx = 0;
193  m_sumwy = 0;
194  m_sumwz = 0;
195 
196  // statistics
197  double sumw = h.sumBinHeights();
198  double sumw2 = 0;
199  if ( std::abs( h.equivalentBinEntries() ) > std::numeric_limits<double>::epsilon() )
200  sumw2 = ( sumw * sumw ) / h.equivalentBinEntries();
201  double sumwx = h.meanX() * h.sumBinHeights();
202  double sumwx2 = ( h.meanX() * h.meanX() + h.rmsX() * h.rmsX() ) * h.sumBinHeights();
203  double sumwy = h.meanY() * h.sumBinHeights();
204  double sumwy2 = ( h.meanY() * h.meanY() + h.rmsY() * h.rmsY() ) * h.sumBinHeights();
205  double sumwz = h.meanZ() * h.sumBinHeights();
206  double sumwz2 = ( h.meanZ() * h.meanZ() + h.rmsZ() * h.rmsZ() ) * h.sumBinHeights();
207  double sumwxy = 0;
208  double sumwxz = 0;
209  double sumwyz = 0;
210 
211  // copy the contents in (AIDA underflow/overflow are -2,-1)
212  for ( int i = -2; i < xAxis().bins(); ++i ) {
213  for ( int j = -2; j < yAxis().bins(); ++j ) {
214  for ( int k = -2; k < zAxis().bins(); ++k ) {
215  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), rIndexZ( k ), h.binHeight( i, j, k ) );
216  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), rIndexZ( k ), h.binError( i, j, k ) );
217  // calculate statistics
218  if ( i >= 0 && j >= 0 && k >= 0 ) {
219  sumwxy += h.binHeight( i, j, k ) * h.binMeanX( i, j, k ) * h.binMeanY( i, j, k );
220  sumwxz += h.binHeight( i, j, k ) * h.binMeanX( i, j, k ) * h.binMeanZ( i, j, k );
221  sumwyz += h.binHeight( i, j, k ) * h.binMeanY( i, j, k ) * h.binMeanZ( i, j, k );
222  }
223  }
224  }
225  }
226  // need to do set entries after setting contents otherwise root will recalulate them
227  // taking into account how many time SetBinContents() has been called
228  m_rep->SetEntries( h.allEntries() );
229 
230  // fill stat vector
231  std::vector<double> stat( 11 );
232  stat[0] = sumw;
233  stat[1] = sumw2;
234  stat[2] = sumwx;
235  stat[3] = sumwx2;
236  stat[4] = sumwy;
237  stat[5] = sumwy2;
238  stat[6] = sumwxy;
239  stat[7] = sumwz;
240  stat[8] = sumwz2;
241  stat[9] = sumwxz;
242  stat[10] = sumwyz;
243  m_rep->PutStats( &stat.front() );
244 }
Gaudi::Histogram3D::setRms
virtual bool setRms(double rmsX, double rmsY, double rmsZ)
Sets the rms of the histogram.
Definition: H3D.cpp:136
ISvcLocator::monitoringHub
Gaudi::Monitoring::Hub & monitoringHub()
Definition: ISvcLocator.cpp:25
std::string
STL class.
Gaudi::Histogram3D::m_sumwx
double m_sumwx
Definition: H3D.h:59
H3D.h
Gaudi::Generic3D< AIDA::IHistogram3D, TH3D >::m_rep
std::unique_ptr< TH3D > m_rep
Reference to underlying implementation.
Definition: Generic3D.h:230
AtlasMCRecoFullPrecedenceDump.path
path
Definition: AtlasMCRecoFullPrecedenceDump.py:49
MonitoringHub.h
Gaudi::Histogram3D::cast
void * cast(const std::string &className) const override
Introspection method.
Definition: H3D.cpp:127
std::pair
std::vector< double >
ISvcLocator
Definition: ISvcLocator.h:46
Gaudi::Histogram3D::setBinContents
virtual bool setBinContents(int i, int j, int k, int entries, double height, double error, double centreX, double centreY, double centreZ)
Fast filling method for a given bin. It can be also the over/underflow bin.
Definition: H3D.cpp:97
Gaudi::Histogram3D::Histogram3D
Histogram3D()
Standard Constructor.
Definition: H3D.cpp:79
Gaudi::Histogram3D::m_sumwz
double m_sumwz
Definition: H3D.h:61
std::abs
Gaudi::ParticleID abs(const Gaudi::ParticleID &p)
Return the absolute value for a PID.
Definition: ParticleID.h:191
ObjectFactory.h
std::vector::front
T front(T... args)
Gaudi::Histogram3D::fill
bool fill(double x, double y, double z, double weight) override
Fill bin content.
Definition: H3D.cpp:120
HistogramUtility.h
std::vector::push_back
T push_back(T... args)
Gaudi::svcLocator
GAUDI_API ISvcLocator * svcLocator()
Gaudi::Histogram3D::copyFromAida
void copyFromAida(const AIDA::IHistogram3D &h)
Create new histogram from any AIDA based histogram.
Definition: H3D.cpp:164
ProduceConsume.j
j
Definition: ProduceConsume.py:104
std::string::c_str
T c_str(T... args)
AlgSequencer.h
h
Definition: AlgSequencer.py:31
Gaudi::createH3D
std::pair< DataObject *, AIDA::IHistogram3D * > createH3D(ISvcLocator *svcLocator, const std::string &path, const AIDA::IHistogram3D &hist)
Copy constructor.
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:140
std::runtime_error
STL class.
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
std::numeric_limits::epsilon
T epsilon(T... args)
Gaudi::Generic3D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
Gaudi::Histogram3D::m_sumwy
double m_sumwy
Definition: H3D.h:60
DataObject.h
GaudiPI.h
Gaudi::Edges
std::vector< double > Edges
Definition: GaudiPI.h:28
Gaudi::Histogram3D::reset
bool reset() override
Definition: H3D.cpp:111
Gaudi::Generic3D< AIDA::IHistogram3D, TH3D >::setTitle
bool setTitle(const std::string &title) override
Set the title of the object.
Definition: Generic3D.h:238
Generic3D.h
Gaudi::Generic3D< AIDA::IHistogram3D, TH3D >
std::numeric_limits