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