Loading [MathJax]/extensions/tex2jax.js
The Gaudi Framework  v31r0 (aeb156f0)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
H2D.cpp
Go to the documentation of this file.
1 #ifdef __ICC
2 // disable icc remark #2259: non-pointer conversion from "X" to "Y" may lose significant bits
3 // TODO: To be removed, since it comes from ROOT TMathBase.h
4 # pragma warning( disable : 2259 )
5 #endif
6 #ifdef WIN32
7 // Disable warning
8 // warning C4996: 'sprintf': This function or variable may be unsafe.
9 // coming from TString.h
10 # pragma warning( disable : 4996 )
11 #endif
12 #include "GaudiPI.h"
13 #include <GaudiCommonSvc/H1D.h>
14 #include <GaudiCommonSvc/H2D.h>
16 #include <GaudiCommonSvc/P1D.h>
18 #include <TH2D.h>
19 #include <TProfile.h>
20 #include <array>
21 
22 namespace {
23  using AIDA::IHistogram1D;
24  using AIDA::IHistogram2D;
25  using AIDA::IProfile1D;
26 } // namespace
27 
28 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( const std::string& title, int binsX, double iminX, double imaxX,
29  int binsY, double iminY, double imaxY ) {
30  auto p = new Histogram2D( new TH2D( title.c_str(), title.c_str(), binsX, iminX, imaxX, binsY, iminY, imaxY ) );
31  return {p, p};
32 }
33 
35  auto p = new Histogram2D(
36  new TH2D( title.c_str(), title.c_str(), eX.size() - 1, &eX.front(), eY.size() - 1, &eY.front() ) );
37  return {p, p};
38 }
39 
41  auto p = new Histogram2D( rep );
42  return {p, p};
43 }
44 
45 std::pair<DataObject*, IHistogram2D*> Gaudi::createH2D( const IHistogram2D& hist ) {
46  TH2D* h = getRepresentation<AIDA::IHistogram2D, TH2D>( hist );
47  Histogram2D* n = h ? new Histogram2D( new TH2D( *h ) ) : nullptr;
48  return {n, n};
49 }
50 
51 std::pair<DataObject*, IHistogram1D*> Gaudi::slice1DX( const std::string& nam, const IHistogram2D& hist, int first,
52  int last ) {
53  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
54  TH1D* t = r ? r->ProjectionX( "_px", first, last, "e" ) : nullptr;
55  if ( t ) t->SetName( nam.c_str() );
56  Histogram1D* p = ( t ? new Histogram1D( t ) : nullptr );
57  return {p, p};
58 }
59 
60 std::pair<DataObject*, IHistogram1D*> Gaudi::slice1DY( const std::string& nam, const IHistogram2D& hist, int first,
61  int last ) {
62  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
63  TH1D* t = r ? r->ProjectionY( "_py", first, last, "e" ) : nullptr;
64  if ( t ) t->SetName( nam.c_str() );
65  Histogram1D* p = ( t ? new Histogram1D( t ) : nullptr );
66  return {p, p};
67 }
68 
69 std::pair<DataObject*, IHistogram1D*> Gaudi::project1DY( const std::string& nam, const IHistogram2D& hist, int first,
70  int last ) {
71  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
72  TH1D* t = r ? r->ProjectionY( "_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*, IProfile1D*> Gaudi::profile1DX( const std::string& nam, const IHistogram2D& hist, int first,
79  int last ) {
80  TH2* r = Gaudi::getRepresentation<IHistogram2D, TH2>( hist );
81  TProfile* t = r ? r->ProfileX( "_pfx", first, last, "e" ) : nullptr;
82  if ( t ) t->SetName( nam.c_str() );
83  Profile1D* p = ( t ? new Profile1D( t ) : nullptr );
84  return {p, p};
85 }
86 
87 std::pair<DataObject*, IProfile1D*> Gaudi::profile1DY( const std::string& nam, const IHistogram2D& hist, int first,
88  int last ) {
89  TH2* r = getRepresentation<IHistogram2D, TH2>( hist );
90  TProfile* t = r ? r->ProfileY( "_pfx", first, last, "e" ) : nullptr;
91  if ( t ) t->SetName( nam.c_str() );
92  Profile1D* p = ( t ? new Profile1D( t ) : nullptr );
93  return {p, p};
94 }
95 
96 namespace Gaudi {
97  template <>
98  void* Generic2D<IHistogram2D, TH2D>::cast( const std::string& className ) const {
99  if ( className == "AIDA::IHistogram2D" )
100  return (IHistogram2D*)this;
101  else if ( className == "AIDA::IHistogram" )
102  return (IHistogram*)this;
103  return nullptr;
104  }
105 
106  template <>
107  int Generic2D<IHistogram2D, TH2D>::binEntries( int indexX, int indexY ) const {
108  if ( binHeight( indexX, indexY ) <= 0 ) return 0;
109  double xx = binHeight( indexX, indexY ) / binError( indexX, indexY );
110  return int( xx * xx + 0.5 );
111  }
112 
113  template <>
115  TH2D* imp = dynamic_cast<TH2D*>( rep );
116  if ( !imp ) throw std::runtime_error( "Cannot adopt native histogram representation." );
117  m_rep.reset( imp );
118  m_xAxis.initialize( m_rep->GetXaxis(), true );
119  m_yAxis.initialize( m_rep->GetYaxis(), true );
120  const TArrayD* a = m_rep->GetSumw2();
121  if ( !a || ( a && a->GetSize() == 0 ) ) m_rep->Sumw2();
122  setTitle( m_rep->GetTitle() );
123  }
124 } // namespace Gaudi
125 
127  m_rep->Sumw2();
128  m_sumwx = m_sumwy = 0;
129  setTitle( "" );
130  m_rep->SetDirectory( nullptr );
131 }
132 
134  adoptRepresentation( rep );
135  m_sumwx = m_sumwy = 0;
136  m_rep->SetDirectory( nullptr );
137 }
138 
139 bool Gaudi::Histogram2D::setBinContents( int i, int j, int entries, double height, double error, double centreX,
140  double centreY ) {
141  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), height );
142  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), error );
143  // accumulate sumwx for in range bins
144  if ( i >= 0 && j >= 0 ) {
145  m_sumwx += centreX * height;
146  m_sumwy += centreY * height;
147  }
149  return true;
150 }
151 
153  m_sumwx = 0;
154  m_sumwy = 0;
155  return Base::reset();
156 }
157 
158 #ifdef __ICC
159 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
160 // The comparison are meant
161 # pragma warning( push )
162 # pragma warning( disable : 1572 )
163 #endif
164 bool Gaudi::Histogram2D::fill( double x, double y, double weight ) {
165  // avoid race conditiosn when filling the histogram
167  ( weight == 1. ) ? m_rep->Fill( x, y ) : m_rep->Fill( x, y, weight );
168  return true;
169 }
170 
171 bool Gaudi::Histogram2D::setRms( double rmsX, double rmsY ) {
172  m_rep->SetEntries( m_sumEntries );
173  std::vector<double> stat( 11 );
174  stat[0] = sumBinHeights();
175  stat[1] = 0;
176  if ( equivalentBinEntries() != 0 ) stat[1] = ( sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
177  stat[2] = m_sumwx;
178  double meanX = 0;
179  if ( sumBinHeights() != 0 ) meanX = m_sumwx / sumBinHeights();
180  stat[3] = ( meanX * meanX + rmsX * rmsX ) * sumBinHeights();
181  stat[4] = m_sumwy;
182  double meanY = 0;
183  if ( sumBinHeights() != 0 ) meanY = m_sumwy / sumBinHeights();
184  stat[5] = ( meanY * meanY + rmsY * rmsY ) * sumBinHeights();
185  stat[6] = 0;
186  m_rep->PutStats( &stat.front() );
187  return true;
188 }
189 
190 void Gaudi::Histogram2D::copyFromAida( const IHistogram2D& h ) {
191  // implement here the copy
192  const char* title = h.title().c_str();
193  if ( h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning() )
194  m_rep.reset( new TH2D( title, title, h.xAxis().bins(), h.xAxis().lowerEdge(), h.xAxis().upperEdge(),
195  h.yAxis().bins(), h.yAxis().lowerEdge(), h.yAxis().upperEdge() ) );
196  else {
197  Edges eX, eY;
198  for ( int i = 0; i < h.xAxis().bins(); ++i ) eX.push_back( h.xAxis().binLowerEdge( i ) );
199  // add also upperedges at the end
200  eX.push_back( h.xAxis().upperEdge() );
201  for ( int i = 0; i < h.yAxis().bins(); ++i ) eY.push_back( h.yAxis().binLowerEdge( i ) );
202  // add also upperedges at the end
203  eY.push_back( h.yAxis().upperEdge() );
204  m_rep.reset( new TH2D( title, title, eX.size() - 1, &eX.front(), eY.size() - 1, &eY.front() ) );
205  }
206  m_xAxis.initialize( m_rep->GetXaxis(), true );
207  m_yAxis.initialize( m_rep->GetYaxis(), true );
208  m_rep->Sumw2();
209  m_sumEntries = 0;
210  m_sumwx = 0;
211  m_sumwy = 0;
212  // statistics
213  double sumw = h.sumBinHeights();
214  double sumw2 = 0;
215  if ( h.equivalentBinEntries() != 0 ) sumw2 = ( sumw * sumw ) / h.equivalentBinEntries();
216  double sumwx = h.meanX() * h.sumBinHeights();
217  double sumwx2 = ( h.meanX() * h.meanX() + h.rmsX() * h.rmsX() ) * h.sumBinHeights();
218  double sumwy = h.meanY() * h.sumBinHeights();
219  double sumwy2 = ( h.meanY() * h.meanY() + h.rmsY() * h.rmsY() ) * h.sumBinHeights();
220  double sumwxy = 0;
221 
222  // copy the contents in (AIDA underflow/overflow are -2,-1)
223  for ( int i = -2; i < xAxis().bins(); ++i ) {
224  for ( int j = -2; j < yAxis().bins(); ++j ) {
225  // root binning starts from one !
226  m_rep->SetBinContent( rIndexX( i ), rIndexY( j ), h.binHeight( i, j ) );
227  m_rep->SetBinError( rIndexX( i ), rIndexY( j ), h.binError( i, j ) );
228  // calculate statistics
229  if ( i >= 0 && j >= 0 ) { sumwxy += h.binHeight( i, j ) * h.binMeanX( i, j ) * h.binMeanY( i, j ); }
230  }
231  }
232  // need to do set entries after setting contents otherwise root will recalculate them
233  // taking into account how many time SetBinContents() has been called
234  m_rep->SetEntries( h.allEntries() );
235  // fill stat vector
236  std::array<double, 11> stat = {{sumw, sumw2, sumwx, sumwx2, sumwy, sumwy2, sumwxy}};
237  m_rep->PutStats( stat.data() );
238 }
239 
240 #ifdef __ICC
241 // re-enable icc remark #1572
242 # pragma warning( pop )
243 #endif
int m_sumEntries
cache sumEntries (allEntries) when setting contents since Root can&#39;t compute by himself ...
Definition: Generic2D.h:155
bool setRms(double rmsX, double rmsY)
Sets the rms of the histogram.
Definition: H2D.cpp:171
virtual double equivalentBinEntries() const
Number of equivalent entries, i.e. SUM[ weight ] ^ 2 / SUM[ weight^2 ]
Definition: Generic2D.h:303
std::pair< DataObject *, AIDA::IHistogram2D * > createH2D(const AIDA::IHistogram2D &hist)
Copy constructor.
double rmsX() const override
Returns the rms of the profile as calculated on filling-time projected on the X axis.
Definition: Generic2D.h:267
T front(T...args)
std::pair< DataObject *, AIDA::IHistogram1D * > slice1DY(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D slice from 2D histogram.
bool fill(double x, double y, double weight=1.) override
Fill the Histogram2D with a value and the.
Definition: H2D.cpp:164
bool reset() override
Definition: H2D.cpp:152
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:139
std::vector< double > Edges
Definition: GaudiPI.h:17
double rmsY() const override
Returns the rms of the profile as calculated on filling-time projected on the Y axis.
Definition: Generic2D.h:272
STL class.
void adoptRepresentation(TObject *rep) override
Adopt ROOT histogram representation.
virtual int rIndexX(int index) const
operator methods
Definition: Generic2D.h:68
Histogram2D()
Standard Constructor.
Definition: H2D.cpp:126
double sumBinHeights() const override
Get the sum of in range bin heights in the IProfile.
Definition: Generic2D.h:194
T push_back(T...args)
T data(T...args)
Axis m_xAxis
X axis member.
Definition: Generic2D.h:145
int binEntries(int indexX, int indexY) const override
The number of entries (ie the number of times fill was called for this bin).
T reset(T...args)
std::pair< DataObject *, AIDA::IHistogram1D * > slice1DX(const std::string &name, const AIDA::IHistogram2D &h, int firstbin, int lastbin)
Create 1D slice from 2D histogram.
void initialize(TAxis *itaxi, bool)
Definition: Axis.h:63
bool setTitle(const std::string &title) override
Set the title of the object.
Definition: Generic2D.h:159
double meanX() const override
Returns the mean of the profile, as calculated on filling-time projected on the X axis...
Definition: Generic2D.h:257
std::mutex m_fillSerialization
Definition: H2D.h:46
T c_str(T...args)
Axis m_yAxis
Y axis member.
Definition: Generic2D.h:147
void * cast(const std::string &className) const override
Introspection method.
STL class.
std::unique_ptr< IMPLEMENTATION > m_rep
Reference to underlying implementation.
Definition: Generic2D.h:151
int entries() const override
Get the number or all the entries.
Definition: Generic2D.h:174
void copyFromAida(const AIDA::IHistogram2D &h)
Create new histogram from any AIDA based histogram.
Definition: H2D.cpp:190
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.
virtual int rIndexY(int index) const
operator methods
Definition: Generic2D.h:70
double m_sumwy
Definition: H2D.h:43
double m_sumwx
Definition: H2D.h:42
double meanY() const override
Returns the mean of the profile, as calculated on filling-time projected on the Y axis...
Definition: Generic2D.h:262
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.
Helper functions to set/get the application return code.
Definition: __init__.py:1
const AIDA::IAxis & yAxis() const override
Return the Y axis.
Definition: Generic2D.h:66
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.
const AIDA::IAxis & xAxis() const override
Return the X axis.
Definition: Generic2D.h:64
Common AIDA implementation stuff for histograms and profiles using ROOT implementations.
Definition: Generic2D.h:35