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