| The Gaudi Framework v32r1 (f65d50dc)

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