The Gaudi Framework  v37r0 (b608885e)
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 #ifdef __ICC
129 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
130 // The comparison are meant
131 # pragma warning( push )
132 # pragma warning( disable : 1572 )
133 #endif
134 bool Gaudi::Histogram1D::setRms( double rms ) {
135  m_rep->SetEntries( m_sumEntries );
136  std::vector<double> stat( 11 );
137  // sum weights
138  stat[0] = sumBinHeights();
139  stat[1] = 0;
140  if ( equivalentBinEntries() != 0 ) stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
141  stat[2] = m_sumwx;
142  double mean = 0;
143  if ( sumBinHeights() != 0 ) mean = m_sumwx / sumBinHeights();
144  stat[3] = ( mean * mean + rms * rms ) * sumBinHeights();
145  m_rep->PutStats( &stat.front() );
146  return true;
147 }
148 
149 // set histogram statistics
150 bool Gaudi::Histogram1D::setStatistics( int allEntries, double eqBinEntries, double mean, double rms ) {
151  m_rep->SetEntries( allEntries );
152  // fill statistcal vector for Root
153  std::vector<double> stat( 11 );
154  // sum weights
155  stat[0] = sumBinHeights();
156  // sum weights **2
157  stat[1] = 0;
158  if ( eqBinEntries != 0 ) stat[1] = ( sumBinHeights() * sumBinHeights() ) / eqBinEntries;
159  // sum weights * x
160  stat[2] = mean * sumBinHeights();
161  // sum weight * x **2
162  stat[3] = ( mean * mean + rms * rms ) * sumBinHeights();
163  m_rep->PutStats( &stat.front() );
164  return true;
165 }
166 
167 bool Gaudi::Histogram1D::fill( double x, double weight ) {
168  // avoid race conditions when filling the histogram
169  auto guard = std::scoped_lock{ m_fillSerialization };
170  ( weight == 1. ) ? m_rep->Fill( x ) : m_rep->Fill( x, weight );
171  return true;
172 }
173 
174 void Gaudi::Histogram1D::copyFromAida( const AIDA::IHistogram1D& h ) {
175  // implement here the copy
176  std::string title = h.title() + "Copy";
177  if ( h.axis().isFixedBinning() ) {
178  m_rep.reset(
179  new TH1D( title.c_str(), title.c_str(), h.axis().bins(), h.axis().lowerEdge(), h.axis().upperEdge() ) );
180  } else {
181  Edges e;
182  for ( int i = 0; i < h.axis().bins(); ++i ) { e.push_back( h.axis().binLowerEdge( i ) ); }
183  // add also upperedges at the end
184  e.push_back( h.axis().upperEdge() );
185  m_rep.reset( new TH1D( title.c_str(), title.c_str(), e.size() - 1, &e.front() ) );
186  }
187  m_axis.initialize( m_rep->GetXaxis(), false );
188  m_rep->Sumw2();
189  m_sumEntries = 0;
190  m_sumwx = 0;
191  // sumw
192  double sumw = h.sumBinHeights();
193  // sumw2
194  double sumw2 = 0;
195  if ( h.equivalentBinEntries() != 0 ) sumw2 = ( sumw * sumw ) / h.equivalentBinEntries();
196 
197  double sumwx = h.mean() * h.sumBinHeights();
198  double sumwx2 = ( h.mean() * h.mean() + h.rms() * h.rms() ) * h.sumBinHeights();
199 
200  // copy the contents in
201  for ( int i = -2; i < axis().bins(); ++i ) {
202  // root binning starts from one !
203  m_rep->SetBinContent( rIndex( i ), h.binHeight( i ) );
204  m_rep->SetBinError( rIndex( i ), h.binError( i ) );
205  }
206  // need to do set entries after setting contents otherwise root will recalulate them
207  // taking into account how many time SetBinContents() has been called
208  m_rep->SetEntries( h.allEntries() );
209  // stat vector
210  std::vector<double> stat( 11 );
211  stat[0] = sumw;
212  stat[1] = sumw2;
213  stat[2] = sumwx;
214  stat[3] = sumwx2;
215  m_rep->PutStats( &stat.front() );
216 }
217 
218 #ifdef __ICC
219 // re-enable icc remark #1572
220 # pragma warning( pop )
221 #endif
222 
223 StreamBuffer& Gaudi::Histogram1D::serialize( StreamBuffer& s ) {
224  // DataObject::serialize(s);
225  std::string title;
226  int size;
227  s >> size;
228  for ( int j = 0; j < size; j++ ) {
229  std::string key, value;
230  s >> key >> value;
231  if ( !annotation().addItem( key, value ) ) { annotation().setValue( key, value ); };
232  if ( "Title" == key ) { title = value; }
233  }
234  double lowerEdge, upperEdge, binHeight, binError;
235  int isFixedBinning, bins;
236  s >> isFixedBinning >> bins;
237 
238  if ( isFixedBinning ) {
239  s >> lowerEdge >> upperEdge;
240  m_rep.reset( new TH1D( title.c_str(), title.c_str(), bins, lowerEdge, upperEdge ) );
241  } else {
242  Edges edges;
243  edges.resize( bins );
244  for ( int i = 0; i <= bins; ++i ) s >> edges[i];
245  m_rep.reset( new TH1D( title.c_str(), title.c_str(), edges.size() - 1, &edges.front() ) );
246  }
247  m_axis.initialize( m_rep->GetXaxis(), true );
248  m_rep->Sumw2();
249  m_sumEntries = 0;
250  m_sumwx = 0;
251 
252  for ( int i = 0; i <= bins + 1; ++i ) {
253  s >> binHeight >> binError;
254  m_rep->SetBinContent( i, binHeight );
255  m_rep->SetBinError( i, binError );
256  }
257  Stat_t allEntries;
258  s >> allEntries;
259  m_rep->SetEntries( allEntries );
260  Stat_t stats[4]; // stats array
261  s >> stats[0] >> stats[1] >> stats[2] >> stats[3];
262  m_rep->PutStats( stats );
263  return s;
264 }
265 
266 StreamBuffer& Gaudi::Histogram1D::serialize( StreamBuffer& s ) const {
267  // DataObject::serialize(s);
268  s << static_cast<int>( annotation().size() );
269  for ( int i = 0; i < annotation().size(); i++ ) {
270  s << annotation().key( i );
271  s << annotation().value( i );
272  }
273  const AIDA::IAxis& axis( this->axis() );
274  const int isFixedBinning = axis.isFixedBinning();
275  const int bins = axis.bins();
276  s << isFixedBinning << bins;
277  if ( isFixedBinning ) {
278  s << axis.lowerEdge();
279  } else {
280  for ( int i = 0; i < bins; ++i ) s << axis.binLowerEdge( i );
281  }
282  s << axis.upperEdge();
283  for ( int i = 0; i <= bins + 1; ++i ) s << m_rep->GetBinContent( i ) << m_rep->GetBinError( i );
284 
285  s << m_rep->GetEntries();
286  Stat_t stats[4]; // stats array
287  m_rep->GetStats( stats );
288  s << stats[0] << stats[1] << stats[2] << stats[3];
289  return s;
290 }
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:167
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:223
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:150
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:134
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:92
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
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:174
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