Loading [MathJax]/extensions/tex2jax.js
The Gaudi Framework  v36r13 (995e4364)
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
H2D.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/H2D.h>
26 #include <GaudiCommonSvc/HistogramUtility.h>
27 #include <GaudiCommonSvc/P1D.h>
28 #include <GaudiKernel/ObjectFactory.h>
29 #include <TH2D.h>
30 #include <TProfile.h>
31 #include <array>
32 
33 #include <Gaudi/Histograming/Sink/Utils.h>
34 
35 namespace {
36  using AIDA::IHistogram1D;
37  using AIDA::IHistogram2D;
38  using AIDA::IProfile1D;
39 } // namespace
40 
41 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( ISvcLocator* svcLocator, const std::string& path,
42  const std::string& title, int binsX, double iminX, double imaxX,
43  int binsY, double iminY, double imaxY ) {
44  auto p = new Histogram2D( new TH2D( title.c_str(), title.c_str(), binsX, iminX, imaxX, binsY, iminY, imaxY ) );
45  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
46  return { p, p };
47 }
48 
49 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( ISvcLocator* svcLocator, const std::string& path,
50  const std::string& title, const Edges& eX, const Edges& eY ) {
51  auto p = new Histogram2D(
52  new TH2D( title.c_str(), title.c_str(), eX.size() - 1, &eX.front(), eY.size() - 1, &eY.front() ) );
53  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
54  return { p, p };
55 }
56 
57 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( ISvcLocator* svcLocator, const std::string& path, TH2D* rep ) {
58  auto p = new Histogram2D( rep );
59  svcLocator->monitoringHub().registerEntity( "", path, "histogram:Histogram:double", *p );
60  return { p, p };
61 }
62 
63 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( ISvcLocator* svcLocator, const std::string& path,
64  const IHistogram2D& hist ) {
65  TH2D* h = getRepresentation<AIDA::IHistogram2D, TH2D>( hist );
66  Histogram2D* n = h ? new Histogram2D( new TH2D( *h ) ) : nullptr;
67  if ( n ) { svcLocator->monitoringHub().registerEntity( path, h->GetName(), "histogram:Histogram:double", *n ); }
68  return { n, n };
69 }
70 
71 std::pair<DataObject*, IHistogram1D*> Gaudi::slice1DX( const std::string& nam, const IHistogram2D& hist, int first,
72  int last ) {
73  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
74  TH1D* t = r ? r->ProjectionX( "_px", first, last, "e" ) : nullptr;
75  if ( t ) t->SetName( nam.c_str() );
76  Histogram1D* p = ( t ? new Histogram1D( t ) : nullptr );
77  return { p, p };
78 }
79 
80 std::pair<DataObject*, IHistogram1D*> Gaudi::slice1DY( const std::string& nam, const IHistogram2D& hist, int first,
81  int last ) {
82  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
83  TH1D* t = r ? r->ProjectionY( "_py", first, last, "e" ) : nullptr;
84  if ( t ) t->SetName( nam.c_str() );
85  Histogram1D* p = ( t ? new Histogram1D( t ) : nullptr );
86  return { p, p };
87 }
88 
89 std::pair<DataObject*, IHistogram1D*> Gaudi::project1DY( const std::string& nam, const IHistogram2D& hist, int first,
90  int last ) {
91  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
92  TH1D* t = r ? r->ProjectionY( "_px", first, last, "e" ) : nullptr;
93  if ( t ) t->SetName( nam.c_str() );
94  Histogram1D* p = ( t ? new Histogram1D( t ) : nullptr );
95  return { p, p };
96 }
97 
98 std::pair<DataObject*, IProfile1D*> Gaudi::profile1DX( const std::string& nam, const IHistogram2D& hist, int first,
99  int last ) {
100  TH2* r = Gaudi::getRepresentation<IHistogram2D, TH2>( hist );
101  TProfile* t = r ? r->ProfileX( "_pfx", first, last, "e" ) : nullptr;
102  if ( t ) t->SetName( nam.c_str() );
103  Profile1D* p = ( t ? new Profile1D( t ) : nullptr );
104  return { p, p };
105 }
106 
107 std::pair<DataObject*, IProfile1D*> Gaudi::profile1DY( const std::string& nam, const IHistogram2D& hist, int first,
108  int last ) {
109  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
110  TProfile* t = r ? r->ProfileY( "_pfx", first, last, "e" ) : nullptr;
111  if ( t ) t->SetName( nam.c_str() );
112  Profile1D* p = ( t ? new Profile1D( t ) : nullptr );
113  return { p, p };
114 }
115 
116 namespace Gaudi {
117  template <>
118  void* Generic2D<IHistogram2D, TH2D>::cast( const std::string& className ) const {
119  if ( className == "AIDA::IHistogram2D" )
120  return (IHistogram2D*)this;
121  else if ( className == "AIDA::IHistogram" )
122  return (IHistogram*)this;
123  return nullptr;
124  }
125 
126  template <>
127  int Generic2D<IHistogram2D, TH2D>::binEntries( int indexX, int indexY ) const {
128  if ( binHeight( indexX, indexY ) <= 0 ) return 0;
129  double xx = binHeight( indexX, indexY ) / binError( indexX, indexY );
130  return int( xx * xx + 0.5 );
131  }
132 
133  template <>
134  void Generic2D<AIDA::IHistogram2D, TH2D>::adoptRepresentation( TObject* rep ) {
135  TH2D* imp = dynamic_cast<TH2D*>( rep );
136  if ( !imp ) throw std::runtime_error( "Cannot adopt native histogram representation." );
137  m_rep.reset( imp );
138  m_xAxis.initialize( m_rep->GetXaxis(), true );
139  m_yAxis.initialize( m_rep->GetYaxis(), true );
140  const TArrayD* a = m_rep->GetSumw2();
141  if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
142  setTitle( m_rep->GetTitle() );
143  }
144 } // namespace Gaudi
145 
146 Gaudi::Histogram2D::Histogram2D() : Base( new TH2D() ) {
147  m_rep->Sumw2();
148  m_sumwx = m_sumwy = 0;
149  setTitle( "" );
150  m_rep->SetDirectory( nullptr );
151 }
152 
153 Gaudi::Histogram2D::Histogram2D( TH2D* rep ) {
154  adoptRepresentation( rep );
155  m_sumwx = m_sumwy = 0;
156  m_rep->SetDirectory( nullptr );
157 }
158 
159 bool Gaudi::Histogram2D::setBinContents( int i, int j, int entries, double height, double error, double centreX,
160  double centreY ) {
161  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), height );
162  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), error );
163  // accumulate sumwx for in range bins
164  if ( i >= 0 && j >= 0 ) {
165  m_sumwx += centreX * height;
166  m_sumwy += centreY * height;
167  }
168  m_sumEntries += entries;
169  return true;
170 }
171 
172 bool Gaudi::Histogram2D::reset() {
173  m_sumwx = 0;
174  m_sumwy = 0;
175  return Base::reset();
176 }
177 
178 nlohmann::json Gaudi::Histogram2D::toJSON() const { return *m_rep.get(); }
179 
180 #ifdef __ICC
181 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
182 // The comparison are meant
183 # pragma warning( push )
184 # pragma warning( disable : 1572 )
185 #endif
186 bool Gaudi::Histogram2D::fill( double x, double y, double weight ) {
187  // avoid race conditiosn when filling the histogram
188  auto guard = std::scoped_lock{ m_fillSerialization };
189  ( weight == 1. ) ? m_rep->Fill( x, y ) : m_rep->Fill( x, y, weight );
190  return true;
191 }
192 
193 bool Gaudi::Histogram2D::setRms( double rmsX, double rmsY ) {
194  m_rep->SetEntries( m_sumEntries );
195  std::vector<double> stat( 11 );
196  stat[0] = sumBinHeights();
197  stat[1] = 0;
198  if ( equivalentBinEntries() != 0 ) stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
199  stat[2] = m_sumwx;
200  double meanX = 0;
201  if ( sumBinHeights() != 0 ) meanX = m_sumwx / sumBinHeights();
202  stat[3] = ( meanX * meanX + rmsX * rmsX ) * sumBinHeights();
203  stat[4] = m_sumwy;
204  double meanY = 0;
205  if ( sumBinHeights() != 0 ) meanY = m_sumwy / sumBinHeights();
206  stat[5] = ( meanY * meanY + rmsY * rmsY ) * sumBinHeights();
207  stat[6] = 0;
208  m_rep->PutStats( &stat.front() );
209  return true;
210 }
211 
212 void Gaudi::Histogram2D::copyFromAida( const IHistogram2D& h ) {
213  // implement here the copy
214  std::string titlestr = h.title();
215  const char* title = titlestr.c_str();
216  if ( h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning() )
217  m_rep.reset( new TH2D( title, title, h.xAxis().bins(), h.xAxis().lowerEdge(), h.xAxis().upperEdge(),
218  h.yAxis().bins(), h.yAxis().lowerEdge(), h.yAxis().upperEdge() ) );
219  else {
220  Edges eX, eY;
221  for ( int i = 0; i < h.xAxis().bins(); ++i ) eX.push_back( h.xAxis().binLowerEdge( i ) );
222  // add also upperedges at the end
223  eX.push_back( h.xAxis().upperEdge() );
224  for ( int i = 0; i < h.yAxis().bins(); ++i ) eY.push_back( h.yAxis().binLowerEdge( i ) );
225  // add also upperedges at the end
226  eY.push_back( h.yAxis().upperEdge() );
227  m_rep.reset( new TH2D( title, title, eX.size() - 1, &eX.front(), eY.size() - 1, &eY.front() ) );
228  }
229  m_xAxis.initialize( m_rep->GetXaxis(), true );
230  m_yAxis.initialize( m_rep->GetYaxis(), true );
231  m_rep->Sumw2();
232  m_sumEntries = 0;
233  m_sumwx = 0;
234  m_sumwy = 0;
235  // statistics
236  double sumw = h.sumBinHeights();
237  double sumw2 = 0;
238  if ( h.equivalentBinEntries() != 0 ) sumw2 = ( sumw * sumw ) / h.equivalentBinEntries();
239  double sumwx = h.meanX() * h.sumBinHeights();
240  double sumwx2 = ( h.meanX() * h.meanX() + h.rmsX() * h.rmsX() ) * h.sumBinHeights();
241  double sumwy = h.meanY() * h.sumBinHeights();
242  double sumwy2 = ( h.meanY() * h.meanY() + h.rmsY() * h.rmsY() ) * h.sumBinHeights();
243  double sumwxy = 0;
244 
245  // copy the contents in (AIDA underflow/overflow are -2,-1)
246  for ( int i = -2; i < xAxis().bins(); ++i ) {
247  for ( int j = -2; j < yAxis().bins(); ++j ) {
248  // root binning starts from one !
249  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), h.binHeight( i, j ) );
250  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), h.binError( i, j ) );
251  // calculate statistics
252  if ( i >= 0 && j >= 0 ) { sumwxy += h.binHeight( i, j ) * h.binMeanX( i, j ) * h.binMeanY( i, j ); }
253  }
254  }
255  // need to do set entries after setting contents otherwise root will recalculate them
256  // taking into account how many time SetBinContents() has been called
257  m_rep->SetEntries( h.allEntries() );
258  // fill stat vector
259  std::array<double, 11> stat = { { sumw, sumw2, sumwx, sumwx2, sumwy, sumwy2, sumwxy } };
260  m_rep->PutStats( stat.data() );
261 }
262 
263 #ifdef __ICC
264 // re-enable icc remark #1572
265 # pragma warning( pop )
266 #endif
Gaudi::Histogram2D::reset
bool reset() override
Definition: H2D.cpp:172
Gaudi::profile1DX
std::pair< DataObject *, AIDA::IProfile1D * > profile1DX(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D profile in X from 2D histogram.
Gaudi::profile1DY
std::pair< DataObject *, AIDA::IProfile1D * > profile1DY(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D profile in Y from 2D histogram.
Gaudi::Histogram2D::m_sumwx
double m_sumwx
Definition: H2D.h:56
ISvcLocator::monitoringHub
Gaudi::Monitoring::Hub & monitoringHub()
Definition: ISvcLocator.cpp:26
std::string
STL class.
Gaudi::Histogram2D::toJSON
nlohmann::json toJSON() const
dumps Histogram to json data
Definition: H2D.cpp:178
AtlasMCRecoFullPrecedenceDump.path
path
Definition: AtlasMCRecoFullPrecedenceDump.py:49
MonitoringHub.h
std::pair
Gaudi::Histogram2D::m_sumwy
double m_sumwy
Definition: H2D.h:57
std::vector< double >
ISvcLocator
Definition: ISvcLocator.h:46
jsonFromLHCbLog.json
json
Definition: jsonFromLHCbLog.py:87
P1D.h
Gaudi::slice1DY
std::pair< DataObject *, AIDA::IHistogram1D * > slice1DY(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D slice from 2D histogram.
ObjectFactory.h
Gaudi::Generic2D< AIDA::IHistogram2D, TH2D >::setTitle
bool setTitle(const std::string &title) override
Set the title of the object.
Definition: Generic2D.h:170
std::vector::front
T front(T... args)
Gaudi::Generic2D::cast
void * cast(const std::string &className) const override
Introspection method.
Utils.h
HistogramUtility.h
Gaudi::Generic2D::binEntries
int binEntries(int indexX, int indexY) const override
The number of entries (ie the number of times fill was called for this bin).
std::vector::push_back
T push_back(T... args)
bug_34121.t
t
Definition: bug_34121.py:30
Gaudi::svcLocator
GAUDI_API ISvcLocator * svcLocator()
Gaudi::Generic2D::adoptRepresentation
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
Gaudi::Histogram2D::setBinContents
virtual bool setBinContents(int binIndexX, int binIndexY, int entries, double height, double error, double centreX, double centreY)
Fast filling method for a given bin. It can be also the over/underflow bin.
Definition: H2D.cpp:159
HistoDumpEx.r
r
Definition: HistoDumpEx.py:20
ProduceConsume.j
j
Definition: ProduceConsume.py:101
std::string::c_str
T c_str(T... args)
AlgSequencer.h
h
Definition: AlgSequencer.py:32
std::array
STL class.
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:163
std::runtime_error
STL class.
H1D.h
Gaudi::Histogram2D::copyFromAida
void copyFromAida(const AIDA::IHistogram2D &h)
Create new histogram from any AIDA based histogram.
Definition: H2D.cpp:212
Gaudi
Header file for std:chrono::duration-based Counters.
Definition: __init__.py:1
GaudiPluginService.cpluginsvc.n
n
Definition: cpluginsvc.py:235
H2D.h
Gaudi::Histogram2D::Histogram2D
Histogram2D()
Standard Constructor.
Definition: H2D.cpp:146
Gaudi::Generic2D< AIDA::IHistogram2D, TH2D >::m_rep
std::unique_ptr< TH2D > m_rep
Reference to underlying implementation.
Definition: Generic2D.h:162
GaudiPI.h
Gaudi::Edges
std::vector< double > Edges
Definition: GaudiPI.h:28
Gaudi::slice1DX
std::pair< DataObject *, AIDA::IHistogram1D * > slice1DX(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D slice from 2D histogram.
Gaudi::project1DY
std::pair< DataObject *, AIDA::IHistogram1D * > project1DY(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D projection in Y from 2D histogram.
Gaudi::Generic2D< AIDA::IHistogram2D, TH2D >
Gaudi::Histogram2D::setRms
bool setRms(double rmsX, double rmsY)
Sets the rms of the histogram.
Definition: H2D.cpp:193
Gaudi::createH2D
std::pair< DataObject *, AIDA::IHistogram2D * > createH2D(ISvcLocator *svcLocator, const std::string &path, const AIDA::IHistogram2D &hist)
Copy constructor.
std::array::data
T data(T... args)
Gaudi::Histogram2D::fill
bool fill(double x, double y, double weight=1.) override
Fill the Histogram2D with a value and the.
Definition: H2D.cpp:186