| The Gaudi Framework v32r1 (f65d50dc)

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
13 #include <GaudiCommonSvc/Generic3D.h>
14 #include <GaudiCommonSvc/H3D.h>
15 #include <GaudiCommonSvc/HistogramUtility.h>
16 #include <GaudiKernel/DataObject.h>
17 #include <GaudiKernel/ObjectFactory.h>
18
19 #include "GaudiPI.h"
20 #include "TH3D.h"
21
22 namespace Gaudi {
23 template <>
24 void Generic3D<AIDA::IHistogram3D, TH3D>::adoptRepresentation( TObject* rep ) {
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
47 std::pair<DataObject*, AIDA::IHistogram3D*> Gaudi::createH3D( const std::string& title, const Edges& eX,
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
60 Gaudi::Histogram3D::Histogram3D() : Base( new TH3D() ) {
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
69 Gaudi::Histogram3D::Histogram3D( TH3D* rep ) {
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 }
88 m_sumEntries += entries;
89 return true;
90 }
91
92 bool Gaudi::Histogram3D::reset() {
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
103 std::lock_guard<std::mutex> guard( m_fillSerialization );
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
std::vector::front
T front(T... args)
ObjectFactory.h
H3D.h
Gaudi::Edges
std::vector< double > Edges
Definition: GaudiPI.h:17
Gaudi::Histogram3D::m_sumwx
double m_sumwx
Definition: H3D.h:43
GaudiPI.h
Gaudi::Histogram3D::cast
void * cast(const std::string &className) const override
Introspection method.
Definition: H3D.cpp:108
std::string
STL class.
GaudiPluginService.cpluginsvc.n
n
Definition: cpluginsvc.py:210
Gaudi::Histogram3D::m_sumwy
double m_sumwy
Definition: H3D.h:44
std::vector::push_back
T push_back(T... args)
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:78
Gaudi::createH3D
std::pair< DataObject *, AIDA::IHistogram3D * > createH3D(const AIDA::IHistogram3D &hist)
Copy constructor.
Gaudi::Histogram3D::Histogram3D
Histogram3D()
Standard Constructor.
Definition: H3D.cpp:60
std::runtime_error
STL class.
Gaudi::Generic3D::setTitle
bool setTitle(const std::string &title) override
Set the title of the object.
Definition: Generic3D.h:227
Gaudi::Histogram3D::fill
bool fill(double x, double y, double z, double weight) override
Fill bin content.
Definition: H3D.cpp:101
Gaudi::Histogram3D::copyFromAida
void copyFromAida(const AIDA::IHistogram3D &h)
Create new histogram from any AIDA based histogram.
Definition: H3D.cpp:149
Gaudi::Generic3D::m_rep
std::unique_ptr< IMPLEMENTATION > m_rep
Reference to underlying implementation.
Definition: Generic3D.h:219
std::pair
std::vector< double >
std::lock_guard
std::string::c_str
T c_str(T... args)
Gaudi::Histogram3D::setRms
virtual bool setRms(double rmsX, double rmsY, double rmsZ)
Sets the rms of the histogram.
Definition: H3D.cpp:123
Gaudi::Generic3D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
Generic3D.h
HistogramUtility.h
Gaudi::Histogram3D::reset
bool reset() override
Definition: H3D.cpp:92
Gaudi
Header file for std:chrono::duration-based Counters.
Definition: __init__.py:1
Gaudi::Histogram3D::m_sumwz
double m_sumwz
Definition: H3D.h:45
DataObject.h
Gaudi::Generic3D
Common AIDA implementation stuff for histograms and profiles using ROOT implementations.
Definition: Generic3D.h:35