laura is hosted by Hepforge, IPPP Durham
Laura++  v3r1
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 
28 { }
29 
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 
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 
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  const TAxis* xaxis = h.GetXaxis();
80  const TAxis* yaxis = h.GetYaxis();
81  // verify that all bins have same size
82  for (Int_t i = 1; i < nBinsX; ++i) {
83  if (std::abs(xaxis->GetBinWidth(i) / binSizeX - 1.) > 1e-9) {
84  std::cerr << "ERROR in Lau2DCubicSpline constructor : the histogram has variable bin sizes." << std::endl;
85  gSystem->Exit(EXIT_FAILURE);
86  }
87  }
88  for (Int_t i = 1; i < nBinsY; ++i) {
89  if (std::abs(yaxis->GetBinWidth(i) / binSizeY - 1.) > 1e-9) {
90  std::cerr << "ERROR in Lau2DCubicSpline constructor : the histogram has variable bin sizes." << std::endl;
91  gSystem->Exit(EXIT_FAILURE);
92  }
93  }
94  // ok, go through histogram to precalculate the interpolation coefficients
95  // in rectangles between bin centres
96  //
97  // for that purpose, we map each of those rectangles to the unit square
98  for (Int_t j = -1; j < nBinsY - 1; ++j) {
99  for (Int_t i = -1; i < nBinsX - 1; ++i) {
100  const Double_t rhs[NCoeff] = {
101  // function values in bin centres
102  histcont(h, i, j),
103  histcont(h, i + 1, j),
104  histcont(h, i, j + 1),
105  histcont(h, i + 1, j + 1),
106  // df/dx in bin centres (finite difference approximation)
107  dhistdx(h, i, j),
108  dhistdx(h, i + 1, j),
109  dhistdx(h, i, j + 1),
110  dhistdx(h, i + 1, j + 1),
111  // df/dy in bin centres (finite difference approximation)
112  dhistdy(h, i, j),
113  dhistdy(h, i + 1, j),
114  dhistdy(h, i, j + 1),
115  dhistdy(h, i + 1, j + 1),
116  // d^2f/dxdy in bin centres (finite difference approximation)
117  d2histdxdy(h, i, j),
118  d2histdxdy(h, i + 1, j),
119  d2histdxdy(h, i, j + 1),
120  d2histdxdy(h, i + 1, j + 1)
121  };
122  // work out solution - strange array placement is due to the fact
123  // that terms with x/y to high powers can be small, so they should
124  // be added up first during evaluation to avoid cancellation
125  // issues; at the same time you want to access them in order to
126  // not confuse the CPU cache, so they're stored back to front
127  //
128  // a_00 ... a_30
129  coeff(1 + i, 1 + j, 15) = rhs[0];
130  coeff(1 + i, 1 + j, 14) = rhs[4];
131  coeff(1 + i, 1 + j, 13) =
132  3. * (-rhs[0] + rhs[1]) - 2. * rhs[4] - rhs[5];
133  coeff(1 + i, 1 + j, 12) =
134  2. * (rhs[0] - rhs[1]) + rhs[4] + rhs[5];
135  // a_31 ... a_31
136  coeff(1 + i, 1 + j, 11) = rhs[8];
137  coeff(1 + i, 1 + j, 10) = rhs[12];
138  coeff(1 + i, 1 + j, 9) =
139  3. * (-rhs[8] + rhs[9]) - 2. * rhs[12] - rhs[13];
140  coeff(1 + i, 1 + j, 8) =
141  2. * (rhs[8] - rhs[9]) + rhs[12] + rhs[13];
142  // a_02 ... a_32
143  coeff(1 + i, 1 + j, 7) =
144  3. * (-rhs[0] + rhs[2]) - 2. * rhs[8] - rhs[10];
145  coeff(1 + i, 1 + j, 6) =
146  3. * (-rhs[4] + rhs[6]) - 2. * rhs[12] - rhs[14];
147  coeff(1 + i, 1 + j, 5) =
148  9. * (rhs[0] - rhs[1] - rhs[2] + rhs[3]) +
149  6. * (rhs[4] - rhs[6] + rhs[8] - rhs[9]) + 4. * rhs[12] +
150  3. * (rhs[5] - rhs[7] + rhs[10] - rhs[11]) +
151  2. * (rhs[13] + rhs[14]) + rhs[15];
152  coeff(1 + i, 1 + j, 4) =
153  6. * (-rhs[0] + rhs[1] + rhs[2] - rhs[3]) +
154  4. * (-rhs[8] + rhs[9]) +
155  3. * (-rhs[4] - rhs[5] + rhs[6] + rhs[7]) +
156  2. * (-rhs[10] + rhs[11] - rhs[12] - rhs[13]) -
157  rhs[14] - rhs[15];
158  // a_03 ... a_33
159  coeff(1 + i, 1 + j, 3) =
160  2. * (rhs[0] - rhs[2]) + rhs[8] + rhs[10];
161  coeff(1 + i, 1 + j, 2) =
162  2. * (rhs[4] - rhs[6]) + rhs[12] + rhs[14];
163  coeff(1 + i, 1 + j, 1) =
164  6. * (-rhs[0] + rhs[1] + rhs[2] - rhs[3]) +
165  4. * (-rhs[4] + rhs[6]) +
166  3. * (-rhs[8] + rhs[9] - rhs[10] + rhs[11]) +
167  2. * (- rhs[5] + rhs[7] - rhs[12] - rhs[14]) -
168  rhs[13] - rhs[15];
169  coeff(1 + i, 1 + j, 0) =
170  4. * (rhs[0] - rhs[1] - rhs[2] + rhs[3]) +
171  2. * (rhs[4] + rhs[5] - rhs[6] - rhs[7] +
172  rhs[8] - rhs[9] + 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) /
180  Double_t((4 - (k % 4)) * (4 - (k / 4)));
181  coeff(1 + i, 1 + j, NCoeff) = sum;
182  }
183  }
184 }
185 
186 Double_t Lau2DCubicSpline::evaluate(Double_t x, Double_t y) const
187 {
188  // protect against NaN and out of range
189  if (x <= xmin || x >= xmax || y <= ymin || y >= ymax || x != x || y != y)
190  return 0.;
191  // find the bin in question
192  const Int_t binx = Int_t(Double_t(nBinsX) * (x - xmin) / (xmax - xmin));
193  const Int_t biny = Int_t(Double_t(nBinsY) * (y - ymin) / (ymax - ymin));
194  // get low edge of bin
195  const Double_t xlo = Double_t(nBinsX - binx) / Double_t(nBinsX) * xmin +
196  Double_t(binx) / Double_t(nBinsX) * xmax;
197  const Double_t ylo = Double_t(nBinsY - biny) / Double_t(nBinsY) * ymin +
198  Double_t(biny) / Double_t(nBinsY) * ymax;
199  // normalise to coordinates in unit sqare
200  const Double_t hx = (x - xlo) / binSizeX;
201  const Double_t hy = (y - ylo) / binSizeY;
202  // monomials
203  const Double_t hxton[4] = { hx * hx * hx, hx * hx, hx, 1. };
204  const Double_t hyton[4] = { hy * hy * hy, hy * hy, hy, 1. };
205  // sum up
206  Double_t retVal = 0.;
207  for (Int_t k = 0; k < NCoeff; ++k)
208  retVal += coeff(binx, biny, k) * hxton[k % 4] * hyton[k / 4];
209 
210  return retVal;
211 }
212 
214 {
215  return evalXY(xmin,xmax,ymin,ymax);
216 }
217 
218 Double_t Lau2DCubicSpline::analyticalIntegral(Double_t x1, Double_t x2, Double_t y1, Double_t y2) const
219 {
220  if(y1==y2) return evalX(x1, x2, y1);
221  if(x1==x2) return evalY(x1, y1, y2);
222  return evalXY(x1, x2, y1, y2);
223 }
224 
225 Double_t Lau2DCubicSpline::evalX(Double_t x1, Double_t x2, Double_t y) const
226 {
227  // protect against NaN
228  if (x1 != x1 || x2 != x2 || y != y) return 0.;
229  // find the bin in question
230  const Int_t biny = Int_t(Double_t(nBinsY) * (y - ymin) / (ymax - ymin));
231  // get low edge of bin
232  const Double_t ylo = Double_t(nBinsY - biny) / Double_t(nBinsY) * ymin +
233  Double_t(biny) / Double_t(nBinsY) * ymax;
234  // normalise to coordinates in unit sqare
235  const Double_t hy = (y - ylo) / binSizeY;
236  // monomials
237  const Double_t hyton[4] = { hy * hy * hy, hy * hy, hy, 1. };
238  // integral
239  Double_t sum = 0.;
240  for (Int_t binx = 0; binx < nBinsX; ++binx) {
241  // get low/high edge of bin
242  const Double_t xhi = Double_t(nBinsX - binx - 1) / Double_t(nBinsX) * xmin +
243  Double_t(binx + 1) / Double_t(nBinsX) * xmax;
244  if (xhi < x1) continue;
245  const Double_t xlo = Double_t(nBinsX - binx) / Double_t(nBinsX) * xmin +
246  Double_t(binx) / Double_t(nBinsX) * xmax;
247  if (xlo > x2) break;
248  // work out integration range
249  const Double_t a = ((xlo > x1) ? 0. : (x1 - xlo)) / binSizeX;
250  const Double_t b = ((xhi < x2) ? binSizeX : (x2 - xlo)) / binSizeX;
251  // integrated monomials
252  const Double_t hxton[4] = { 0.25 * (b * b * b * b - a * a * a * a),
253  (b * b * b - a * a * a) / 3., 0.5 * (b * b - a * a), b - a };
254  Double_t lsum = 0.;
255  for (Int_t k = 0; k < NCoeff; ++k)
256  lsum += coeff(binx, biny, k) * hxton[k % 4] * hyton[k / 4];
257  sum += lsum;
258  }
259  // move from unit square coordinates to user coordinates
260  return sum * binSizeX;
261 }
262 
263 Double_t Lau2DCubicSpline::evalY(Double_t x, Double_t y1, Double_t y2) const
264 {
265  // protect against NaN
266  if (x != x || y1 != y1 || y2 != y2) return 0.;
267  // find the bin in question
268  const Int_t binx = Int_t(Double_t(nBinsX) * (x - xmin) / (xmax - xmin));
269  // get low edge of bin
270  const Double_t xlo = Double_t(nBinsX - binx) / Double_t(nBinsX) * xmin +
271  Double_t(binx) / Double_t(nBinsX) * xmax;
272  // normalise to coordinates in unit sqare
273  const Double_t hx = (x - xlo) / binSizeX;
274  // monomials
275  const Double_t hxton[4] = { hx * hx * hx, hx * hx, hx, 1. };
276  // integral
277  Double_t sum = 0.;
278  for (Int_t biny = 0; biny < nBinsY; ++biny) {
279  // get low/high edge of bin
280  const Double_t yhi = Double_t(nBinsY - biny - 1) / Double_t(nBinsY) * ymin +
281  Double_t(biny + 1) / Double_t(nBinsY) * ymax;
282  if (yhi < y1) continue;
283  const Double_t ylo = Double_t(nBinsY - biny) / Double_t(nBinsY) * ymin +
284  Double_t(biny) / Double_t(nBinsY) * ymax;
285  if (ylo > y2) break;
286  // work out integration range
287  const Double_t a = ((ylo > y1) ? 0. : (y1 - ylo)) / binSizeY;
288  const Double_t b = ((yhi < y2) ? binSizeY : (y2 - ylo)) / binSizeY;
289  // integrated monomials
290  const Double_t hyton[4] = { 0.25 * (b * b * b * b - a * a * a * a),
291  (b * b * b - a * a * a) / 3., 0.5 * (b * b - a * a), b - a };
292  Double_t lsum = 0.;
293  for (Int_t k = 0; k < NCoeff; ++k)
294  lsum += coeff(binx, biny, k) * hxton[k % 4] * hyton[k / 4];
295  sum += lsum;
296  }
297  // move from unit square coordinates to user coordinates
298  return sum * binSizeY;
299 }
300 
302  Double_t x1, Double_t x2, Double_t y1, Double_t y2) const
303 {
304  // protect against NaN
305  if (x1 != x1 || x2 != x2 || y1 != y1 || y2 != y2) return 0.;
306  // integral
307  Double_t sum = 0.;
308  for (Int_t biny = 0; biny < nBinsY; ++biny) {
309  // get low/high edge of bin
310  const Double_t yhi = Double_t(nBinsY - biny - 1) / Double_t(nBinsY) * ymin +
311  Double_t(biny + 1) / Double_t(nBinsY) * ymax;
312  if (yhi < y1) continue;
313  const Double_t ylo = Double_t(nBinsY - biny) / Double_t(nBinsY) * ymin +
314  Double_t(biny) / Double_t(nBinsY) * ymax;
315  if (ylo > y2) break;
316  // work out integration range
317  const Double_t ay = ((ylo > y1) ? 0. : (y1 - ylo)) / binSizeY;
318  const Double_t by = ((yhi < y2) ? binSizeY : (y2 - ylo)) / binSizeY;
319  const Bool_t yFullyContained = std::abs(by - ay - 1.0) < 1e-15;
320  // integrated monomials
321  const Double_t hyton[4] = {
322  0.25 * (by * by * by * by - ay * ay * ay * ay),
323  (by * by * by - ay * ay * ay) / 3., 0.5 * (by * by - ay * ay),
324  by - ay };
325  for (Int_t binx = 0; binx < nBinsX; ++binx) {
326  // get low/high edge of bin
327  const Double_t xhi = Double_t(nBinsX - binx - 1) / Double_t(nBinsX) * xmin +
328  Double_t(binx + 1) / Double_t(nBinsX) * xmax;
329  if (xhi < x1) continue;
330  const Double_t xlo = Double_t(nBinsX - binx) / Double_t(nBinsX) * xmin +
331  Double_t(binx) / Double_t(nBinsX) * xmax;
332  if (xlo > x2) break;
333  // work out integration range
334  const Double_t ax = ((xlo > x1) ? 0. : (x1 - xlo)) / binSizeX;
335  const Double_t bx = ((xhi < x2) ? binSizeX : (x2 - xlo)) / binSizeX;
336  const Bool_t xFullyContained = std::abs(bx - ax - 1.0) < 1e-15;
337  if (xFullyContained && yFullyContained) {
338  // for fully contained bins, we have cached the integral
339  sum += coeff(binx, biny, NCoeff);
340  continue;
341  }
342  // integrated monomials
343  const Double_t hxton[4] = {
344  0.25 * (bx * bx * bx * bx - ax * ax * ax * ax),
345  (bx * bx * bx - ax * ax * ax) / 3., 0.5 * (bx * bx - ax * ax),
346  bx - ax };
347  // integrate over bin
348 
349  Double_t lsum = 0.;
350  for (Int_t k = 0; k < NCoeff; ++k)
351  lsum += coeff(binx, biny, k) * hxton[k % 4] * hyton[k / 4];
352  sum += lsum;
353  }
354  }
355  // move from unit square coordinates to user coordinates
356  return sum * binSizeX * binSizeY;
357 }
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.
Double_t histcont(const TH2 &h, Int_t xbin, Int_t ybin) const
Get contents of a given bin from a histogram.
Double_t xmax
Maximum x value.
virtual Double_t analyticalIntegral() const
Evaluate analytical integral across the whole function.
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.
Int_t nBinsX
Number of bins in x axis.
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)
Int_t nBinsY
Number of bins in y axis.
virtual Double_t evaluate(Double_t x, Double_t y) const
Evaluate the function at given point.
Double_t xmin
Minimum x value.
Lau2DCubicSpline(const TH2 &h)
Constructor from histogram.
Double_t binSizeX
Bin size in x.
virtual ~Lau2DCubicSpline()
Destructor.
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)
Double_t ymax
Maximum y value.
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.
Double_t ymin
Minimum y value.
Double_t binSizeY
Bin size in y.
File containing declaration of Lau2DCubicSpline class.
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)
const Double_t & coeff(Int_t binx, Int_t biny, Int_t theCoeff) const
Const access to coefficients.