| # 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
H2D.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 #include "H2D.h"
13 #include "H1D.h"
14 #include "P1D.h"
15 #include "TH1D.h"
16 #include "TProfile.h"
17 #include "GaudiKernel/ObjectFactory.h"
18
19 std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(const std::string & title,int binsX,double iminX,double imaxX,int binsY,double iminY,double imaxY) {
20 Histogram2D* p = new Histogram2D(new TH2D(title.c_str(),title.c_str(),binsX, iminX, imaxX, binsY, iminY, imaxY));
21 return std::pair<DataObject*,IHistogram2D*>(p,p);
22 }
23
24 std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(const std::string & title, const Edges& eX, const Edges& eY) {
25 Histogram2D* p = new Histogram2D(new TH2D(title.c_str(),title.c_str(),eX.size()-1,&eX.front(),eY.size()-1,&eY.front()));
26 return std::pair<DataObject*,IHistogram2D*>(p,p);
27 }
28
29 std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(TH2D* rep) {
30 Histogram2D* p = new Histogram2D(rep);
31 return std::pair<DataObject*,IHistogram2D*>(p,p);
32 }
33
34 std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(const IHistogram2D& hist) {
35 TH2D *h = getRepresentation<AIDA::IHistogram2D,TH2D>(hist);
36 Histogram2D *n = h ? new Histogram2D(new TH2D(*h)) : 0;
37 return std::pair<DataObject*,IHistogram2D*>(n,n);
38 }
39
40 std::pair<DataObject*,IHistogram1D*>
41 Gaudi::slice1DX(const std::string& nam,const IHistogram2D& hist,int first, int last) {
42 TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
43 TH1D *t = r ? r->ProjectionX("_px",first,last,"e") : 0;
44 if ( t ) t->SetName(nam.c_str());
45 Histogram1D* p = t ? new Histogram1D(t) : 0;
46 return std::pair<DataObject*,IHistogram1D*>(p,p);
47 }
48
49 std::pair<DataObject*,IHistogram1D*>
50 Gaudi::slice1DY(const std::string& nam,const IHistogram2D& hist,int first, int last) {
51 TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
52 TH1D *t = r ? r->ProjectionY("_py",first,last,"e") : 0;
53 if ( t ) t->SetName(nam.c_str());
54 Histogram1D* p = t ? new Histogram1D(t) : 0;
55 return std::pair<DataObject*,IHistogram1D*>(p,p);
56 }
57
58 std::pair<DataObject*,IHistogram1D*>
59 Gaudi::project1DY(const std::string& nam,const IHistogram2D& hist,int first,int last) {
60 TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
61 TH1D *t = r ? r->ProjectionY("_px",first,last,"e") : 0;
62 if ( t ) t->SetName(nam.c_str());
63 Histogram1D* p = t ? new Histogram1D(t) : 0;
64 return std::pair<DataObject*,IHistogram1D*>(p,p);
65 }
66
67 std::pair<DataObject*,IProfile1D*>
68 Gaudi::profile1DX(const std::string& nam,const IHistogram2D& hist,int first,int last) {
69 TH2 *r = Gaudi::getRepresentation<IHistogram2D,TH2>(hist);
70 TProfile *t = r ? r->ProfileX("_pfx",first,last,"e") : 0;
71 if ( t ) t->SetName(nam.c_str());
72 Profile1D* p = t ? new Profile1D(t) : 0;
73 return std::pair<DataObject*,IProfile1D*>(p,p);
74 }
75
76 std::pair<DataObject*,IProfile1D*>
77 Gaudi::profile1DY(const std::string& nam, const IHistogram2D& hist,int first, int last) {
78 TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
79 TProfile *t = r ? r->ProfileY("_pfx",first,last,"e") : 0;
80 if ( t ) t->SetName(nam.c_str());
81 Profile1D* p = t ? new Profile1D(t) : 0;
82 return std::pair<DataObject*,IProfile1D*>(p,p);
83 }
84
85
86 namespace Gaudi {
87 template <>
88 void * Generic2D<IHistogram2D,TH2D>::cast(const std::string & className) const {
89 if (className == "AIDA::IHistogram2D")
90 return (IHistogram2D*)this;
91 else if (className == "AIDA::IHistogram")
92 return (IHistogram*)this;
93 return 0;
94 }
95
96 template <>
97 int Generic2D<IHistogram2D,TH2D>::binEntries(int indexX,int indexY) const {
98 if (binHeight(indexX, indexY)<=0) return 0;
99 double xx = binHeight(indexX, indexY)/binError(indexX, indexY);
100 return int(xx*xx+0.5);
101 }
102
103 template <>
104 void Generic2D<AIDA::IHistogram2D,TH2D>::adoptRepresentation(TObject* rep) {
105 TH2D* imp = dynamic_cast<TH2D*>(rep);
106 if ( imp ) {
107 if ( m_rep ) delete m_rep;
108 m_rep = imp;
109 m_xAxis.initialize(m_rep->GetXaxis(),true);
110 m_yAxis.initialize(m_rep->GetYaxis(),true);
111 const TArrayD* a = m_rep->GetSumw2();
112 if ( 0 == a || (a && a->GetSize()==0) ) m_rep->Sumw2();
113 setTitle(m_rep->GetTitle());
114 return;
115 }
116 throw std::runtime_error("Cannot adopt native histogram representation.");
117 }
118 }
119
120 Gaudi::Histogram2D::Histogram2D() {
121 m_rep = new TH2D();
122 m_rep->Sumw2();
123 m_sumEntries = 0;
124 m_sumwx = m_sumwy = 0;
125 setTitle("");
126 m_rep->SetDirectory(0);
127 }
128
129 Gaudi::Histogram2D::Histogram2D(TH2D* rep) {
130 m_rep = 0;
131 adoptRepresentation(rep);
132 m_sumEntries = 0;
133 m_sumwx = m_sumwy = 0;
134 m_rep->SetDirectory(0);
135 }
136
137 bool Gaudi::Histogram2D::setBinContents( int i,int j,int entries,double height,double error,double centreX,double centreY) {
138 m_rep->SetBinContent(rIndexX(i), rIndexY(j), height);
139 m_rep->SetBinError(rIndexX(i), rIndexY(j), error);
140 // accumulate sumwx for in range bins
141 if (i >=0 && j >= 0) {
142 m_sumwx += centreX*height;
143 m_sumwy += centreY*height;
144 }
145 m_sumEntries += entries;
146 return true;
147 }
148
149 bool Gaudi::Histogram2D::reset() {
150 m_sumwx = 0;
151 m_sumwy = 0;
152 return Base::reset();
153 }
154
155 #ifdef __ICC
156 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
157 // The comparison are meant
158 #pragma warning(push)
159 #pragma warning(disable:1572)
160 #endif
161 bool Gaudi::Histogram2D::fill ( double x,double y,double weight) {
162 (weight==1.) ? m_rep->Fill(x,y) : m_rep->Fill(x,y,weight );
163 return true;
164 }
165
166 bool Gaudi::Histogram2D::setRms(double rmsX,double rmsY) {
167 m_rep->SetEntries(m_sumEntries);
168 std::vector<double> stat(11);
169 stat[0] = sumBinHeights();
170 stat[1] = 0;
171 if(equivalentBinEntries() != 0)
172 stat[1] = (sumBinHeights() * sumBinHeights()) / equivalentBinEntries();
173 stat[2] = m_sumwx;
174 double meanX = 0;
175 if(sumBinHeights() != 0) meanX = m_sumwx/ sumBinHeights();
176 stat[3] = (meanX*meanX + rmsX*rmsX) * sumBinHeights();
177 stat[4] = m_sumwy;
178 double meanY = 0;
179 if(sumBinHeights() != 0) meanY = m_sumwy/ sumBinHeights();
180 stat[5] = (meanY*meanY + rmsY*rmsY) * sumBinHeights();
181 stat[6] = 0;
182 m_rep->PutStats(&stat.front());
183 return true;
184 }
185
186 void Gaudi::Histogram2D::copyFromAida(const IHistogram2D& h) {
187 // implement here the copy
188 delete m_rep;
189 const char* tit = h.title().c_str();
190 if (h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning() )
191 m_rep = new TH2D(tit,tit,
192 h.xAxis().bins(),h.xAxis().lowerEdge(),h.xAxis().upperEdge(),
193 h.yAxis().bins(),h.yAxis().lowerEdge(),h.yAxis().upperEdge() );
194 else {
195 Edges eX, eY;
196 for (int i =0; i < h.xAxis().bins(); ++i)
197 eX.push_back(h.xAxis().binLowerEdge(i));
198 // add also upperedges at the end
199 eX.push_back(h.xAxis().upperEdge() );
200 for (int i =0; i < h.yAxis().bins(); ++i)
201 eY.push_back(h.yAxis().binLowerEdge(i));
202 // add also upperedges at the end
203 eY.push_back(h.yAxis().upperEdge() );
204 m_rep = new TH2D(tit,tit,eX.size()-1,&eX.front(),eY.size()-1,&eY.front());
205 }
206 m_xAxis.initialize(m_rep->GetXaxis(),true);
207 m_yAxis.initialize(m_rep->GetYaxis(),true);
208 m_rep->Sumw2();
209 m_sumEntries = 0;
210 m_sumwx = 0;
211 m_sumwy = 0;
212 // statistics
213 double sumw = h.sumBinHeights();
214 double sumw2 = 0;
215 if (h.equivalentBinEntries() != 0)
216 sumw2 = ( sumw * sumw ) /h.equivalentBinEntries();
217 double sumwx = h.meanX()*h.sumBinHeights();
218 double sumwx2 = (h.meanX()*h.meanX() + h.rmsX()*h.rmsX() )*h.sumBinHeights();
219 double sumwy = h.meanY()*h.sumBinHeights();
220 double sumwy2 = (h.meanY()*h.meanY() + h.rmsY()*h.rmsY() )*h.sumBinHeights();
221 double sumwxy = 0;
222
223 // copy the contents in (AIDA underflow/overflow are -2,-1)
224 for (int i=-2; i < xAxis().bins(); ++i) {
225 for (int j=-2; j < yAxis().bins(); ++j) {
226 // root binning starts from one !
227 m_rep->SetBinContent(rIndexX(i), rIndexY(j), h.binHeight(i,j) );
228 m_rep->SetBinError(rIndexX(i), rIndexY(j), h.binError(i,j) );
229 // calculate statistics
230 if ( i >= 0 && j >= 0) {
231 sumwxy += h.binHeight(i,j)*h.binMeanX(i,j)*h.binMeanY(i,j);
232 }
233 }
234 }
235 // need to do set entries after setting contents otherwise root will recalculate them
236 // taking into account how many time SetBinContents() has been called
237 m_rep->SetEntries(h.allEntries());
238 // fill stat vector
239 std::vector<double> stat(11);
240 stat[0] = sumw;
241 stat[1] = sumw2;
242 stat[2] = sumwx;
243 stat[3] = sumwx2;
244 stat[4] = sumwy;
245 stat[5] = sumwy2;
246 stat[6] = sumwxy;
247 m_rep->PutStats(&stat.front());
248 }
249
250 typedef Gaudi::Histogram2D H2D;
251 DECLARE_DATAOBJECT_FACTORY(H2D)
252
253 #ifdef __ICC
254 // re-enable icc remark #1572
255 #pragma warning(pop)
256 #endif