The Gaudi Framework  v37r1 (a7f61348)
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 <Gaudi/MonitoringHub.h>
24 #include <GaudiCommonSvc/H1D.h>
25 #include <GaudiCommonSvc/HistogramUtility.h>
26 #include <GaudiKernel/ObjectFactory.h>
27 #include <GaudiKernel/StreamBuffer.h>
28 
29 std::pair<DataObject*, AIDA::IHistogram1D*> Gaudi::createH1D( ISvcLocator* svcLocator, const std::string& path,
30  const std::string& title, int nBins, double xlow,
31  double xup ) {
32  auto p = new Histogram1D( new TH1D( title.c_str(), title.c_str(), nBins, xlow, xup ) );
33  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
34  return { p, p };
35 }
36 
37 std::pair<DataObject*, AIDA::IHistogram1D*> Gaudi::createH1D( ISvcLocator* svcLocator, const std::string& path,
38  const std::string& title, const Edges& e ) {
39  auto p = new Histogram1D( new TH1D( title.c_str(), title.c_str(), e.size() - 1, &e.front() ) );
40  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
41  return { p, p };
42 }
43 
44 std::pair<DataObject*, AIDA::IHistogram1D*> Gaudi::createH1D( ISvcLocator* svcLocator, const std::string& path,
45  const AIDA::IHistogram1D& hist ) {
46  TH1D* h = getRepresentation<AIDA::IHistogram1D, TH1D>( hist );
47  auto n = ( h ? new Histogram1D( new TH1D( *h ) ) : nullptr );
48  if ( n ) { svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *n ); }
49  return { n, n };
50 }
51 
52 namespace Gaudi {
53 
54  template <>
55  void* Generic1D<AIDA::IHistogram1D, TH1D>::cast( const std::string& className ) const {
56  if ( className == "AIDA::IHistogram1D" ) return const_cast<AIDA::IHistogram1D*>( (AIDA::IHistogram1D*)this );
57  if ( className == "AIDA::IHistogram" ) return const_cast<AIDA::IHistogram*>( (AIDA::IHistogram*)this );
58  return nullptr;
59  }
60 
61  template <>
62  int Generic1D<AIDA::IHistogram1D, TH1D>::binEntries( int index ) const {
63  if ( binHeight( index ) <= 0 ) return 0;
64  double xx = binHeight( index ) / binError( index );
65  return int( xx * xx + 0.5 );
66  }
67 
68  template <>
69  void Generic1D<AIDA::IHistogram1D, TH1D>::adoptRepresentation( TObject* rep ) {
70  TH1D* imp = dynamic_cast<TH1D*>( rep );
71  if ( !imp ) throw std::runtime_error( "Cannot adopt native histogram representation." );
72  m_rep.reset( imp );
73  }
74 } // namespace Gaudi
75 
76 Gaudi::Histogram1D::Histogram1D() : Base( new TH1D() ) { init( "", false ); }
77 
78 Gaudi::Histogram1D::Histogram1D( TH1D* rep ) {
79  m_rep.reset( rep );
80  init( m_rep->GetTitle() );
81  initSums();
82 }
83 
84 void Gaudi::Histogram1D::init( const std::string& title, bool initialize_axis ) {
85  m_classType = "IHistogram1D";
86  if ( initialize_axis ) { m_axis.initialize( m_rep->GetXaxis(), false ); }
87  const TArrayD* a = m_rep->GetSumw2();
88  if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
89  setTitle( title );
90  m_rep->SetDirectory( nullptr );
91  m_sumEntries = 0;
92  m_sumwx = 0;
93 }
94 
95 void Gaudi::Histogram1D::initSums() {
96  m_sumwx = 0;
97  m_sumEntries = 0;
98  for ( int i = 1, n = m_rep->GetNbinsX(); i <= n; ++i ) {
99  m_sumwx += m_rep->GetBinContent( i ) * m_rep->GetBinCenter( i );
100  m_sumEntries += m_rep->GetBinContent( i );
101  }
102 }
103 
104 bool Gaudi::Histogram1D::reset() {
105  m_sumwx = 0;
106  m_sumEntries = 0;
107  return Base::reset();
108 }
109 
111 void Gaudi::Histogram1D::adoptRepresentation( TObject* rep ) {
112  Gaudi::Generic1D<AIDA::IHistogram1D, TH1D>::adoptRepresentation( rep );
113  if ( m_rep ) {
114  init( m_rep->GetTitle() );
115  initSums();
116  }
117 }
118 
119 bool Gaudi::Histogram1D::setBinContents( int i, int entries, double height, double error, double centre ) {
120  m_rep->SetBinContent( rIndex( i ), height );
121  m_rep->SetBinError( rIndex( i ), error );
122  // accumulate sumwx for in range bins
123  if ( i != AIDA::IAxis::UNDERFLOW_BIN && i != AIDA::IAxis::OVERFLOW_BIN ) m_sumwx += centre * height;
124  m_sumEntries += entries;
125  return true;
126 }
127 
128 bool Gaudi::Histogram1D::setRms( double rms ) {
129  m_rep->SetEntries( m_sumEntries );
130  std::vector<double> stat( 11 );
131  // sum weights
132  stat[0] = sumBinHeights();
133  stat[1] = 0;
134  if ( abs( equivalentBinEntries() ) > std::numeric_limits<double>::epsilon() )
135  stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
136  stat[2] = m_sumwx;
137  double mean = 0;
138  if ( abs( sumBinHeights() ) > std::numeric_limits<double>::epsilon() ) mean = m_sumwx / sumBinHeights();
139  stat[3] = ( mean * mean + rms * rms ) * sumBinHeights();
140  m_rep->PutStats( &stat.front() );
141  return true;
142 }
143 
144 // set histogram statistics
145 bool Gaudi::Histogram1D::setStatistics( int allEntries, double eqBinEntries, double mean, double rms ) {
146  m_rep->SetEntries( allEntries );
147  // fill statistcal vector for Root
148  std::vector<double> stat( 11 );
149  // sum weights
150  stat[0] = sumBinHeights();
151  // sum weights **2
152  stat[1] = 0;
153  if ( abs( eqBinEntries ) > std::numeric_limits<double>::epsilon() )
154  stat[1] = ( sumBinHeights() * sumBinHeights() ) / eqBinEntries;
155  // sum weights * x
156  stat[2] = mean * sumBinHeights();
157  // sum weight * x **2
158  stat[3] = ( mean * mean + rms * rms ) * sumBinHeights();
159  m_rep->PutStats( &stat.front() );
160  return true;
161 }
162 
163 bool Gaudi::Histogram1D::fill( double x, double weight ) {
164  // avoid race conditions when filling the histogram
165  auto guard = std::scoped_lock{ m_fillSerialization };
166  m_rep->Fill( x, weight );
167  return true;
168 }
169 
170 void Gaudi::Histogram1D::copyFromAida( const AIDA::IHistogram1D& h ) {
171  // implement here the copy
172  std::string title = h.title() + "Copy";
173  if ( h.axis().isFixedBinning() ) {
174  m_rep.reset(
175  new TH1D( title.c_str(), title.c_str(), h.axis().bins(), h.axis().lowerEdge(), h.axis().upperEdge() ) );
176  } else {
177  Edges e;
178  for ( int i = 0; i < h.axis().bins(); ++i ) { e.push_back( h.axis().binLowerEdge( i ) ); }
179  // add also upperedges at the end
180  e.push_back( h.axis().upperEdge() );
181  m_rep.reset( new TH1D( title.c_str(), title.c_str(), e.size() - 1, &e.front() ) );
182  }
183  m_axis.initialize( m_rep->GetXaxis(), false );
184  m_rep->Sumw2();
185  m_sumEntries = 0;
186  m_sumwx = 0;
187  // sumw
188  double sumw = h.sumBinHeights();
189  // sumw2
190  double sumw2 = 0;
191  if ( abs( h.equivalentBinEntries() ) > std::numeric_limits<double>::epsilon() )
192  sumw2 = ( sumw * sumw ) / h.equivalentBinEntries();
193 
194  double sumwx = h.mean() * h.sumBinHeights();
195  double sumwx2 = ( h.mean() * h.mean() + h.rms() * h.rms() ) * h.sumBinHeights();
196 
197  // copy the contents in
198  for ( int i = -2; i < axis().bins(); ++i ) {
199  // root binning starts from one !
200  m_rep->SetBinContent( rIndex( i ), h.binHeight( i ) );
201  m_rep->SetBinError( rIndex( i ), h.binError( i ) );
202  }
203  // need to do set entries after setting contents otherwise root will recalulate them
204  // taking into account how many time SetBinContents() has been called
205  m_rep->SetEntries( h.allEntries() );
206  // stat vector
207  std::vector<double> stat( 11 );
208  stat[0] = sumw;
209  stat[1] = sumw2;
210  stat[2] = sumwx;
211  stat[3] = sumwx2;
212  m_rep->PutStats( &stat.front() );
213 }
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 }
std::vector::resize
T resize(T... args)
ISvcLocator::monitoringHub
Gaudi::Monitoring::Hub & monitoringHub()
Definition: ISvcLocator.cpp:26
std::string
STL class.
details::size
constexpr auto size(const T &, Args &&...) noexcept
Definition: AnyDataWrapper.h:22
GaudiAlg.HistoUtils.path
path
Definition: HistoUtils.py:961
MonitoringHub.h
std::pair
Gaudi::Histogram1D::Histogram1D
Histogram1D()
Standard constructor.
Definition: H1D.cpp:76
gaudirun.s
string s
Definition: gaudirun.py:348
Gaudi::Histogram1D::fill
bool fill(double x, double weight) override
Fill the Profile1D with a value and the corresponding weight.
Definition: H1D.cpp:163
std::vector< double >
Gaudi::Generic1D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
ISvcLocator
Definition: ISvcLocator.h:46
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:145
Gaudi::svcLocator
GAUDI_API ISvcLocator * svcLocator()
Gaudi::Monitoring::reset
void reset(T &)
default (empty) implementation of reset method for types stored into an entity
Definition: MonitoringHub.h:75
ProduceConsume.j
j
Definition: ProduceConsume.py:101
Gaudi::Histogram1D::setRms
bool setRms(double rms)
Update histogram RMS.
Definition: H1D.cpp:128
std::string::c_str
T c_str(T... args)
AlgSequencer.h
h
Definition: AlgSequencer.py:32
Gaudi::Histogram1D::reset
bool reset() override
Definition: H1D.cpp:104
GaudiPython.Bindings.nullptr
nullptr
Definition: Bindings.py:88
Gaudi::Monitoring::Hub::registerEntity
void registerEntity(std::string c, std::string n, std::string t, T &ent)
Definition: MonitoringHub.h:174
std::runtime_error
STL class.
Gaudi::Histogram1D::init
void init(const std::string &title, bool initialize_axis=true)
Definition: H1D.cpp:84
GaudiAlg.HistoUtils.rms
rms
Definition: HistoUtils.py:908
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:119
Gaudi
Header file for std:chrono::duration-based Counters.
Definition: __init__.py:1
GaudiPluginService.cpluginsvc.n
n
Definition: cpluginsvc.py:235
std::numeric_limits::epsilon
T epsilon(T... args)
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:28
Gaudi::Histogram1D::copyFromAida
void copyFromAida(const AIDA::IHistogram1D &h)
Create new histogram from any AIDA based histogram.
Definition: H1D.cpp:170
Gaudi::Histogram1D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
Definition: H1D.cpp:111
Gaudi::Generic1D
Definition: Generic1D.h:46
GaudiAlg.HistoUtils.mean
mean
Definition: HistoUtils.py:904
Gaudi::createH1D
std::pair< DataObject *, AIDA::IHistogram1D * > createH1D(ISvcLocator *svcLocator, const std::string &path, const AIDA::IHistogram1D &hist)
Copy constructor.
Gaudi::Histogram1D::initSums
void initSums()
Definition: H1D.cpp:95
ProduceConsume.key
key
Definition: ProduceConsume.py:81
std::numeric_limits