Loading [MathJax]/extensions/tex2jax.js
The Gaudi Framework  v36r16 (ea80daf8)
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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 #include <Gaudi/Histograming/Sink/Utils.h>
30 
31 std::pair<DataObject*, AIDA::IHistogram1D*> Gaudi::createH1D( ISvcLocator* svcLocator, const std::string& path,
32  const std::string& title, int nBins, double xlow,
33  double xup ) {
34  auto p = new Histogram1D( new TH1D( title.c_str(), title.c_str(), nBins, xlow, xup ) );
35  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
36  return { p, p };
37 }
38 
39 std::pair<DataObject*, AIDA::IHistogram1D*> Gaudi::createH1D( ISvcLocator* svcLocator, const std::string& path,
40  const std::string& title, const Edges& e ) {
41  auto p = new Histogram1D( new TH1D( title.c_str(), title.c_str(), e.size() - 1, &e.front() ) );
42  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
43  return { p, p };
44 }
45 
46 std::pair<DataObject*, AIDA::IHistogram1D*> Gaudi::createH1D( ISvcLocator* svcLocator, const std::string& path,
47  const AIDA::IHistogram1D& hist ) {
48  TH1D* h = getRepresentation<AIDA::IHistogram1D, TH1D>( hist );
49  auto n = ( h ? new Histogram1D( new TH1D( *h ) ) : nullptr );
50  if ( n ) { svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *n ); }
51  return { n, n };
52 }
53 
54 namespace Gaudi {
55 
56  template <>
57  void* Generic1D<AIDA::IHistogram1D, TH1D>::cast( const std::string& className ) const {
58  if ( className == "AIDA::IHistogram1D" ) return const_cast<AIDA::IHistogram1D*>( (AIDA::IHistogram1D*)this );
59  if ( className == "AIDA::IHistogram" ) return const_cast<AIDA::IHistogram*>( (AIDA::IHistogram*)this );
60  return nullptr;
61  }
62 
63  template <>
64  int Generic1D<AIDA::IHistogram1D, TH1D>::binEntries( int index ) const {
65  if ( binHeight( index ) <= 0 ) return 0;
66  double xx = binHeight( index ) / binError( index );
67  return int( xx * xx + 0.5 );
68  }
69 
70  template <>
71  void Generic1D<AIDA::IHistogram1D, TH1D>::adoptRepresentation( TObject* rep ) {
72  TH1D* imp = dynamic_cast<TH1D*>( rep );
73  if ( !imp ) throw std::runtime_error( "Cannot adopt native histogram representation." );
74  m_rep.reset( imp );
75  }
76 } // namespace Gaudi
77 
78 Gaudi::Histogram1D::Histogram1D() : Base( new TH1D() ) { init( "", false ); }
79 
80 Gaudi::Histogram1D::Histogram1D( TH1D* rep ) {
81  m_rep.reset( rep );
82  init( m_rep->GetTitle() );
83  initSums();
84 }
85 
86 void Gaudi::Histogram1D::init( const std::string& title, bool initialize_axis ) {
87  m_classType = "IHistogram1D";
88  if ( initialize_axis ) { m_axis.initialize( m_rep->GetXaxis(), false ); }
89  const TArrayD* a = m_rep->GetSumw2();
90  if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
91  setTitle( title );
92  m_rep->SetDirectory( nullptr );
93  m_sumEntries = 0;
94  m_sumwx = 0;
95 }
96 
97 void Gaudi::Histogram1D::initSums() {
98  m_sumwx = 0;
99  m_sumEntries = 0;
100  for ( int i = 1, n = m_rep->GetNbinsX(); i <= n; ++i ) {
101  m_sumwx += m_rep->GetBinContent( i ) * m_rep->GetBinCenter( i );
102  m_sumEntries += m_rep->GetBinContent( i );
103  }
104 }
105 
106 bool Gaudi::Histogram1D::reset() {
107  m_sumwx = 0;
108  m_sumEntries = 0;
109  return Base::reset();
110 }
111 
112 nlohmann::json Gaudi::Histogram1D::toJSON() const { return *m_rep.get(); }
113 
115 void Gaudi::Histogram1D::adoptRepresentation( TObject* rep ) {
116  Gaudi::Generic1D<AIDA::IHistogram1D, TH1D>::adoptRepresentation( rep );
117  if ( m_rep ) {
118  init( m_rep->GetTitle() );
119  initSums();
120  }
121 }
122 
123 bool Gaudi::Histogram1D::setBinContents( int i, int entries, double height, double error, double centre ) {
124  m_rep->SetBinContent( rIndex( i ), height );
125  m_rep->SetBinError( rIndex( i ), error );
126  // accumulate sumwx for in range bins
127  if ( i != AIDA::IAxis::UNDERFLOW_BIN && i != AIDA::IAxis::OVERFLOW_BIN ) m_sumwx += centre * height;
128  m_sumEntries += entries;
129  return true;
130 }
131 
132 #ifdef __ICC
133 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
134 // The comparison are meant
135 # pragma warning( push )
136 # pragma warning( disable : 1572 )
137 #endif
138 bool Gaudi::Histogram1D::setRms( double rms ) {
139  m_rep->SetEntries( m_sumEntries );
140  std::vector<double> stat( 11 );
141  // sum weights
142  stat[0] = sumBinHeights();
143  stat[1] = 0;
144  if ( equivalentBinEntries() != 0 ) stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
145  stat[2] = m_sumwx;
146  double mean = 0;
147  if ( sumBinHeights() != 0 ) mean = m_sumwx / sumBinHeights();
148  stat[3] = ( mean * mean + rms * rms ) * sumBinHeights();
149  m_rep->PutStats( &stat.front() );
150  return true;
151 }
152 
153 // set histogram statistics
154 bool Gaudi::Histogram1D::setStatistics( int allEntries, double eqBinEntries, double mean, double rms ) {
155  m_rep->SetEntries( allEntries );
156  // fill statistcal vector for Root
157  std::vector<double> stat( 11 );
158  // sum weights
159  stat[0] = sumBinHeights();
160  // sum weights **2
161  stat[1] = 0;
162  if ( eqBinEntries != 0 ) stat[1] = ( sumBinHeights() * sumBinHeights() ) / eqBinEntries;
163  // sum weights * x
164  stat[2] = mean * sumBinHeights();
165  // sum weight * x **2
166  stat[3] = ( mean * mean + rms * rms ) * sumBinHeights();
167  m_rep->PutStats( &stat.front() );
168  return true;
169 }
170 
171 bool Gaudi::Histogram1D::fill( double x, double weight ) {
172  // avoid race conditions when filling the histogram
173  auto guard = std::scoped_lock{ m_fillSerialization };
174  ( weight == 1. ) ? m_rep->Fill( x ) : m_rep->Fill( x, weight );
175  return true;
176 }
177 
178 void Gaudi::Histogram1D::copyFromAida( const AIDA::IHistogram1D& h ) {
179  // implement here the copy
180  std::string title = h.title() + "Copy";
181  if ( h.axis().isFixedBinning() ) {
182  m_rep.reset(
183  new TH1D( title.c_str(), title.c_str(), h.axis().bins(), h.axis().lowerEdge(), h.axis().upperEdge() ) );
184  } else {
185  Edges e;
186  for ( int i = 0; i < h.axis().bins(); ++i ) { e.push_back( h.axis().binLowerEdge( i ) ); }
187  // add also upperedges at the end
188  e.push_back( h.axis().upperEdge() );
189  m_rep.reset( new TH1D( title.c_str(), title.c_str(), e.size() - 1, &e.front() ) );
190  }
191  m_axis.initialize( m_rep->GetXaxis(), false );
192  m_rep->Sumw2();
193  m_sumEntries = 0;
194  m_sumwx = 0;
195  // sumw
196  double sumw = h.sumBinHeights();
197  // sumw2
198  double sumw2 = 0;
199  if ( h.equivalentBinEntries() != 0 ) sumw2 = ( sumw * sumw ) / h.equivalentBinEntries();
200 
201  double sumwx = h.mean() * h.sumBinHeights();
202  double sumwx2 = ( h.mean() * h.mean() + h.rms() * h.rms() ) * h.sumBinHeights();
203 
204  // copy the contents in
205  for ( int i = -2; i < axis().bins(); ++i ) {
206  // root binning starts from one !
207  m_rep->SetBinContent( rIndex( i ), h.binHeight( i ) );
208  m_rep->SetBinError( rIndex( i ), h.binError( i ) );
209  }
210  // need to do set entries after setting contents otherwise root will recalulate them
211  // taking into account how many time SetBinContents() has been called
212  m_rep->SetEntries( h.allEntries() );
213  // stat vector
214  std::vector<double> stat( 11 );
215  stat[0] = sumw;
216  stat[1] = sumw2;
217  stat[2] = sumwx;
218  stat[3] = sumwx2;
219  m_rep->PutStats( &stat.front() );
220 }
221 
222 #ifdef __ICC
223 // re-enable icc remark #1572
224 # pragma warning( pop )
225 #endif
226 
227 StreamBuffer& Gaudi::Histogram1D::serialize( StreamBuffer& s ) {
228  // DataObject::serialize(s);
229  std::string title;
230  int size;
231  s >> size;
232  for ( int j = 0; j < size; j++ ) {
233  std::string key, value;
234  s >> key >> value;
235  if ( !annotation().addItem( key, value ) ) { annotation().setValue( key, value ); };
236  if ( "Title" == key ) { title = value; }
237  }
238  double lowerEdge, upperEdge, binHeight, binError;
239  int isFixedBinning, bins;
240  s >> isFixedBinning >> bins;
241 
242  if ( isFixedBinning ) {
243  s >> lowerEdge >> upperEdge;
244  m_rep.reset( new TH1D( title.c_str(), title.c_str(), bins, lowerEdge, upperEdge ) );
245  } else {
246  Edges edges;
247  edges.resize( bins );
248  for ( int i = 0; i <= bins; ++i ) s >> edges[i];
249  m_rep.reset( new TH1D( title.c_str(), title.c_str(), edges.size() - 1, &edges.front() ) );
250  }
251  m_axis.initialize( m_rep->GetXaxis(), true );
252  m_rep->Sumw2();
253  m_sumEntries = 0;
254  m_sumwx = 0;
255 
256  for ( int i = 0; i <= bins + 1; ++i ) {
257  s >> binHeight >> binError;
258  m_rep->SetBinContent( i, binHeight );
259  m_rep->SetBinError( i, binError );
260  }
261  Stat_t allEntries;
262  s >> allEntries;
263  m_rep->SetEntries( allEntries );
264  Stat_t stats[4]; // stats array
265  s >> stats[0] >> stats[1] >> stats[2] >> stats[3];
266  m_rep->PutStats( stats );
267  return s;
268 }
269 
270 StreamBuffer& Gaudi::Histogram1D::serialize( StreamBuffer& s ) const {
271  // DataObject::serialize(s);
272  s << static_cast<int>( annotation().size() );
273  for ( int i = 0; i < annotation().size(); i++ ) {
274  s << annotation().key( i );
275  s << annotation().value( i );
276  }
277  const AIDA::IAxis& axis( this->axis() );
278  const int isFixedBinning = axis.isFixedBinning();
279  const int bins = axis.bins();
280  s << isFixedBinning << bins;
281  if ( isFixedBinning ) {
282  s << axis.lowerEdge();
283  } else {
284  for ( int i = 0; i < bins; ++i ) s << axis.binLowerEdge( i );
285  }
286  s << axis.upperEdge();
287  for ( int i = 0; i <= bins + 1; ++i ) s << m_rep->GetBinContent( i ) << m_rep->GetBinError( i );
288 
289  s << m_rep->GetEntries();
290  Stat_t stats[4]; // stats array
291  m_rep->GetStats( stats );
292  s << stats[0] << stats[1] << stats[2] << stats[3];
293  return s;
294 }
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
AtlasMCRecoFullPrecedenceDump.path
path
Definition: AtlasMCRecoFullPrecedenceDump.py:49
MonitoringHub.h
std::pair
Gaudi::Histogram1D::Histogram1D
Histogram1D()
Standard constructor.
Definition: H1D.cpp:78
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:171
std::vector< double >
Gaudi::Generic1D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
ISvcLocator
Definition: ISvcLocator.h:46
jsonFromLHCbLog.json
json
Definition: jsonFromLHCbLog.py:87
StreamBuffer.h
ObjectFactory.h
Gaudi::Histogram1D::serialize
StreamBuffer & serialize(StreamBuffer &s)
Serialization mechanism, Serialize the object for reading.
Definition: H1D.cpp:227
std::vector::front
T front(T... args)
StreamBuffer
Definition: StreamBuffer.h:52
Utils.h
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:154
Gaudi::svcLocator
GAUDI_API ISvcLocator * svcLocator()
ProduceConsume.j
j
Definition: ProduceConsume.py:101
Gaudi::Histogram1D::setRms
bool setRms(double rms)
Update histogram RMS.
Definition: H1D.cpp:138
std::string::c_str
T c_str(T... args)
AlgSequencer.h
h
Definition: AlgSequencer.py:32
Gaudi::Histogram1D::reset
bool reset() override
need to overwrite reset to reset the sums
Definition: H1D.cpp:106
GaudiPython.Bindings.nullptr
nullptr
Definition: Bindings.py:92
GaudiPython.HistoUtils.rms
rms
Definition: HistoUtils.py:908
Gaudi::Monitoring::Hub::registerEntity
void registerEntity(std::string c, std::string n, std::string t, T &ent)
Definition: MonitoringHub.h:167
std::runtime_error
STL class.
Gaudi::Histogram1D::init
void init(const std::string &title, bool initialize_axis=true)
Definition: H1D.cpp:86
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:123
GaudiPython.HistoUtils.mean
mean
Definition: HistoUtils.py:904
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::toJSON
nlohmann::json toJSON() const
dumps Histogram to json data
Definition: H1D.cpp:112
Gaudi::Histogram1D::copyFromAida
void copyFromAida(const AIDA::IHistogram1D &h)
Create new histogram from any AIDA based histogram.
Definition: H1D.cpp:178
Gaudi::Histogram1D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
Definition: H1D.cpp:115
Gaudi::Generic1D
Definition: Generic1D.h:46
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:97
ProduceConsume.key
key
Definition: ProduceConsume.py:81