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