| # Gaudi Framework, version v25r0 |

Home |  | Generated: Mon Feb 17 2014

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

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups 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
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 virtual ~Histogram3D() {}
36 bool fill ( double x, double y, double z, double weight);
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 virtual bool reset();
44 void* cast(const std::string & className) const;
46 void copyFromAida(const AIDA::IHistogram3D & h);
48 virtual const CLID& clID() const { 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;
54 double m_sumwy;
55 double m_sumwz;
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 ) {
64 if ( m_rep ) delete m_rep;
65 m_rep = imp;
66 m_xAxis.initialize(m_rep->GetXaxis(),true);
67 m_yAxis.initialize(m_rep->GetYaxis(),true);
68 m_zAxis.initialize(m_rep->GetZaxis(),true);
69 const TArrayD* a = m_rep->GetSumw2();
70 if ( 0 == a || (a && a->GetSize()==0) ) m_rep->Sumw2();
71 setTitle(m_rep->GetTitle());
72 return;
73 }
74 throw std::runtime_error("Cannot adopt native histogram representation.");
75 }
76 }
77
79 std::pair<DataObject*,AIDA::IHistogram3D*>
80 Gaudi::createH3D(const std::string& title,int nBinsX,double xlow,double xup,int nBinsY,double ylow,double yup,int nBinsZ,double zlow,double zup) {
81 Histogram3D* p = new Histogram3D(new TH3D(title.c_str(),title.c_str(),nBinsX,xlow,xup,nBinsY,ylow,yup,nBinsZ,zlow,zup));
82 return std::pair<DataObject*,AIDA::IHistogram3D*>(p,p);
83 }
84
86 std::pair<DataObject*,AIDA::IHistogram3D*>
87 Gaudi::createH3D(const std::string& title, const Edges& eX,const Edges& eY,const Edges& eZ) {
88 Histogram3D* 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()));
89 return std::pair<DataObject*,AIDA::IHistogram3D*>(p,p);
90 }
91
92 std::pair<DataObject*,AIDA::IHistogram3D*> Gaudi::createH3D(const AIDA::IHistogram3D& hist) {
93 TH3D *h = getRepresentation<AIDA::IHistogram3D,TH3D>(hist);
94 Histogram3D *n = h ? new Histogram3D(new TH3D(*h)) : 0;
95 return std::pair<DataObject*,AIDA::IHistogram3D*>(n,n);
96 }
97
98 Gaudi::Histogram3D::Histogram3D() {
99 m_rep = new TH3D();
100 setTitle("");
101 m_rep->Sumw2();
102 m_sumEntries = 0;
103 m_sumwx = 0;
104 m_sumwy = 0;
105 m_sumwz = 0;
106 m_rep->SetDirectory(0);
107 }
108
109 Gaudi::Histogram3D::Histogram3D(TH3D* rep) {
110 m_rep = 0;
111 adoptRepresentation(rep);
112 m_sumEntries = 0;
113 m_sumwx = 0;
114 m_sumwy = 0;
115 m_sumwz = 0;
116 m_rep->SetDirectory(0);
117 }
118
119 // set bin content (entries and centre are not used )
120 bool Gaudi::Histogram3D::setBinContents(int i,int j,int k,int entries,double height,double error,double centreX, double centreY, double centreZ ) {
121 m_rep->SetBinContent(rIndexX(i), rIndexY(j), rIndexZ(k), height);
122 m_rep->SetBinError(rIndexX(i), rIndexY(j), rIndexZ(k), error);
123 // accumulate sum bin centers
124 if (i >=0 && j >= 0 && k >= 0) {
125 m_sumwx += centreX*height;
126 m_sumwy += centreY*height;
127 m_sumwz += centreZ*height;
128 }
129 m_sumEntries += entries;
130 return true;
131 }
132
133 bool Gaudi::Histogram3D::reset ( ) {
134 m_sumwx = 0;
135 m_sumwy = 0;
136 m_sumwz = 0;
137 m_sumEntries = 0;
138 m_rep->Reset ( );
139 return true;
140 }
141
142 bool Gaudi::Histogram3D::fill ( double x, double y, double z, double weight) {
143 m_rep->Fill ( x , y, z, weight );
144 return true;
145 }
146
147 void* Gaudi::Histogram3D::cast(const std::string & className) const {
148 if (className == "AIDA::IHistogram3D") {
149 return (AIDA::IHistogram3D*)this;
150 }
151 else if (className == "AIDA::IHistogram") {
152 return (AIDA::IHistogram*)this;
153 }
154 return 0;
155 }
156
157 #ifdef __ICC
158 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
159 // The comparison are meant
160 #pragma warning(push)
161 #pragma warning(disable:1572)
162 #endif
163 bool Gaudi::Histogram3D::setRms(double rmsX, double rmsY, double rmsZ ) {
164 m_rep->SetEntries(m_sumEntries);
165 std::vector<double> stat(11);
166 // sum weights
167 stat[0] = sumBinHeights();
168 stat[1] = 0;
169 if (equivalentBinEntries() != 0)
170 stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
171 stat[2] = m_sumwx;
172 double meanX = 0;
173 if ( sumBinHeights() != 0 ) meanX = m_sumwx/ sumBinHeights();
174 stat[3] = ( meanX*meanX + rmsX*rmsX )* sumBinHeights();
175 stat[4] = m_sumwy;
176 double meanY = 0;
177 if ( sumBinHeights() != 0 ) meanY = m_sumwy/ sumBinHeights();
178 stat[5] = ( meanY*meanY + rmsY*rmsY )* sumBinHeights();
179 stat[6] = 0;
180 stat[7] = m_sumwz;
181 double meanZ = 0;
182 if ( sumBinHeights() != 0 ) meanZ = m_sumwz/ sumBinHeights();
183 stat[8] = ( meanZ*meanZ + rmsZ*rmsZ )* sumBinHeights();
184 // do not need to use sumwxy sumwxz and sumwyz
185
186 m_rep->PutStats(&stat.front());
187 return true;
188 }
189
190 void Gaudi::Histogram3D::copyFromAida(const AIDA::IHistogram3D & h) {
191 delete m_rep;
192 // implement here the copy
193 const char* tit = h.title().c_str();
194 if (h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning() && h.zAxis().isFixedBinning() ) {
195 m_rep = new TH3D(tit,tit,
196 h.xAxis().bins(), h.xAxis().lowerEdge(), h.xAxis().upperEdge(),
197 h.yAxis().bins(), h.yAxis().lowerEdge(), h.yAxis().upperEdge(),
198 h.zAxis().bins(), h.zAxis().lowerEdge(), h.zAxis().upperEdge() );
199 }
200 else {
201 Edges eX, eY, eZ;
202 for (int i =0; i < h.xAxis().bins(); ++i)
203 eX.push_back(h.xAxis().binLowerEdge(i));
204 // add also upperedges at the end
205 eX.push_back(h.xAxis().upperEdge() );
206 for (int i =0; i < h.yAxis().bins(); ++i)
207 eY.push_back(h.yAxis().binLowerEdge(i));
208 // add also upperedges at the end
209 eY.push_back(h.yAxis().upperEdge() );
210 for (int i =0; i < h.zAxis().bins(); ++i)
211 eZ.push_back(h.zAxis().binLowerEdge(i));
212 // add also upperedges at the end
213 eZ.push_back(h.zAxis().upperEdge() );
214 m_rep = new TH3D(tit,tit,eX.size()-1,&eX.front(),eY.size()-1,&eY.front(),eZ.size()-1,&eZ.front());
215 }
216 m_xAxis.initialize(m_rep->GetXaxis(),true);
217 m_yAxis.initialize(m_rep->GetYaxis(),true);
218 m_zAxis.initialize(m_rep->GetZaxis(),true);
219 const TArrayD* a = m_rep->GetSumw2();
220 if ( 0 == a || (a && a->GetSize()==0) ) m_rep->Sumw2();
221 m_sumEntries = 0;
222 m_sumwx = 0;
223 m_sumwy = 0;
224 m_sumwz = 0;
225
226 // statistics
227 double sumw = h.sumBinHeights();
228 double sumw2 = 0;
229 if (h.equivalentBinEntries() != 0)
230 sumw2 = ( sumw * sumw ) /h.equivalentBinEntries();
231 double sumwx = h.meanX()*h.sumBinHeights();
232 double sumwx2 = (h.meanX()*h.meanX() + h.rmsX()*h.rmsX() )*h.sumBinHeights();
233 double sumwy = h.meanY()*h.sumBinHeights();
234 double sumwy2 = (h.meanY()*h.meanY() + h.rmsY()*h.rmsY() )*h.sumBinHeights();
235 double sumwz = h.meanZ()*h.sumBinHeights();
236 double sumwz2 = (h.meanZ()*h.meanZ() + h.rmsZ()*h.rmsZ() )*h.sumBinHeights();
237 double sumwxy = 0;
238 double sumwxz = 0;
239 double sumwyz = 0;
240
241 // copy the contents in (AIDA underflow/overflow are -2,-1)
242 for (int i=-2; i < xAxis().bins(); ++i) {
243 for (int j=-2; j < yAxis().bins(); ++j) {
244 for (int k=-2; k < zAxis().bins(); ++k) {
245 m_rep->SetBinContent(rIndexX(i), rIndexY(j), rIndexZ(k), h.binHeight(i,j,k) );
246 m_rep->SetBinError(rIndexX(i), rIndexY(j), rIndexZ(k), h.binError(i,j,k) );
247 // calculate statistics
248 if ( i >= 0 && j >= 0 && k >= 0) {
249 sumwxy += h.binHeight(i,j,k)*h.binMeanX(i,j,k)*h.binMeanY(i,j,k);
250 sumwxz += h.binHeight(i,j,k)*h.binMeanX(i,j,k)*h.binMeanZ(i,j,k);
251 sumwyz += h.binHeight(i,j,k)*h.binMeanY(i,j,k)*h.binMeanZ(i,j,k);
252 }
253 }
254 }
255 }
256 // need to do set entries after setting contents otherwise root will recalulate them
257 // taking into account how many time SetBinContents() has been called
258 m_rep->SetEntries(h.allEntries());
259
260 // fill stat vector
261 std::vector<double> stat(11);
262 stat[0] = sumw;
263 stat[1] = sumw2;
264 stat[2] = sumwx;
265 stat[3] = sumwx2;
266 stat[4] = sumwy;
267 stat[5] = sumwy2;
268 stat[6] = sumwxy;
269 stat[7] = sumwz;
270 stat[8] = sumwz2;
271 stat[9] = sumwxz;
272 stat[10] = sumwyz;
273 m_rep->PutStats(&stat.front());
274 }
275
276 #ifdef __ICC
277 // re-enable icc remark #1572
278 #pragma warning(pop)
279 #endif
280
281 typedef Gaudi::Histogram3D H3D;
282 DECLARE_DATAOBJECT_FACTORY(H3D)