| The Gaudi Framework v32r1 (f65d50dc)

H1D.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/HistogramUtility.h>
15 #include <GaudiKernel/ObjectFactory.h>
16 #include <GaudiKernel/StreamBuffer.h>
17
18 std::pair<DataObject*, AIDA::IHistogram1D*> Gaudi::createH1D( const std::string& title, int nBins, double xlow,
19 double xup ) {
20 auto p = new Histogram1D( new TH1D( title.c_str(), title.c_str(), nBins, xlow, xup ) );
21 return {p, p};
22 }
23
24 std::pair<DataObject*, AIDA::IHistogram1D*> Gaudi::createH1D( const std::string& title, const Edges& e ) {
25 auto p = new Histogram1D( new TH1D( title.c_str(), title.c_str(), e.size() - 1, &e.front() ) );
26 return {p, p};
27 }
28
29 std::pair<DataObject*, AIDA::IHistogram1D*> Gaudi::createH1D( const AIDA::IHistogram1D& hist ) {
30 TH1D* h = getRepresentation<AIDA::IHistogram1D, TH1D>( hist );
31 auto n = ( h ? new Histogram1D( new TH1D( *h ) ) : nullptr );
32 return {n, n};
33 }
34 namespace Gaudi {
35
36 template <>
37 void* Generic1D<AIDA::IHistogram1D, TH1D>::cast( const std::string& className ) const {
38 if ( className == "AIDA::IHistogram1D" ) return const_cast<AIDA::IHistogram1D*>( (AIDA::IHistogram1D*)this );
39 if ( className == "AIDA::IHistogram" ) return const_cast<AIDA::IHistogram*>( (AIDA::IHistogram*)this );
40 return nullptr;
41 }
42
43 template <>
44 int Generic1D<AIDA::IHistogram1D, TH1D>::binEntries( int index ) const {
45 if ( binHeight( index ) <= 0 ) return 0;
46 double xx = binHeight( index ) / binError( index );
47 return int( xx * xx + 0.5 );
48 }
49
50 template <>
51 void Generic1D<AIDA::IHistogram1D, TH1D>::adoptRepresentation( TObject* rep ) {
52 TH1D* imp = dynamic_cast<TH1D*>( rep );
53 if ( !imp ) throw std::runtime_error( "Cannot adopt native histogram representation." );
54 m_rep.reset( imp );
55 }
56 } // namespace Gaudi
57
58 Gaudi::Histogram1D::Histogram1D() : Base( new TH1D() ) { init( "", false ); }
59
60 Gaudi::Histogram1D::Histogram1D( TH1D* rep ) {
61 m_rep.reset( rep );
62 init( m_rep->GetTitle() );
63 initSums();
64 }
65
66 void Gaudi::Histogram1D::init( const std::string& title, bool initialize_axis ) {
67 m_classType = "IHistogram1D";
68 if ( initialize_axis ) { m_axis.initialize( m_rep->GetXaxis(), false ); }
69 const TArrayD* a = m_rep->GetSumw2();
70 if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
71 setTitle( title );
72 m_rep->SetDirectory( nullptr );
73 m_sumEntries = 0;
74 m_sumwx = 0;
75 }
76
77 void Gaudi::Histogram1D::initSums() {
78 m_sumwx = 0;
79 m_sumEntries = 0;
80 for ( int i = 1, n = m_rep->GetNbinsX(); i <= n; ++i ) {
81 m_sumwx += m_rep->GetBinContent( i ) * m_rep->GetBinCenter( i );
82 m_sumEntries += m_rep->GetBinContent( i );
83 }
84 }
85
86 bool Gaudi::Histogram1D::reset() {
87 m_sumwx = 0;
88 m_sumEntries = 0;
89 return Base::reset();
90 }
91
93 void Gaudi::Histogram1D::adoptRepresentation( TObject* rep ) {
94 Gaudi::Generic1D<AIDA::IHistogram1D, TH1D>::adoptRepresentation( rep );
95 if ( m_rep ) {
96 init( m_rep->GetTitle() );
97 initSums();
98 }
99 }
100
101 bool Gaudi::Histogram1D::setBinContents( int i, int entries, double height, double error, double centre ) {
102 m_rep->SetBinContent( rIndex( i ), height );
103 m_rep->SetBinError( rIndex( i ), error );
104 // accumulate sumwx for in range bins
105 if ( i != AIDA::IAxis::UNDERFLOW_BIN && i != AIDA::IAxis::OVERFLOW_BIN ) m_sumwx += centre * height;
106 m_sumEntries += entries;
107 return true;
108 }
109
110 #ifdef __ICC
111 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
112 // The comparison are meant
113 # pragma warning( push )
114 # pragma warning( disable : 1572 )
115 #endif
116 bool Gaudi::Histogram1D::setRms( double rms ) {
117 m_rep->SetEntries( m_sumEntries );
118 std::vector<double> stat( 11 );
119 // sum weights
120 stat[0] = sumBinHeights();
121 stat[1] = 0;
122 if ( equivalentBinEntries() != 0 ) stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
123 stat[2] = m_sumwx;
124 double mean = 0;
125 if ( sumBinHeights() != 0 ) mean = m_sumwx / sumBinHeights();
126 stat[3] = ( mean * mean + rms * rms ) * sumBinHeights();
127 m_rep->PutStats( &stat.front() );
128 return true;
129 }
130
131 // set histogram statistics
132 bool Gaudi::Histogram1D::setStatistics( int allEntries, double eqBinEntries, double mean, double rms ) {
133 m_rep->SetEntries( allEntries );
134 // fill statistcal vector for Root
135 std::vector<double> stat( 11 );
136 // sum weights
137 stat[0] = sumBinHeights();
138 // sum weights **2
139 stat[1] = 0;
140 if ( eqBinEntries != 0 ) stat[1] = ( sumBinHeights() * sumBinHeights() ) / eqBinEntries;
141 // sum weights * x
142 stat[2] = mean * sumBinHeights();
143 // sum weight * x **2
144 stat[3] = ( mean * mean + rms * rms ) * sumBinHeights();
145 m_rep->PutStats( &stat.front() );
146 return true;
147 }
148
149 bool Gaudi::Histogram1D::fill( double x, double weight ) {
150 // avoid race conditions when filling the histogram
151 std::lock_guard<std::mutex> guard( m_fillSerialization );
152 ( weight == 1. ) ? m_rep->Fill( x ) : m_rep->Fill( x, weight );
153 return true;
154 }
155
156 void Gaudi::Histogram1D::copyFromAida( const AIDA::IHistogram1D& h ) {
157 // implement here the copy
158 std::string title = h.title() + "Copy";
159 if ( h.axis().isFixedBinning() ) {
160 m_rep.reset(
161 new TH1D( title.c_str(), title.c_str(), h.axis().bins(), h.axis().lowerEdge(), h.axis().upperEdge() ) );
162 } else {
163 Edges e;
164 for ( int i = 0; i < h.axis().bins(); ++i ) { e.push_back( h.axis().binLowerEdge( i ) ); }
165 // add also upperedges at the end
166 e.push_back( h.axis().upperEdge() );
167 m_rep.reset( new TH1D( title.c_str(), title.c_str(), e.size() - 1, &e.front() ) );
168 }
169 m_axis.initialize( m_rep->GetXaxis(), false );
170 m_rep->Sumw2();
171 m_sumEntries = 0;
172 m_sumwx = 0;
173 // sumw
174 double sumw = h.sumBinHeights();
175 // sumw2
176 double sumw2 = 0;
177 if ( h.equivalentBinEntries() != 0 ) sumw2 = ( sumw * sumw ) / h.equivalentBinEntries();
178
179 double sumwx = h.mean() * h.sumBinHeights();
180 double sumwx2 = ( h.mean() * h.mean() + h.rms() * h.rms() ) * h.sumBinHeights();
181
182 // copy the contents in
183 for ( int i = -2; i < axis().bins(); ++i ) {
184 // root binning starts from one !
185 m_rep->SetBinContent( rIndex( i ), h.binHeight( i ) );
186 m_rep->SetBinError( rIndex( i ), h.binError( i ) );
187 }
188 // need to do set entries after setting contents otherwise root will recalulate them
189 // taking into account how many time SetBinContents() has been called
190 m_rep->SetEntries( h.allEntries() );
191 // stat vector
192 std::vector<double> stat( 11 );
193 stat[0] = sumw;
194 stat[1] = sumw2;
195 stat[2] = sumwx;
196 stat[3] = sumwx2;
197 m_rep->PutStats( &stat.front() );
198 }
199
200 #ifdef __ICC
201 // re-enable icc remark #1572
202 # pragma warning( pop )
203 #endif
204
205 StreamBuffer& Gaudi::Histogram1D::serialize( StreamBuffer& s ) {
206 // DataObject::serialize(s);
207 std::string title;
208 int size;
209 s >> size;
210 for ( int j = 0; j < size; j++ ) {
211 std::string key, value;
212 s >> key >> value;
213 annotation().addItem( key, value );
214 if ( "Title" == key ) { title = value; }
215 }
216 double lowerEdge, upperEdge, binHeight, binError;
217 int isFixedBinning, bins;
218 s >> isFixedBinning >> bins;
219
220 if ( isFixedBinning ) {
221 s >> lowerEdge >> upperEdge;
222 m_rep.reset( new TH1D( title.c_str(), title.c_str(), bins, lowerEdge, upperEdge ) );
223 } else {
224 Edges edges;
225 edges.resize( bins );
226 for ( int i = 0; i <= bins; ++i ) s >> edges[i];
227 m_rep.reset( new TH1D( title.c_str(), title.c_str(), edges.size() - 1, &edges.front() ) );
228 }
229 m_axis.initialize( m_rep->GetXaxis(), true );
230 m_rep->Sumw2();
231 m_sumEntries = 0;
232 m_sumwx = 0;
233
234 for ( int i = 0; i <= bins + 1; ++i ) {
235 s >> binHeight >> binError;
236 m_rep->SetBinContent( i, binHeight );
237 m_rep->SetBinError( i, binError );
238 }
239 Stat_t allEntries;
240 s >> allEntries;
241 m_rep->SetEntries( allEntries );
242 Stat_t stats[4]; // stats array
243 s >> stats[0] >> stats[1] >> stats[2] >> stats[3];
244 m_rep->PutStats( stats );
245 return s;
246 }
247
248 StreamBuffer& Gaudi::Histogram1D::serialize( StreamBuffer& s ) const {
249 // DataObject::serialize(s);
250 s << static_cast<int>( annotation().size() );
251 for ( int i = 0; i < annotation().size(); i++ ) {
252 s << annotation().key( i );
253 s << annotation().value( i );
254 }
255 const AIDA::IAxis& axis( this->axis() );
256 const int isFixedBinning = axis.isFixedBinning();
257 const int bins = axis.bins();
258 s << isFixedBinning << bins;
259 if ( isFixedBinning ) {
260 s << axis.lowerEdge();
261 } else {
262 for ( int i = 0; i < bins; ++i ) s << axis.binLowerEdge( i );
263 }
264 s << axis.upperEdge();
265 for ( int i = 0; i <= bins + 1; ++i ) s << m_rep->GetBinContent( i ) << m_rep->GetBinError( i );
266
267 s << m_rep->GetEntries();
268 Stat_t stats[4]; // stats array
269 m_rep->GetStats( stats );
270 s << stats[0] << stats[1] << stats[2] << stats[3];
271 return s;
272 }
Gaudi::Histogram1D::fill
bool fill(double x, double weight) override
Fill the Profile1D with a value and the corresponding weight.
Definition: H1D.cpp:149
details::size
constexpr auto size(const T &, Args &&...) noexcept
Definition: AnyDataWrapper.h:14
Gaudi::Histogram1D::initSums
void initSums()
Definition: H1D.cpp:77
Gaudi::Histogram1D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
Definition: H1D.cpp:93
Gaudi::Generic1D::cast
void * cast(const std::string &cl) const override
Manual cast by class name.
Gaudi::Histogram1D::setStatistics
virtual bool setStatistics(int allEntries, double eqBinEntries, double mean, double rms)
set histogram statistics
Definition: H1D.cpp:132
std::vector::front
T front(T... args)
StreamBuffer
The stream buffer is a small object collecting object data.
Definition: StreamBuffer.h:41
Gaudi::Histogram1D::serialize
StreamBuffer & serialize(StreamBuffer &s)
Serialization mechanism, Serialize the object for reading.
Definition: H1D.cpp:205
ObjectFactory.h
Gaudi::Histogram1D::reset
bool reset() override
need to overwrite reset to reset the sums
Definition: H1D.cpp:86
Gaudi::Edges
std::vector< double > Edges
Definition: GaudiPI.h:17
Gaudi::Histogram1D::init
void init(const std::string &title, bool initialize_axis=true)
Definition: H1D.cpp:66
Gaudi::Generic1D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
StreamBuffer.h
GaudiPI.h
std::vector::resize
T resize(T... args)
std::string
STL class.
GaudiPluginService.cpluginsvc.n
n
Definition: cpluginsvc.py:210
std::vector::push_back
T push_back(T... args)
GaudiPython.HistoUtils.mean
mean
Definition: HistoUtils.py:875
std::runtime_error
STL class.
Gaudi::Histogram1D::copyFromAida
void copyFromAida(const AIDA::IHistogram1D &h)
Create new histogram from any AIDA based histogram.
Definition: H1D.cpp:156
std::pair
std::vector< double >
std::lock_guard
Gaudi::Histogram1D::setRms
bool setRms(double rms)
Update histogram RMS.
Definition: H1D.cpp:116
H1D.h
std::string::c_str
T c_str(T... args)
gaudirun.s
string s
Definition: gaudirun.py:316
Gaudi::Generic1D::binEntries
int binEntries(int index) const override
Number of entries in the corresponding bin (ie the number of times fill was called for this bin).
HistogramUtility.h
Gaudi::Histogram1D::setBinContents
virtual bool setBinContents(int i, int entries, double height, double error, double centre)
set bin content (entries and centre are not used )
Definition: H1D.cpp:101
Gaudi::Histogram1D::Histogram1D
Histogram1D()
Standard constructor.
Definition: H1D.cpp:58
Gaudi
Header file for std:chrono::duration-based Counters.
Definition: __init__.py:1
Gaudi::createH1D
std::pair< DataObject *, AIDA::IHistogram1D * > createH1D(const AIDA::IHistogram1D &hist)
Copy constructor.
GaudiPython.HistoUtils.rms
rms
Definition: HistoUtils.py:879
Gaudi::Generic1D
Common AIDA implementation stuff for histograms and profiles using ROOT implementations.
Definition: Generic1D.h:35