Loading [MathJax]/extensions/tex2jax.js
The Gaudi Framework  v31r0 (aeb156f0)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
H3D.cpp
Go to the documentation of this file.
1 #ifdef __ICC
2 // disable icc remark #2259: non-pointer conversion from "X" to "Y" may lose significant bits
3 // TODO: To be removed, since it comes from ROOT TMathBase.h
4 # pragma warning( disable : 2259 )
5 #endif
6 #ifdef WIN32
7 // Disable warning
8 // warning C4996: 'sprintf': This function or variable may be unsafe.
9 // coming from TString.h
10 # pragma warning( disable : 4996 )
11 #endif
12 
14 #include <GaudiCommonSvc/H3D.h>
16 #include <GaudiKernel/DataObject.h>
18 
19 #include "GaudiPI.h"
20 #include "TH3D.h"
21 
22 namespace Gaudi {
23  template <>
25  TH3D* imp = dynamic_cast<TH3D*>( rep );
26  if ( !imp ) throw std::runtime_error( "Cannot adopt native histogram representation." );
27  m_rep.reset( imp );
28  m_xAxis.initialize( m_rep->GetXaxis(), true );
29  m_yAxis.initialize( m_rep->GetYaxis(), true );
30  m_zAxis.initialize( m_rep->GetZaxis(), true );
31  const TArrayD* a = m_rep->GetSumw2();
32  if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
33  setTitle( m_rep->GetTitle() );
34  }
35 } // namespace Gaudi
36 
38 std::pair<DataObject*, AIDA::IHistogram3D*> Gaudi::createH3D( const std::string& title, int nBinsX, double xlow,
39  double xup, int nBinsY, double ylow, double yup,
40  int nBinsZ, double zlow, double zup ) {
41  auto p = new Histogram3D(
42  new TH3D( title.c_str(), title.c_str(), nBinsX, xlow, xup, nBinsY, ylow, yup, nBinsZ, zlow, zup ) );
43  return {p, p};
44 }
45 
48  const Edges& eY, const Edges& eZ ) {
49  auto p = new Histogram3D( new TH3D( title.c_str(), title.c_str(), eX.size() - 1, &eX.front(), eY.size() - 1,
50  &eY.front(), eZ.size() - 1, &eZ.front() ) );
51  return {p, p};
52 }
53 
54 std::pair<DataObject*, AIDA::IHistogram3D*> Gaudi::createH3D( const AIDA::IHistogram3D& hist ) {
55  TH3D* h = getRepresentation<AIDA::IHistogram3D, TH3D>( hist );
56  Histogram3D* n = h ? new Histogram3D( new TH3D( *h ) ) : nullptr;
57  return {n, n};
58 }
59 
61  setTitle( "" );
62  m_rep->Sumw2();
63  m_sumwx = 0;
64  m_sumwy = 0;
65  m_sumwz = 0;
66  m_rep->SetDirectory( nullptr );
67 }
68 
70  adoptRepresentation( rep );
71  m_sumwx = 0;
72  m_sumwy = 0;
73  m_sumwz = 0;
74  m_rep->SetDirectory( nullptr );
75 }
76 
77 // set bin content (entries and centre are not used )
78 bool Gaudi::Histogram3D::setBinContents( int i, int j, int k, int entries, double height, double error, double centreX,
79  double centreY, double centreZ ) {
80  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), rIndexZ( k ), height );
81  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), rIndexZ( k ), error );
82  // accumulate sum bin centers
83  if ( i >= 0 && j >= 0 && k >= 0 ) {
84  m_sumwx += centreX * height;
85  m_sumwy += centreY * height;
86  m_sumwz += centreZ * height;
87  }
89  return true;
90 }
91 
93  m_sumwx = 0;
94  m_sumwy = 0;
95  m_sumwz = 0;
96  m_sumEntries = 0;
97  m_rep->Reset();
98  return true;
99 }
100 
101 bool Gaudi::Histogram3D::fill( double x, double y, double z, double weight ) {
102  // avoid race conditiosn when filling the histogram
104  m_rep->Fill( x, y, z, weight );
105  return true;
106 }
107 
108 void* Gaudi::Histogram3D::cast( const std::string& className ) const {
109  if ( className == "AIDA::IHistogram3D" ) {
110  return (AIDA::IHistogram3D*)this;
111  } else if ( className == "AIDA::IHistogram" ) {
112  return (AIDA::IHistogram*)this;
113  }
114  return nullptr;
115 }
116 
117 #ifdef __ICC
118 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
119 // The comparison are meant
120 # pragma warning( push )
121 # pragma warning( disable : 1572 )
122 #endif
123 bool Gaudi::Histogram3D::setRms( double rmsX, double rmsY, double rmsZ ) {
124  m_rep->SetEntries( m_sumEntries );
125  std::vector<double> stat( 11 );
126  // sum weights
127  stat[0] = sumBinHeights();
128  stat[1] = 0;
129  if ( equivalentBinEntries() != 0 ) stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
130  stat[2] = m_sumwx;
131  double meanX = 0;
132  if ( sumBinHeights() != 0 ) meanX = m_sumwx / sumBinHeights();
133  stat[3] = ( meanX * meanX + rmsX * rmsX ) * sumBinHeights();
134  stat[4] = m_sumwy;
135  double meanY = 0;
136  if ( sumBinHeights() != 0 ) meanY = m_sumwy / sumBinHeights();
137  stat[5] = ( meanY * meanY + rmsY * rmsY ) * sumBinHeights();
138  stat[6] = 0;
139  stat[7] = m_sumwz;
140  double meanZ = 0;
141  if ( sumBinHeights() != 0 ) meanZ = m_sumwz / sumBinHeights();
142  stat[8] = ( meanZ * meanZ + rmsZ * rmsZ ) * sumBinHeights();
143  // do not need to use sumwxy sumwxz and sumwyz
144 
145  m_rep->PutStats( &stat.front() );
146  return true;
147 }
148 
149 void Gaudi::Histogram3D::copyFromAida( const AIDA::IHistogram3D& h ) {
150  // implement here the copy
151  const char* title = h.title().c_str();
152  if ( h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning() && h.zAxis().isFixedBinning() ) {
153  m_rep.reset( new TH3D( title, title, h.xAxis().bins(), h.xAxis().lowerEdge(), h.xAxis().upperEdge(),
154  h.yAxis().bins(), h.yAxis().lowerEdge(), h.yAxis().upperEdge(), h.zAxis().bins(),
155  h.zAxis().lowerEdge(), h.zAxis().upperEdge() ) );
156  } else {
157  Edges eX, eY, eZ;
158  for ( int i = 0; i < h.xAxis().bins(); ++i ) eX.push_back( h.xAxis().binLowerEdge( i ) );
159  // add also upperedges at the end
160  eX.push_back( h.xAxis().upperEdge() );
161  for ( int i = 0; i < h.yAxis().bins(); ++i ) eY.push_back( h.yAxis().binLowerEdge( i ) );
162  // add also upperedges at the end
163  eY.push_back( h.yAxis().upperEdge() );
164  for ( int i = 0; i < h.zAxis().bins(); ++i ) eZ.push_back( h.zAxis().binLowerEdge( i ) );
165  // add also upperedges at the end
166  eZ.push_back( h.zAxis().upperEdge() );
167  m_rep.reset(
168  new TH3D( title, title, eX.size() - 1, &eX.front(), eY.size() - 1, &eY.front(), eZ.size() - 1, &eZ.front() ) );
169  }
170  m_xAxis.initialize( m_rep->GetXaxis(), true );
171  m_yAxis.initialize( m_rep->GetYaxis(), true );
172  m_zAxis.initialize( m_rep->GetZaxis(), true );
173  const TArrayD* a = m_rep->GetSumw2();
174  if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
175  m_sumEntries = 0;
176  m_sumwx = 0;
177  m_sumwy = 0;
178  m_sumwz = 0;
179 
180  // statistics
181  double sumw = h.sumBinHeights();
182  double sumw2 = 0;
183  if ( h.equivalentBinEntries() != 0 ) sumw2 = ( sumw * sumw ) / h.equivalentBinEntries();
184  double sumwx = h.meanX() * h.sumBinHeights();
185  double sumwx2 = ( h.meanX() * h.meanX() + h.rmsX() * h.rmsX() ) * h.sumBinHeights();
186  double sumwy = h.meanY() * h.sumBinHeights();
187  double sumwy2 = ( h.meanY() * h.meanY() + h.rmsY() * h.rmsY() ) * h.sumBinHeights();
188  double sumwz = h.meanZ() * h.sumBinHeights();
189  double sumwz2 = ( h.meanZ() * h.meanZ() + h.rmsZ() * h.rmsZ() ) * h.sumBinHeights();
190  double sumwxy = 0;
191  double sumwxz = 0;
192  double sumwyz = 0;
193 
194  // copy the contents in (AIDA underflow/overflow are -2,-1)
195  for ( int i = -2; i < xAxis().bins(); ++i ) {
196  for ( int j = -2; j < yAxis().bins(); ++j ) {
197  for ( int k = -2; k < zAxis().bins(); ++k ) {
198  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), rIndexZ( k ), h.binHeight( i, j, k ) );
199  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), rIndexZ( k ), h.binError( i, j, k ) );
200  // calculate statistics
201  if ( i >= 0 && j >= 0 && k >= 0 ) {
202  sumwxy += h.binHeight( i, j, k ) * h.binMeanX( i, j, k ) * h.binMeanY( i, j, k );
203  sumwxz += h.binHeight( i, j, k ) * h.binMeanX( i, j, k ) * h.binMeanZ( i, j, k );
204  sumwyz += h.binHeight( i, j, k ) * h.binMeanY( i, j, k ) * h.binMeanZ( i, j, k );
205  }
206  }
207  }
208  }
209  // need to do set entries after setting contents otherwise root will recalulate them
210  // taking into account how many time SetBinContents() has been called
211  m_rep->SetEntries( h.allEntries() );
212 
213  // fill stat vector
214  std::vector<double> stat( 11 );
215  stat[0] = sumw;
216  stat[1] = sumw2;
217  stat[2] = sumwx;
218  stat[3] = sumwx2;
219  stat[4] = sumwy;
220  stat[5] = sumwy2;
221  stat[6] = sumwxy;
222  stat[7] = sumwz;
223  stat[8] = sumwz2;
224  stat[9] = sumwxz;
225  stat[10] = sumwyz;
226  m_rep->PutStats( &stat.front() );
227 }
228 
229 #ifdef __ICC
230 // re-enable icc remark #1572
231 # pragma warning( pop )
232 #endif
int rIndexZ(int index) const
Definition: Generic3D.h:83
int rIndexX(int index) const
Definition: Generic3D.h:81
Gaudi::Axis m_xAxis
Definition: Generic3D.h:213
double meanZ() const override
The mean of the IHistogram3D along the z axis.
Definition: Generic3D.h:162
Gaudi::Axis m_yAxis
Definition: Generic3D.h:214
int entries() const override
Get the number or all the entries.
Definition: Generic3D.h:241
T front(T...args)
std::mutex m_fillSerialization
Definition: H3D.h:48
const AIDA::IAxis & yAxis() const override
Get the y axis of the IHistogram3D.
Definition: Generic3D.h:172
std::vector< double > Edges
Definition: GaudiPI.h:17
const AIDA::IAxis & zAxis() const override
Get the z axis of the IHistogram3D.
Definition: Generic3D.h:174
double m_sumwx
Definition: H3D.h:43
double rmsZ() const override
The RMS of the IHistogram3D along the z axis.
Definition: Generic3D.h:168
STL class.
double m_sumwy
Definition: H3D.h:44
T push_back(T...args)
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:78
std::pair< DataObject *, AIDA::IHistogram3D * > createH3D(const AIDA::IHistogram3D &hist)
Copy constructor.
Histogram3D()
Standard Constructor.
Definition: H3D.cpp:60
AIDA implementation for 3 D histograms using ROOT THD2.
Definition: H3D.h:17
double sumBinHeights() const override
Get the sum of in range bin heights in the IProfile.
Definition: Generic3D.h:261
T reset(T...args)
Gaudi::Axis m_zAxis
Definition: Generic3D.h:215
bool setTitle(const std::string &title) override
Set the title of the object.
Definition: Generic3D.h:227
bool fill(double x, double y, double z, double weight) override
Fill bin content.
Definition: H3D.cpp:101
int rIndexY(int index) const
Definition: Generic3D.h:82
void copyFromAida(const AIDA::IHistogram3D &h)
Create new histogram from any AIDA based histogram.
Definition: H3D.cpp:149
std::unique_ptr< IMPLEMENTATION > m_rep
Reference to underlying implementation.
Definition: Generic3D.h:219
void initialize(TAxis *itaxi, bool)
Definition: Axis.h:63
T c_str(T...args)
const AIDA::IAxis & xAxis() const override
Get the x axis of the IHistogram3D.
Definition: Generic3D.h:170
virtual bool setRms(double rmsX, double rmsY, double rmsZ)
Sets the rms of the histogram.
Definition: H3D.cpp:123
double rmsY() const override
The RMS of the IHistogram3D along the y axis.
Definition: Generic3D.h:166
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
double rmsX() const override
The RMS of the IHistogram3D along the x axis.
Definition: Generic3D.h:164
double meanY() const override
The mean of the IHistogram3D along the y axis.
Definition: Generic3D.h:160
void * cast(const std::string &className) const override
Introspection method.
Definition: H3D.cpp:108
double equivalentBinEntries() const override
Number of equivalent entries, i.e. SUM[ weight ] ^ 2 / SUM[ weight^2 ]
Definition: Generic3D.h:271
bool reset() override
Definition: H3D.cpp:92
Helper functions to set/get the application return code.
Definition: __init__.py:1
double meanX() const override
The mean of the IHistogram3D along the x axis.
Definition: Generic3D.h:157
double m_sumwz
Definition: H3D.h:45
Common AIDA implementation stuff for histograms and profiles using ROOT implementations.
Definition: Generic3D.h:35