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