laura is hosted by Hepforge, IPPP Durham
Laura++  3.6.0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
Lau2DCubicSpline.cc
Go to the documentation of this file.
1 
2 /*
3 Copyright 2013 University of Warwick
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 */
17 
18 /*
19 Laura++ package authors:
20 John Back
21 Paul Harrison
22 Thomas Latham
23 */
24 
32 #include "Lau2DCubicSpline.hh"
33 
34 #include <cmath>
35 #include <cstdlib>
36 #include <iostream>
37 
38 #include <TH2.h>
39 #include <TSystem.h>
40 
41 Lau2DCubicSpline::~Lau2DCubicSpline()
42 {
43 }
44 
45 inline Double_t Lau2DCubicSpline::histcont( const TH2& h, Int_t xbin, Int_t ybin ) const
46 {
47  //reflect until we're in range
48  while ( xbin < 0 || xbin >= nBinsX - 1 ) {
49  if ( xbin < 0 )
50  xbin = -xbin - 1;
51  if ( xbin >= nBinsX - 1 )
52  xbin = 2 * ( nBinsX - 1 ) - xbin - 1;
53  }
54  while ( ybin < 0 || ybin >= nBinsY - 1 ) {
55  if ( ybin < 0 )
56  ybin = -ybin - 1;
57  if ( ybin >= nBinsY - 1 )
58  ybin = 2 * ( nBinsY - 1 ) - ybin - 1;
59  }
60 
61  return h.GetBinContent( 1 + xbin, 1 + ybin );
62 }
63 
64 inline Double_t Lau2DCubicSpline::dhistdx( const TH2& h, Int_t xbin, Int_t ybin ) const
65 {
66  return 0.5 * ( histcont( h, xbin + 1, ybin ) - histcont( h, xbin - 1, ybin ) );
67 }
68 
69 inline Double_t Lau2DCubicSpline::dhistdy( const TH2& h, Int_t xbin, Int_t ybin ) const
70 {
71  return 0.5 * ( histcont( h, xbin, ybin + 1 ) - histcont( h, xbin, ybin - 1 ) );
72 }
73 
74 inline Double_t Lau2DCubicSpline::d2histdxdy( const TH2& h, Int_t xbin, Int_t ybin ) const
75 {
76  return 0.5 * ( histcont( h, xbin - 1, ybin - 1 ) - histcont( h, xbin + 1, ybin - 1 ) +
77  histcont( h, xbin + 1, ybin + 1 ) - histcont( h, xbin - 1, ybin + 1 ) );
78 }
79 
80 Lau2DCubicSpline::Lau2DCubicSpline( const TH2& h ) :
81  nBinsX( 1 + h.GetNbinsX() ),
82  nBinsY( 1 + h.GetNbinsY() ),
83  binSizeX( h.GetXaxis()->GetBinWidth( 1 ) ),
84  binSizeY( h.GetYaxis()->GetBinWidth( 1 ) ),
85  xmin( h.GetXaxis()->GetBinCenter( 1 ) - binSizeX ),
86  xmax( h.GetXaxis()->GetBinCenter( nBinsX - 1 ) + binSizeX ),
87  ymin( h.GetYaxis()->GetBinCenter( 1 ) - binSizeY ),
88  ymax( h.GetYaxis()->GetBinCenter( nBinsY - 1 ) + binSizeY ),
89  coeffs( CoeffRecLen * nBinsX * nBinsY )
90 {
91  const TAxis* xaxis = h.GetXaxis();
92  const TAxis* yaxis = h.GetYaxis();
93  // verify that all bins have same size
94  for ( Int_t i = 1; i < nBinsX; ++i ) {
95  if ( std::abs( xaxis->GetBinWidth( i ) / binSizeX - 1. ) > 1e-9 ) {
96  std::cerr << "ERROR in Lau2DCubicSpline constructor : the histogram has variable bin sizes."
97  << std::endl;
98  gSystem->Exit( EXIT_FAILURE );
99  }
100  }
101  for ( Int_t i = 1; i < nBinsY; ++i ) {
102  if ( std::abs( yaxis->GetBinWidth( i ) / binSizeY - 1. ) > 1e-9 ) {
103  std::cerr << "ERROR in Lau2DCubicSpline constructor : the histogram has variable bin sizes."
104  << std::endl;
105  gSystem->Exit( EXIT_FAILURE );
106  }
107  }
108  // ok, go through histogram to precalculate the interpolation coefficients
109  // in rectangles between bin centres
110  //
111  // for that purpose, we map each of those rectangles to the unit square
112  for ( Int_t j = -1; j < nBinsY - 1; ++j ) {
113  for ( Int_t i = -1; i < nBinsX - 1; ++i ) {
114  const Double_t rhs[NCoeff] = { // function values in bin centres
115  histcont( h, i, j ),
116  histcont( h, i + 1, j ),
117  histcont( h, i, j + 1 ),
118  histcont( h, i + 1, j + 1 ),
119  // df/dx in bin centres (finite difference approximation)
120  dhistdx( h, i, j ),
121  dhistdx( h, i + 1, j ),
122  dhistdx( h, i, j + 1 ),
123  dhistdx( h, i + 1, j + 1 ),
124  // df/dy in bin centres (finite difference approximation)
125  dhistdy( h, i, j ),
126  dhistdy( h, i + 1, j ),
127  dhistdy( h, i, j + 1 ),
128  dhistdy( h, i + 1, j + 1 ),
129  // d^2f/dxdy in bin centres (finite difference approximation)
130  d2histdxdy( h, i, j ),
131  d2histdxdy( h, i + 1, j ),
132  d2histdxdy( h, i, j + 1 ),
133  d2histdxdy( h, i + 1, j + 1 ) };
134  // work out solution - strange array placement is due to the fact
135  // that terms with x/y to high powers can be small, so they should
136  // be added up first during evaluation to avoid cancellation
137  // issues; at the same time you want to access them in order to
138  // not confuse the CPU cache, so they're stored back to front
139  //
140  // a_00 ... a_30
141  coeff( 1 + i, 1 + j, 15 ) = rhs[0];
142  coeff( 1 + i, 1 + j, 14 ) = rhs[4];
143  coeff( 1 + i, 1 + j, 13 ) = 3. * ( -rhs[0] + rhs[1] ) - 2. * rhs[4] - rhs[5];
144  coeff( 1 + i, 1 + j, 12 ) = 2. * ( rhs[0] - rhs[1] ) + rhs[4] + rhs[5];
145  // a_31 ... a_31
146  coeff( 1 + i, 1 + j, 11 ) = rhs[8];
147  coeff( 1 + i, 1 + j, 10 ) = rhs[12];
148  coeff( 1 + i, 1 + j, 9 ) = 3. * ( -rhs[8] + rhs[9] ) - 2. * rhs[12] - rhs[13];
149  coeff( 1 + i, 1 + j, 8 ) = 2. * ( rhs[8] - rhs[9] ) + rhs[12] + rhs[13];
150  // a_02 ... a_32
151  coeff( 1 + i, 1 + j, 7 ) = 3. * ( -rhs[0] + rhs[2] ) - 2. * rhs[8] - rhs[10];
152  coeff( 1 + i, 1 + j, 6 ) = 3. * ( -rhs[4] + rhs[6] ) - 2. * rhs[12] - rhs[14];
153  coeff( 1 + i, 1 + j, 5 ) = 9. * ( rhs[0] - rhs[1] - rhs[2] + rhs[3] ) +
154  6. * ( rhs[4] - rhs[6] + rhs[8] - rhs[9] ) + 4. * rhs[12] +
155  3. * ( rhs[5] - rhs[7] + rhs[10] - rhs[11] ) +
156  2. * ( rhs[13] + rhs[14] ) + rhs[15];
157  coeff( 1 + i, 1 + j, 4 ) = 6. * ( -rhs[0] + rhs[1] + rhs[2] - rhs[3] ) +
158  4. * ( -rhs[8] + rhs[9] ) +
159  3. * ( -rhs[4] - rhs[5] + rhs[6] + rhs[7] ) +
160  2. * ( -rhs[10] + rhs[11] - rhs[12] - rhs[13] ) - rhs[14] -
161  rhs[15];
162  // a_03 ... a_33
163  coeff( 1 + i, 1 + j, 3 ) = 2. * ( rhs[0] - rhs[2] ) + rhs[8] + rhs[10];
164  coeff( 1 + i, 1 + j, 2 ) = 2. * ( rhs[4] - rhs[6] ) + rhs[12] + rhs[14];
165  coeff( 1 + i, 1 + j, 1 ) = 6. * ( -rhs[0] + rhs[1] + rhs[2] - rhs[3] ) +
166  4. * ( -rhs[4] + rhs[6] ) +
167  3. * ( -rhs[8] + rhs[9] - rhs[10] + rhs[11] ) +
168  2. * ( -rhs[5] + rhs[7] - rhs[12] - rhs[14] ) - rhs[13] -
169  rhs[15];
170  coeff( 1 + i, 1 + j, 0 ) = 4. * ( rhs[0] - rhs[1] - rhs[2] + rhs[3] ) +
171  2. * ( rhs[4] + rhs[5] - rhs[6] - rhs[7] + rhs[8] - rhs[9] +
172  rhs[10] - rhs[11] ) +
173  rhs[12] + rhs[13] + rhs[14] + rhs[15];
174  // coeff(1 + i, 1 + j, 17) contains integral of function over the
175  // square in "unit square coordinates", i.e. neglecting bin widths
176  // this is done to help speed up calculations of 2D integrals
177  Double_t sum = 0.;
178  for ( Int_t k = 0; k < NCoeff; ++k )
179  sum += coeff( 1 + i, 1 + j, k ) / Double_t( ( 4 - ( k % 4 ) ) * ( 4 - ( k / 4 ) ) );
180  coeff( 1 + i, 1 + j, NCoeff ) = sum;
181  }
182  }
183 }
184 
185 Double_t Lau2DCubicSpline::evaluate( Double_t x, Double_t y ) const
186 {
187  // protect against NaN and out of range
188  if ( x <= xmin || x >= xmax || y <= ymin || y >= ymax || x != x || y != y )
189  return 0.;
190  // find the bin in question
191  const Int_t binx = Int_t( Double_t( nBinsX ) * ( x - xmin ) / ( xmax - xmin ) );
192  const Int_t biny = Int_t( Double_t( nBinsY ) * ( y - ymin ) / ( ymax - ymin ) );
193  // get low edge of bin
194  const Double_t xlo = Double_t( nBinsX - binx ) / Double_t( nBinsX ) * xmin +
195  Double_t( binx ) / Double_t( nBinsX ) * xmax;
196  const Double_t ylo = Double_t( nBinsY - biny ) / Double_t( nBinsY ) * ymin +
197  Double_t( biny ) / Double_t( nBinsY ) * ymax;
198  // normalise to coordinates in unit sqare
199  const Double_t hx = ( x - xlo ) / binSizeX;
200  const Double_t hy = ( y - ylo ) / binSizeY;
201  // monomials
202  const Double_t hxton[4] = { hx * hx * hx, hx * hx, hx, 1. };
203  const Double_t hyton[4] = { hy * hy * hy, hy * hy, hy, 1. };
204  // sum up
205  Double_t retVal = 0.;
206  for ( Int_t k = 0; k < NCoeff; ++k )
207  retVal += coeff( binx, biny, k ) * hxton[k % 4] * hyton[k / 4];
208 
209  return retVal;
210 }
211 
212 Double_t Lau2DCubicSpline::analyticalIntegral() const
213 {
214  return evalXY( xmin, xmax, ymin, ymax );
215 }
216 
217 Double_t Lau2DCubicSpline::analyticalIntegral( Double_t x1, Double_t x2, Double_t y1, Double_t y2 ) const
218 {
219  if ( y1 == y2 )
220  return evalX( x1, x2, y1 );
221  if ( x1 == x2 )
222  return evalY( x1, y1, y2 );
223  return evalXY( x1, x2, y1, y2 );
224 }
225 
226 Double_t Lau2DCubicSpline::evalX( Double_t x1, Double_t x2, Double_t y ) const
227 {
228  // protect against NaN
229  if ( x1 != x1 || x2 != x2 || y != y )
230  return 0.;
231  // find the bin in question
232  const Int_t biny = Int_t( Double_t( nBinsY ) * ( y - ymin ) / ( ymax - ymin ) );
233  // get low edge of bin
234  const Double_t ylo = Double_t( nBinsY - biny ) / Double_t( nBinsY ) * ymin +
235  Double_t( biny ) / Double_t( nBinsY ) * ymax;
236  // normalise to coordinates in unit sqare
237  const Double_t hy = ( y - ylo ) / binSizeY;
238  // monomials
239  const Double_t hyton[4] = { hy * hy * hy, hy * hy, hy, 1. };
240  // integral
241  Double_t sum = 0.;
242  for ( Int_t binx = 0; binx < nBinsX; ++binx ) {
243  // get low/high edge of bin
244  const Double_t xhi = Double_t( nBinsX - binx - 1 ) / Double_t( nBinsX ) * xmin +
245  Double_t( binx + 1 ) / Double_t( nBinsX ) * xmax;
246  if ( xhi < x1 )
247  continue;
248  const Double_t xlo = Double_t( nBinsX - binx ) / Double_t( nBinsX ) * xmin +
249  Double_t( binx ) / Double_t( nBinsX ) * xmax;
250  if ( xlo > x2 )
251  break;
252  // work out integration range
253  const Double_t a = ( ( xlo > x1 ) ? 0. : ( x1 - xlo ) ) / binSizeX;
254  const Double_t b = ( ( xhi < x2 ) ? binSizeX : ( x2 - xlo ) ) / binSizeX;
255  // integrated monomials
256  const Double_t hxton[4] = { 0.25 * ( b * b * b * b - a * a * a * a ),
257  ( b * b * b - a * a * a ) / 3.,
258  0.5 * ( b * b - a * a ),
259  b - a };
260  Double_t lsum = 0.;
261  for ( Int_t k = 0; k < NCoeff; ++k )
262  lsum += coeff( binx, biny, k ) * hxton[k % 4] * hyton[k / 4];
263  sum += lsum;
264  }
265  // move from unit square coordinates to user coordinates
266  return sum * binSizeX;
267 }
268 
269 Double_t Lau2DCubicSpline::evalY( Double_t x, Double_t y1, Double_t y2 ) const
270 {
271  // protect against NaN
272  if ( x != x || y1 != y1 || y2 != y2 )
273  return 0.;
274  // find the bin in question
275  const Int_t binx = Int_t( Double_t( nBinsX ) * ( x - xmin ) / ( xmax - xmin ) );
276  // get low edge of bin
277  const Double_t xlo = Double_t( nBinsX - binx ) / Double_t( nBinsX ) * xmin +
278  Double_t( binx ) / Double_t( nBinsX ) * xmax;
279  // normalise to coordinates in unit sqare
280  const Double_t hx = ( x - xlo ) / binSizeX;
281  // monomials
282  const Double_t hxton[4] = { hx * hx * hx, hx * hx, hx, 1. };
283  // integral
284  Double_t sum = 0.;
285  for ( Int_t biny = 0; biny < nBinsY; ++biny ) {
286  // get low/high edge of bin
287  const Double_t yhi = Double_t( nBinsY - biny - 1 ) / Double_t( nBinsY ) * ymin +
288  Double_t( biny + 1 ) / Double_t( nBinsY ) * ymax;
289  if ( yhi < y1 )
290  continue;
291  const Double_t ylo = Double_t( nBinsY - biny ) / Double_t( nBinsY ) * ymin +
292  Double_t( biny ) / Double_t( nBinsY ) * ymax;
293  if ( ylo > y2 )
294  break;
295  // work out integration range
296  const Double_t a = ( ( ylo > y1 ) ? 0. : ( y1 - ylo ) ) / binSizeY;
297  const Double_t b = ( ( yhi < y2 ) ? binSizeY : ( y2 - ylo ) ) / binSizeY;
298  // integrated monomials
299  const Double_t hyton[4] = { 0.25 * ( b * b * b * b - a * a * a * a ),
300  ( b * b * b - a * a * a ) / 3.,
301  0.5 * ( b * b - a * a ),
302  b - a };
303  Double_t lsum = 0.;
304  for ( Int_t k = 0; k < NCoeff; ++k )
305  lsum += coeff( binx, biny, k ) * hxton[k % 4] * hyton[k / 4];
306  sum += lsum;
307  }
308  // move from unit square coordinates to user coordinates
309  return sum * binSizeY;
310 }
311 
312 Double_t Lau2DCubicSpline::evalXY( Double_t x1, Double_t x2, Double_t y1, Double_t y2 ) const
313 {
314  // protect against NaN
315  if ( x1 != x1 || x2 != x2 || y1 != y1 || y2 != y2 )
316  return 0.;
317  // integral
318  Double_t sum = 0.;
319  for ( Int_t biny = 0; biny < nBinsY; ++biny ) {
320  // get low/high edge of bin
321  const Double_t yhi = Double_t( nBinsY - biny - 1 ) / Double_t( nBinsY ) * ymin +
322  Double_t( biny + 1 ) / Double_t( nBinsY ) * ymax;
323  if ( yhi < y1 )
324  continue;
325  const Double_t ylo = Double_t( nBinsY - biny ) / Double_t( nBinsY ) * ymin +
326  Double_t( biny ) / Double_t( nBinsY ) * ymax;
327  if ( ylo > y2 )
328  break;
329  // work out integration range
330  const Double_t ay = ( ( ylo > y1 ) ? 0. : ( y1 - ylo ) ) / binSizeY;
331  const Double_t by = ( ( yhi < y2 ) ? binSizeY : ( y2 - ylo ) ) / binSizeY;
332  const Bool_t yFullyContained = std::abs( by - ay - 1.0 ) < 1e-15;
333  // integrated monomials
334  const Double_t hyton[4] = { 0.25 * ( by * by * by * by - ay * ay * ay * ay ),
335  ( by * by * by - ay * ay * ay ) / 3.,
336  0.5 * ( by * by - ay * ay ),
337  by - ay };
338  for ( Int_t binx = 0; binx < nBinsX; ++binx ) {
339  // get low/high edge of bin
340  const Double_t xhi = Double_t( nBinsX - binx - 1 ) / Double_t( nBinsX ) * xmin +
341  Double_t( binx + 1 ) / Double_t( nBinsX ) * xmax;
342  if ( xhi < x1 )
343  continue;
344  const Double_t xlo = Double_t( nBinsX - binx ) / Double_t( nBinsX ) * xmin +
345  Double_t( binx ) / Double_t( nBinsX ) * xmax;
346  if ( xlo > x2 )
347  break;
348  // work out integration range
349  const Double_t ax = ( ( xlo > x1 ) ? 0. : ( x1 - xlo ) ) / binSizeX;
350  const Double_t bx = ( ( xhi < x2 ) ? binSizeX : ( x2 - xlo ) ) / binSizeX;
351  const Bool_t xFullyContained = std::abs( bx - ax - 1.0 ) < 1e-15;
352  if ( xFullyContained && yFullyContained ) {
353  // for fully contained bins, we have cached the integral
354  sum += coeff( binx, biny, NCoeff );
355  continue;
356  }
357  // integrated monomials
358  const Double_t hxton[4] = { 0.25 * ( bx * bx * bx * bx - ax * ax * ax * ax ),
359  ( bx * bx * bx - ax * ax * ax ) / 3.,
360  0.5 * ( bx * bx - ax * ax ),
361  bx - ax };
362  // integrate over bin
363 
364  Double_t lsum = 0.;
365  for ( Int_t k = 0; k < NCoeff; ++k )
366  lsum += coeff( binx, biny, k ) * hxton[k % 4] * hyton[k / 4];
367  sum += lsum;
368  }
369  }
370  // move from unit square coordinates to user coordinates
371  return sum * binSizeX * binSizeY;
372 }
Lau2DCubicSpline::binSizeY
Double_t binSizeY
Bin size in y.
Definition: Lau2DCubicSpline.hh:119
Lau2DCubicSpline::dhistdx
Double_t dhistdx(const TH2 &h, Int_t xbin, Int_t ybin) const
Get d/dx finite difference in a given bin from a histogram.
Definition: Lau2DCubicSpline.cc:64
Lau2DCubicSpline::nBinsX
Int_t nBinsX
Number of bins in x axis.
Definition: Lau2DCubicSpline.hh:110
Lau2DCubicSpline::evaluate
virtual Double_t evaluate(Double_t x, Double_t y) const
Evaluate the function at given point.
Definition: Lau2DCubicSpline.cc:185
Lau2DCubicSpline::xmax
Double_t xmax
Maximum x value.
Definition: Lau2DCubicSpline.hh:125
Lau2DCubicSpline::nBinsY
Int_t nBinsY
Number of bins in y axis.
Definition: Lau2DCubicSpline.hh:113
Lau2DCubicSpline::evalY
Double_t evalY(Double_t x, Double_t y1, Double_t y2) const
Evaluate integral over y at a given x from (x, y1) to (x, y2)
Definition: Lau2DCubicSpline.cc:269
Lau2DCubicSpline::~Lau2DCubicSpline
virtual ~Lau2DCubicSpline()
Destructor.
Definition: Lau2DCubicSpline.cc:41
Lau2DCubicSpline::analyticalIntegral
virtual Double_t analyticalIntegral() const
Evaluate analytical integral across the whole function.
Definition: Lau2DCubicSpline.cc:212
Lau2DCubicSpline::Lau2DCubicSpline
Lau2DCubicSpline(const TH2 &h)
Constructor from histogram.
Definition: Lau2DCubicSpline.cc:80
Lau2DCubicSpline::evalX
Double_t evalX(Double_t x1, Double_t x2, Double_t y) const
Evaluate integral over x at a given y from (x1, y) to (x2, y)
Definition: Lau2DCubicSpline.cc:226
Lau2DCubicSpline::binSizeX
Double_t binSizeX
Bin size in x.
Definition: Lau2DCubicSpline.hh:116
Lau2DCubicSpline::coeff
const Double_t & coeff(Int_t binx, Int_t biny, Int_t theCoeff) const
Const access to coefficients.
Definition: Lau2DCubicSpline.hh:179
Lau2DCubicSpline::ymax
Double_t ymax
Maximum y value.
Definition: Lau2DCubicSpline.hh:131
Lau2DCubicSpline::histcont
Double_t histcont(const TH2 &h, Int_t xbin, Int_t ybin) const
Get contents of a given bin from a histogram.
Definition: Lau2DCubicSpline.cc:45
Lau2DCubicSpline::dhistdy
Double_t dhistdy(const TH2 &h, Int_t xbin, Int_t ybin) const
Get d/dy finite difference in a given bin from a histogram.
Definition: Lau2DCubicSpline.cc:69
Lau2DCubicSpline::d2histdxdy
Double_t d2histdxdy(const TH2 &h, Int_t xbin, Int_t ybin) const
Get d^2/dydx finite difference in a given bin from a histogram.
Definition: Lau2DCubicSpline.cc:74
Lau2DCubicSpline.hh
File containing declaration of Lau2DCubicSpline class.
Lau2DCubicSpline::xmin
Double_t xmin
Minimum x value.
Definition: Lau2DCubicSpline.hh:122
Lau2DCubicSpline::ymin
Double_t ymin
Minimum y value.
Definition: Lau2DCubicSpline.hh:128
Lau2DCubicSpline::evalXY
Double_t evalXY(Double_t x1, Double_t x2, Double_t y1, Double_t y2) const
Evaluate integral over x and y from (x1, y1) to (x2, y2)
Definition: Lau2DCubicSpline.cc:312
Generated by   doxygen 1.8.18