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