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