| The Gaudi Framework v27r1

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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 "GaudiKernel/DataObject.h"
14 #include "GaudiKernel/ObjectFactory.h"
15 #include "GaudiPI.h"
16 #include "Generic3D.h"
17 #include "TH3D.h"
18
19 namespace Gaudi {
20
27 class GAUDI_API Histogram3D : public DataObject, public Generic3D<AIDA::IHistogram3D,TH3D> {
28 public:
30 Histogram3D();
32 Histogram3D(TH3D* rep);
34 ~Histogram3D() override = default;
36 bool fill ( double x, double y, double z, double weight) override;
38 virtual bool setBinContents( int i, int j, int k, int entries,double height,double error,double centreX, double centreY, double centreZ );
40 virtual bool setRms(double rmsX, double rmsY, double rmsZ);
41 // overwrite reset
42 bool reset() override;
44 void* cast(const std::string & className) const override;
46 void copyFromAida(const AIDA::IHistogram3D & h);
48 const CLID& clID() const override { return classID(); }
49 static const CLID& classID() { return CLID_H3D; }
50
51 protected:
52 // cache sumwx and sumwy when setting contents since I don't have bin mean
53 double m_sumwx = 0;
54 double m_sumwy = 0;
55 double m_sumwz = 0;
56 };
57 }
58
59 namespace Gaudi {
60 template <>
61 void Generic3D<AIDA::IHistogram3D,TH3D>::adoptRepresentation(TObject* rep) {
62 TH3D* imp = dynamic_cast<TH3D*>(rep);
63 if ( !imp ) throw std::runtime_error("Cannot adopt native histogram representation.");
64 m_rep.reset( imp );
65 m_xAxis.initialize(m_rep->GetXaxis(),true);
66 m_yAxis.initialize(m_rep->GetYaxis(),true);
67 m_zAxis.initialize(m_rep->GetZaxis(),true);
68 const TArrayD* a = m_rep->GetSumw2();
69 if ( !a || (a && a->GetSize()==0) ) m_rep->Sumw2();
70 setTitle(m_rep->GetTitle());
71 }
72 }
73
75 std::pair<DataObject*,AIDA::IHistogram3D*>
76 Gaudi::createH3D(const std::string& title,int nBinsX,double xlow,double xup,int nBinsY,double ylow,double yup,int nBinsZ,double zlow,double zup) {
77 auto p = new Histogram3D(new TH3D(title.c_str(),title.c_str(),nBinsX,xlow,xup,nBinsY,ylow,yup,nBinsZ,zlow,zup));
78 return {p,p};
79 }
80
82 std::pair<DataObject*,AIDA::IHistogram3D*>
83 Gaudi::createH3D(const std::string& title, const Edges& eX,const Edges& eY,const Edges& eZ) {
84 auto p = new Histogram3D(new TH3D(title.c_str(),title.c_str(),eX.size()-1,&eX.front(), eY.size()-1,&eY.front(), eZ.size()-1,&eZ.front()));
85 return {p,p};
86 }
87
88 std::pair<DataObject*,AIDA::IHistogram3D*> Gaudi::createH3D(const AIDA::IHistogram3D& hist) {
89 TH3D *h = getRepresentation<AIDA::IHistogram3D,TH3D>(hist);
90 Histogram3D *n = h ? new Histogram3D(new TH3D(*h)) : nullptr;
91 return {n,n};
92 }
93
94 Gaudi::Histogram3D::Histogram3D()
95 : Base( new TH3D() )
96 {
97 setTitle("");
98 m_rep->Sumw2();
99 m_sumwx = 0;
100 m_sumwy = 0;
101 m_sumwz = 0;
102 m_rep->SetDirectory(nullptr);
103 }
104
105 Gaudi::Histogram3D::Histogram3D(TH3D* rep) {
106 adoptRepresentation(rep);
107 m_sumwx = 0;
108 m_sumwy = 0;
109 m_sumwz = 0;
110 m_rep->SetDirectory(nullptr);
111 }
112
113 // set bin content (entries and centre are not used )
114 bool Gaudi::Histogram3D::setBinContents(int i,int j,int k,int entries,double height,double error,double centreX, double centreY, double centreZ ) {
115 m_rep->SetBinContent(rIndexX(i), rIndexY(j), rIndexZ(k), height);
116 m_rep->SetBinError(rIndexX(i), rIndexY(j), rIndexZ(k), error);
117 // accumulate sum bin centers
118 if (i >=0 && j >= 0 && k >= 0) {
119 m_sumwx += centreX*height;
120 m_sumwy += centreY*height;
121 m_sumwz += centreZ*height;
122 }
123 m_sumEntries += entries;
124 return true;
125 }
126
127 bool Gaudi::Histogram3D::reset ( ) {
128 m_sumwx = 0;
129 m_sumwy = 0;
130 m_sumwz = 0;
131 m_sumEntries = 0;
132 m_rep->Reset ( );
133 return true;
134 }
135
136 bool Gaudi::Histogram3D::fill ( double x, double y, double z, double weight) {
137 m_rep->Fill ( x , y, z, weight );
138 return true;
139 }
140
141 void* Gaudi::Histogram3D::cast(const std::string & className) const {
142 if (className == "AIDA::IHistogram3D") {
143 return (AIDA::IHistogram3D*)this;
144 }
145 else if (className == "AIDA::IHistogram") {
146 return (AIDA::IHistogram*)this;
147 }
148 return nullptr;
149 }
150
151 #ifdef __ICC
152 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
153 // The comparison are meant
154 #pragma warning(push)
155 #pragma warning(disable:1572)
156 #endif
157 bool Gaudi::Histogram3D::setRms(double rmsX, double rmsY, double rmsZ ) {
158 m_rep->SetEntries(m_sumEntries);
159 std::vector<double> stat(11);
160 // sum weights
161 stat[0] = sumBinHeights();
162 stat[1] = 0;
163 if (equivalentBinEntries() != 0)
164 stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
165 stat[2] = m_sumwx;
166 double meanX = 0;
167 if ( sumBinHeights() != 0 ) meanX = m_sumwx/ sumBinHeights();
168 stat[3] = ( meanX*meanX + rmsX*rmsX )* sumBinHeights();
169 stat[4] = m_sumwy;
170 double meanY = 0;
171 if ( sumBinHeights() != 0 ) meanY = m_sumwy/ sumBinHeights();
172 stat[5] = ( meanY*meanY + rmsY*rmsY )* sumBinHeights();
173 stat[6] = 0;
174 stat[7] = m_sumwz;
175 double meanZ = 0;
176 if ( sumBinHeights() != 0 ) meanZ = m_sumwz/ sumBinHeights();
177 stat[8] = ( meanZ*meanZ + rmsZ*rmsZ )* sumBinHeights();
178 // do not need to use sumwxy sumwxz and sumwyz
179
180 m_rep->PutStats(&stat.front());
181 return true;
182 }
183
184 void Gaudi::Histogram3D::copyFromAida(const AIDA::IHistogram3D & h) {
185 // implement here the copy
186 const char* tit = h.title().c_str();
187 if (h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning() && h.zAxis().isFixedBinning() ) {
188 m_rep.reset( new TH3D(tit,tit,
189 h.xAxis().bins(), h.xAxis().lowerEdge(), h.xAxis().upperEdge(),
190 h.yAxis().bins(), h.yAxis().lowerEdge(), h.yAxis().upperEdge(),
191 h.zAxis().bins(), h.zAxis().lowerEdge(), h.zAxis().upperEdge() ) );
192 }
193 else {
194 Edges eX, eY, eZ;
195 for (int i =0; i < h.xAxis().bins(); ++i)
196 eX.push_back(h.xAxis().binLowerEdge(i));
197 // add also upperedges at the end
198 eX.push_back(h.xAxis().upperEdge() );
199 for (int i =0; i < h.yAxis().bins(); ++i)
200 eY.push_back(h.yAxis().binLowerEdge(i));
201 // add also upperedges at the end
202 eY.push_back(h.yAxis().upperEdge() );
203 for (int i =0; i < h.zAxis().bins(); ++i)
204 eZ.push_back(h.zAxis().binLowerEdge(i));
205 // add also upperedges at the end
206 eZ.push_back(h.zAxis().upperEdge() );
207 m_rep.reset( new TH3D(tit,tit,eX.size()-1,&eX.front(),eY.size()-1,&eY.front(),eZ.size()-1,&eZ.front()) );
208 }
209 m_xAxis.initialize(m_rep->GetXaxis(),true);
210 m_yAxis.initialize(m_rep->GetYaxis(),true);
211 m_zAxis.initialize(m_rep->GetZaxis(),true);
212 const TArrayD* a = m_rep->GetSumw2();
213 if ( !a || (a && a->GetSize()==0) ) m_rep->Sumw2();
214 m_sumEntries = 0;
215 m_sumwx = 0;
216 m_sumwy = 0;
217 m_sumwz = 0;
218
219 // statistics
220 double sumw = h.sumBinHeights();
221 double sumw2 = 0;
222 if (h.equivalentBinEntries() != 0)
223 sumw2 = ( sumw * sumw ) /h.equivalentBinEntries();
224 double sumwx = h.meanX()*h.sumBinHeights();
225 double sumwx2 = (h.meanX()*h.meanX() + h.rmsX()*h.rmsX() )*h.sumBinHeights();
226 double sumwy = h.meanY()*h.sumBinHeights();
227 double sumwy2 = (h.meanY()*h.meanY() + h.rmsY()*h.rmsY() )*h.sumBinHeights();
228 double sumwz = h.meanZ()*h.sumBinHeights();
229 double sumwz2 = (h.meanZ()*h.meanZ() + h.rmsZ()*h.rmsZ() )*h.sumBinHeights();
230 double sumwxy = 0;
231 double sumwxz = 0;
232 double sumwyz = 0;
233
234 // copy the contents in (AIDA underflow/overflow are -2,-1)
235 for (int i=-2; i < xAxis().bins(); ++i) {
236 for (int j=-2; j < yAxis().bins(); ++j) {
237 for (int k=-2; k < zAxis().bins(); ++k) {
238 m_rep->SetBinContent(rIndexX(i), rIndexY(j), rIndexZ(k), h.binHeight(i,j,k) );
239 m_rep->SetBinError(rIndexX(i), rIndexY(j), rIndexZ(k), h.binError(i,j,k) );
240 // calculate statistics
241 if ( i >= 0 && j >= 0 && k >= 0) {
242 sumwxy += h.binHeight(i,j,k)*h.binMeanX(i,j,k)*h.binMeanY(i,j,k);
243 sumwxz += h.binHeight(i,j,k)*h.binMeanX(i,j,k)*h.binMeanZ(i,j,k);
244 sumwyz += h.binHeight(i,j,k)*h.binMeanY(i,j,k)*h.binMeanZ(i,j,k);
245 }
246 }
247 }
248 }
249 // need to do set entries after setting contents otherwise root will recalulate them
250 // taking into account how many time SetBinContents() has been called
251 m_rep->SetEntries(h.allEntries());
252
253 // fill stat vector
254 std::vector<double> stat(11);
255 stat[0] = sumw;
256 stat[1] = sumw2;
257 stat[2] = sumwx;
258 stat[3] = sumwx2;
259 stat[4] = sumwy;
260 stat[5] = sumwy2;
261 stat[6] = sumwxy;
262 stat[7] = sumwz;
263 stat[8] = sumwz2;
264 stat[9] = sumwxz;
265 stat[10] = sumwyz;
266 m_rep->PutStats(&stat.front());
267 }
268
269 #ifdef __ICC
270 // re-enable icc remark #1572
271 #pragma warning(pop)
272 #endif
273
274 typedef Gaudi::Histogram3D H3D;
275 DECLARE_DATAOBJECT_FACTORY(H3D)
Gaudi::Histogram3D::classID
static const CLID & classID()
Definition: H3D.cpp:49
std::vector::front
T front(T...args)
ObjectFactory.h
Gaudi::Utils::Histos::fill
GAUDI_API void fill(AIDA::IHistogram1D *histo, const double value, const double weight=1.0)
simple function to fill AIDA::IHistogram1D objects
Definition: Fill.cpp:36
Gaudi::Edges
std::vector< double > Edges
Definition: GaudiPI.h:19
Gaudi::Histogram3D::m_sumwx
double m_sumwx
Definition: H3D.cpp:53
GaudiPluginService.cpluginsvc.n
list n
Definition: cpluginsvc.py:217
GaudiPI.h
Gaudi::Histogram3D::cast
void * cast(const std::string &className) const override
Introspection method.
Definition: H3D.cpp:141
std::string
STL class.
Gaudi::Histogram3D::m_sumwy
double m_sumwy
Definition: H3D.cpp:54
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:114
Gaudi::Histogram3D::Histogram3D
Histogram3D()
Standard Constructor.
Definition: H3D.cpp:94
Gaudi::Histogram3D
AIDA implementation for 2 D histograms using ROOT THD2.
Definition: H3D.cpp:27
std::runtime_error
STL class.
HiveAlgSequencer.h
tuple h
Definition: HiveAlgSequencer.py:24
Gaudi::Histogram3D::clID
const CLID & clID() const override
Retrieve reference to class defininition identifier.
Definition: H3D.cpp:48
Gaudi::Generic3D::setTitle
bool setTitle(const std::string &title) override
Set the title of the object.
Definition: Generic3D.h:228
CLID
unsigned int CLID
Class ID definition.
Definition: ClassID.h:8
Gaudi::Histogram3D::fill
bool fill(double x, double y, double z, double weight) override
Fill bin content.
Definition: H3D.cpp:136
Gaudi::Histogram3D::copyFromAida
void copyFromAida(const AIDA::IHistogram3D &h)
Create new histogram from any AIDA based histogram.
Definition: H3D.cpp:184
Gaudi::Generic3D::m_rep
std::unique_ptr< IMPLEMENTATION > m_rep
Reference to underlying implementation.
Definition: Generic3D.h:220
std::pair
std::vector< double >
Gaudi::createH3D
std::pair< DataObject *, AIDA::IHistogram3D * > createH3D(const AIDA::IHistogram3D &hist)
Copy constructor.
std::string::c_str
T c_str(T...args)
DECLARE_DATAOBJECT_FACTORY
#define DECLARE_DATAOBJECT_FACTORY(x)
Definition: ObjectFactory.h:21
Gaudi::Histogram3D::setRms
virtual bool setRms(double rmsX, double rmsY, double rmsZ)
Sets the rms of the histogram.
Definition: H3D.cpp:157
Gaudi::Generic3D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
Generic3D.h
H3D
Gaudi::Histogram3D H3D
Definition: H3D.cpp:274
GAUDI_API
#define GAUDI_API
Definition: Kernel.h:107
DataObject
A DataObject is the base class of any identifiable object on any data store.
Definition: DataObject.h:30
ana.i
list i
Definition: ana.py:128
Gaudi::Histogram3D::reset
bool reset() override
Definition: H3D.cpp:127
Gaudi
Helper functions to set/get the application return code.
Definition: __init__.py:1
Gaudi::Histogram3D::m_sumwz
double m_sumwz
Definition: H3D.cpp:55
DataObject.h
Gaudi::Generic3D
Common AIDA implementation stuff for histograms and profiles using ROOT implementations.
Definition: Generic3D.h:34