Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
Splines.cxx
1/*
2 * HalSplines.cxx
3 *
4 * Created on: 13 lis 2017
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9
10#include "Splines.h"
11
12#include "Cout.h"
13#include "Std.h"
14#include "StdString.h"
15
16#include <TAxis.h>
17#include <TH1.h>
18#include <TH2.h>
19#include <TH3.h>
20#include <TMath.h>
21#include <TMathBase.h>
22#include <TMatrixDfwd.h>
23#include <TMatrixT.h>
24#include <iostream>
25
26namespace Hal {
27 Spline1D::Spline1D(TH1* h, Double_t begval, Double_t endval) :
28 fXaxis(NULL), fNbins(0), fA(NULL), fB(NULL), fC(NULL), fAe(NULL) {
29 if (h == NULL) return;
30 fXaxis = (TAxis*) h->GetXaxis()->Clone();
31 fNbins = fXaxis->GetNbins();
32 fA = new Double_t[fNbins + 1];
33 fB = new Double_t[fNbins + 1];
34 fC = new Double_t[fNbins + 1];
35 fAe = new Double_t[fNbins + 1];
36 fBe = new Double_t[fNbins + 1];
37 SetFirstBin(h, begval, (Bool_t) begval);
38 SetLastBin(h, endval, (Bool_t) endval);
39 for (int i = 2; i < fXaxis->GetNbins(); i++) {
40 Double_t x1 = h->GetBinCenter(i - 1);
41 Double_t x2 = h->GetBinCenter(i);
42 Double_t x3 = h->GetBinCenter(i + 1);
43 Double_t y1 = h->GetBinContent(i - 1);
44 Double_t y2 = h->GetBinContent(i);
45 Double_t y3 = h->GetBinContent(i + 1);
46 Double_t A, B, C;
47 CalcParams(x1, y1, x2, y2, x3, y3, A, B, C);
48 fA[i] = A;
49 fB[i] = B;
50 fC[i] = C;
51 fAe[i] = TMath::Abs(A * h->GetBinContent(i) / h->GetBinError(i));
52 }
53 }
54
55 Double_t Spline1D::Eval(Double_t x) const {
56 Int_t bin = fXaxis->FindBin(x);
57 return fA[bin] + fB[bin] * x + fC[bin] * x * x;
58 }
59
60 Double_t Spline1D::GetError(Double_t x) const {
61 Int_t bin = fXaxis->FindBin(x);
62 return fAe[bin];
63 }
64
65 Double_t Spline1D::EvalBin(Int_t bin) const {
66 Double_t x = fXaxis->GetBinCenter(bin);
67 return fA[bin] + fB[bin] * x + fC[bin] * x * x;
68 }
69
70 Double_t Spline1D::ErrorBin(Int_t bin) const { return fAe[bin]; }
71
72 void Spline1D::Eval(Double_t x, Double_t& f, Double_t& error) const {
73 Int_t bin = fXaxis->FindBin(x);
74 f = fA[bin] + fB[bin] * x + fC[bin] * x * x;
75 error = fAe[bin];
76 }
77
78 void Spline1D::CalcParams(Double_t x1,
79 Double_t y1,
80 Double_t x2,
81 Double_t y2,
82 Double_t x3,
83 Double_t y3,
84 Double_t& A,
85 Double_t& B,
86 Double_t& C) {
87 TMatrixD AM(3, 3);
88 TMatrixD BM(3, 1);
89 AM[0][0] = 1;
90 AM[0][1] = x1;
91 AM[0][2] = x1 * x1;
92 AM[1][0] = 1;
93 AM[1][1] = x2;
94 AM[1][2] = x2 * x2;
95 AM[2][0] = 1;
96 AM[2][1] = x3;
97 AM[2][2] = x3 * x3;
98
99 BM[0][0] = y1;
100 BM[1][0] = y2;
101 BM[2][0] = y3;
102
103 TMatrixD SOL(3, 1);
104 AM = AM.InvertFast();
105
106 SOL = AM * BM;
107
108 A = SOL[0][0];
109 B = SOL[1][0];
110 C = SOL[2][0];
111 }
112
113 void Spline1D::Refit() {
114 Double_t* tempA = new Double_t[fNbins + 1];
115 Double_t* tempB = new Double_t[fNbins + 1];
116 Double_t* tempC = new Double_t[fNbins + 1];
117 for (int i = 1; i <= fNbins; i++) {
118 Double_t bin_width = fXaxis->GetBinWidth(i);
119 Double_t x2 = fXaxis->GetBinCenter(i);
120 Double_t x1 = x2 - bin_width * 0.5;
121 Double_t x3 = x2 + bin_width * 0.5;
122 Double_t y2 = Eval(x2);
123 Double_t epsilon = bin_width * 0.001;
124 Double_t y1_p = Eval(x1 + epsilon);
125 Double_t y1_m = Eval(x1 - epsilon);
126 if (i == 1) { y1_m = y1_p = Eval(x1 - bin_width * 0.5); }
127 Double_t y3_p = Eval(x3 + epsilon);
128 Double_t y3_m = Eval(x3 - epsilon);
129 if (i == fNbins) { y3_p = y3_m = Eval(x3); }
130 Double_t y1 = (y1_p + y1_m) * 0.5;
131 Double_t y3 = (y3_p + y3_m) * 0.5;
132 Double_t A, B, C;
133 CalcParams(x1, y1, x2, y2, x3, y3, A, B, C);
134 tempA[i] = A;
135 tempB[i] = B;
136 tempC[i] = C;
137 }
138 for (int i = 0; i <= fNbins; i++) {
139 fA[i] = tempA[i];
140 fB[i] = tempB[i];
141 fC[i] = tempC[i];
142 }
143 delete[] tempA;
144 delete[] tempB;
145 delete[] tempC;
146 }
147
148 void Spline1D::SetFirstBin(TH1* h, Double_t begval, Bool_t use_begval) {
149 Double_t y2 = h->GetBinContent(1);
150 Double_t x2 = h->GetBinCenter(1);
151 Double_t x1, x3;
152 Double_t y1, y3;
153 x1 = x2 - h->GetBinWidth(1) * 0.5 - h->GetBinWidth(2) * 0.5;
154 x3 = h->GetBinCenter(2);
155 y3 = h->GetBinContent(2);
156 if (use_begval) {
157 y1 = begval;
158 } else {
159 y1 = y2 - (y3 - y2);
160 }
161 Double_t A, B, C;
162 CalcParams(x1, y1, x2, y2, x3, y3, A, B, C);
163 fA[1] = A;
164 fB[1] = B;
165 fC[1] = C;
166 fAe[1] = TMath::Abs(A * h->GetBinContent(1) / h->GetBinError(1));
167 }
168
169 void Spline1D::SetLastBin(TH1* h, Double_t endval, Bool_t use_endval) {
170 Double_t x2 = h->GetBinCenter(fNbins);
171 Double_t y2 = h->GetBinContent(fNbins);
172 Double_t y1 = h->GetBinContent(fNbins - 1);
173 Double_t x1, x3, y3;
174 if (use_endval) {
175 x1 = h->GetBinCenter(fNbins - 1);
176 x3 = h->GetXaxis()->GetBinUpEdge(fNbins);
177 y3 = endval;
178 } else {
179 x1 = h->GetBinCenter(fNbins - 1);
180 x3 = x2 + h->GetBinWidth(fNbins) * 0.5 + h->GetBinWidth(fNbins - 1) * 0.5;
181 y3 = y1 + (y2 - y1);
182 }
183 Double_t A, B, C;
184 CalcParams(x1, y1, x2, y2, x3, y3, A, B, C);
185 fA[fNbins] = A;
186 fB[fNbins] = B;
187 fC[fNbins] = C;
188 fAe[fNbins] = TMath::Abs(A * h->GetBinContent(fNbins) / h->GetBinError(fNbins));
189 }
190
191 TH1D* Spline1D::GetTHD(TString name) const {
192 TH1D* res = new TH1D(name, name, fXaxis->GetNbins(), fXaxis->GetXmin(), fXaxis->GetXmax());
193 res->SetXTitle(fXaxis->GetTitle());
194 for (int i = 0; i <= fXaxis->GetNbins() + 2; i++) {
195 Double_t x = fXaxis->GetBinCenter(i);
196 // Double_t y = fA[i]+fB[i]*x+fC[i]*x*x;
197 res->SetBinContent(i, Eval(x));
198 res->SetBinError(i, GetError(x));
199 }
200 return res;
201 }
202
203 Spline1D::~Spline1D() {
204 if (fXaxis) {
205 delete fXaxis;
206 delete[] fA;
207 delete[] fB;
208 delete[] fAe;
209 delete[] fC;
210 }
211 }
212
213 Spline2D::Spline2D(TH2* h, Option_t* opt_int) : fXaxis(NULL), fYaxis(NULL), fNbinsX(0), fNbinsY(0), fA(), fAe(NULL) {
214 if (h == nullptr) {
215 Hal::Cout::Text("making copy Spline with null histogram");
216 return;
217 }
218 fXaxis = (TAxis*) h->GetXaxis()->Clone("Xspline");
219 fYaxis = (TAxis*) h->GetYaxis()->Clone("Yspline");
220 fNbinsX = fXaxis->GetNbins();
221 fNbinsY = fYaxis->GetNbins();
222 fAe = new Array_2<Double_t>();
223 fA.MakeBigger(fXaxis->GetNbins() + 2, fYaxis->GetNbins() + 2, 9);
224 fAe->MakeBigger(fXaxis->GetNbins() + 2, fYaxis->GetNbins() + 2);
225 if (fXaxis->IsVariableBinSize() && fYaxis->IsVariableBinSize()) {
226 Cout::PrintInfo("Bot axes have non-fixed bin size 2dim interpolation will probably not "
227 "work",
228 Hal::EInfo::kWarning);
229 }
230 TH2* temp_histo = (TH2*) h->Clone("temp_interpolation");
231 Extrapolate(temp_histo, opt_int);
232 Double_t Params[9];
233 for (int i = 1; i <= fNbinsX; i++) {
234 Double_t x[3];
235 x[0] = fXaxis->GetBinCenter(i - 1);
236 x[1] = fXaxis->GetBinCenter(i);
237 x[2] = fXaxis->GetBinCenter(i + 1);
238 for (int j = 1; j <= fYaxis->GetNbins(); j++) {
239 Double_t y[3];
240 y[0] = fYaxis->GetBinCenter(j - 1);
241 y[1] = fYaxis->GetBinCenter(j);
242 y[2] = fYaxis->GetBinCenter(j + 1);
243 Double_t z[3][3];
244 for (int a = -1; a < 2; a++)
245 for (int b = -1; b < 2; b++)
246 z[a + 1][b + 1] = temp_histo->GetBinContent(i + a, j + b);
247 CalcParams(x, y, z, Params);
248 (*fAe)[i][j] = temp_histo->GetBinError(i, j);
249 for (int p = 0; p < 9; p++)
250 fA[i][j][p] = Params[p];
251 }
252 }
253 // "fix edges" (but no corners)
254 for (int i = 1; i <= fNbinsX; i++) {
255 for (int p = 1; p < 9; p++) {
256 fA[i][0][p] = 0;
257 fA[i][fNbinsY][p] = 0;
258 }
259 fA[i][0][0] = fA[i][1][0];
260 fA[i][fNbinsY][0] = fA[i][fNbinsY - 1][0];
261 }
262
263 for (int i = 1; i <= fNbinsY; i++) {
264 for (int p = 1; p < 9; p++) {
265 fA[0][i][p] = 0;
266 fA[fNbinsX][i][p] = 0;
267 }
268 fA[0][i][0] = fA[1][i][0];
269 fA[fNbinsX][i][0] = fA[fNbinsX - 1][i][0];
270 }
271 delete temp_histo;
272 }
273
274 void Spline2D::Refit() {
275 Double_t Params[9];
276 Array_3<Double_t> tempA(fXaxis->GetNbins() + 2, fYaxis->GetNbins() + 2, 9);
277 for (int i = 1; i <= fNbinsX; i++) {
278 Double_t x[3];
279 Double_t dx = fXaxis->GetBinWidth(i) * 0.5;
280 x[0] = fXaxis->GetBinCenter(i) - dx;
281 x[1] = fXaxis->GetBinCenter(i);
282 x[2] = fXaxis->GetBinCenter(i) + dx;
283
284 Double_t ex = dx * 0.01;
285 Double_t x0p = x[0] + ex;
286 Double_t x0m = x[0] - ex;
287 Double_t x2p = x[2] + ex;
288 Double_t x2m = x[2] - ex;
289 Double_t X = x[1];
290
291 for (int j = 1; j <= fYaxis->GetNbins(); j++) {
292 Double_t y[3];
293 Double_t dy = fYaxis->GetBinWidth(j) * 0.5;
294 Double_t ey = dy * 0.01;
295 y[0] = fYaxis->GetBinCenter(j) - dy;
296 y[1] = fYaxis->GetBinCenter(j);
297 y[2] = fYaxis->GetBinCenter(j) + dy;
298
299 Double_t y0p = y[0] + ey;
300 Double_t y0m = y[0] - ey;
301 Double_t y2p = y[2] + ey;
302 Double_t y2m = y[2] - ey;
303 Double_t Y = y[1];
304
305 Double_t z[3][3];
306 z[0][0] = 0.25 * (Eval(x0m, y0m) + Eval(x0p, y0m) + Eval(x0m, y0p) + Eval(x0p, y0p));
307 z[0][1] = 0.5 * (Eval(x0p, Y) + Eval(x0m, Y));
308 z[0][2] = 0.25 * (Eval(x0m, y2m) + Eval(x0p, y2m) + Eval(x0m, y2p) + Eval(x0p, y2p));
309
310 z[1][0] = 0.5 * (Eval(X, y0m) + Eval(X, y0p));
311 z[1][1] = Eval(X, Y);
312 z[1][2] = 0.5 * (Eval(X, y2m) + Eval(X, y2p));
313
314 z[2][0] = 0.25 * (Eval(x2m, y0p) + Eval(x2m, y0m) + Eval(x2p, y0p) + Eval(x2p, y0m));
315 z[2][1] = 0.5 * (Eval(x2p, Y) + Eval(x2m, Y));
316 z[2][2] = 0.25 * (Eval(x2m, y2m) + Eval(x2p, y2m) + Eval(x2m, y2p) + Eval(x2p, y2p));
317
318 CalcParams(x, y, z, Params);
319 for (int p = 0; p < 9; p++) {
320 tempA[i][j][p] = Params[p];
321 }
322 }
323 }
324 fA = tempA;
325 }
326
327 Double_t Spline2D::Eval(Double_t x, Double_t y) const {
328 Int_t ix = fXaxis->FindFixBin(x);
329 Int_t iy = fYaxis->FindFixBin(y);
330 Double_t x2 = x * x;
331 Double_t y2 = y * y;
332 Double_t res = fA.Get(ix, iy, 0) + fA.Get(ix, iy, 1) * x + fA.Get(ix, iy, 2) * y + fA.Get(ix, iy, 3) * x * y
333 + fA.Get(ix, iy, 4) * x2 + fA.Get(ix, iy, 5) * y2 + fA.Get(ix, iy, 6) * x2 * y + fA.Get(ix, iy, 7) * x * y2
334 + fA.Get(ix, iy, 8) * x2 * y2;
335#ifdef SPLINE_DEBUG
336 if (TMath::IsNaN(res))
337 std::cout << "=>" << res << "(" << x << "," << y << ")\t" << fA.Get(ix, iy, 0) << " " << fA.Get(ix, iy, 1) * x << " "
338 << fA.Get(ix, iy, 2) * y + fA.Get(ix, iy, 3) * x * y << " " << fA.Get(ix, iy, 4) * x2 << " "
339 << fA.Get(ix, iy, 5) * y2 + fA.Get(ix, iy, 6) * x2 * y << " " << fA.Get(ix, iy, 7) * x * y2 << " "
340 << fA.Get(ix, iy, 8) * x2 * y2 << std::endl;
341#endif
342 return res;
343 return fA.Get(ix, iy, 0) + fA.Get(ix, iy, 1) * x + fA.Get(ix, iy, 2) * y + fA.Get(ix, iy, 3) * x * y + fA.Get(ix, iy, 4) * x2
344 + fA.Get(ix, iy, 5) * y2 + fA.Get(ix, iy, 6) * x2 * y + fA.Get(ix, iy, 7) * x * y2 + fA.Get(ix, iy, 8) * x2 * y2;
345 }
346
347 Double_t Spline2D::EvalBin(Int_t ix, Int_t iy) const {
348 Double_t x = fXaxis->GetBinCenter(ix);
349 Double_t y = fYaxis->GetBinCenter(iy);
350 Double_t x2 = x * x;
351 Double_t y2 = y * y;
352 return fA.Get(ix, iy, 0) + fA.Get(ix, iy, 1) * x + fA.Get(ix, iy, 2) * y + fA.Get(ix, iy, 3) * x * y + fA.Get(ix, iy, 4) * x2
353 + fA.Get(ix, iy, 5) * y2 + fA.Get(ix, iy, 6) * x2 * y + fA.Get(ix, iy, 7) * x * y2 + fA.Get(ix, iy, 8) * x2 * y2;
354 }
355
356 Double_t Spline2D::ErrorBin(Int_t bin_x, Int_t bin_y) const { return fAe->Get(bin_x, bin_y); }
357
358 void Spline2D::Eval(Double_t x, Double_t y, Double_t& f, Double_t& error) const {
359 Int_t ix = fXaxis->FindFixBin(x);
360 Int_t iy = fYaxis->FindFixBin(y);
361 Double_t x2 = x * x;
362 Double_t y2 = y * y;
363 f = fA.Get(ix, iy, 0) + fA.Get(ix, iy, 1) * x + fA.Get(ix, iy, 2) * y + fA.Get(ix, iy, 3) * x * y + fA.Get(ix, iy, 4) * x2
364 + fA.Get(ix, iy, 5) * y2 + fA.Get(ix, iy, 6) * x2 * y + fA.Get(ix, iy, 7) * x * y2 + fA.Get(ix, iy, 8) * x2 * y2;
365 error = fAe->Get(ix, iy);
366 }
367
368 void Spline2D::CalcParams(Double_t x[3], Double_t y[3], Double_t z[3][3], Double_t params[9]) {
369 TMatrixD AM(9, 9);
370 TMatrixD BM(9, 1);
371 for (int i = 0; i < 9; i++) {
372 Int_t xi = i % 3;
373 Int_t yi = 0;
374 if (i >= 3) yi = 1;
375 if (i >= 6) yi = 2;
376 Double_t X = x[xi];
377 Double_t Y = y[yi];
378 AM[i][0] = 1; // A
379 AM[i][1] = X; // B
380 AM[i][2] = Y; // C
381 AM[i][3] = X * Y; // D
382 AM[i][4] = X * X; // E
383 AM[i][5] = Y * Y; // F
384 AM[i][6] = X * X * Y; // G
385 AM[i][7] = X * Y * Y; // H
386 AM[i][8] = X * X * Y * Y; // I
387 BM[i][0] = z[xi][yi];
388 }
389
390 TMatrixD SOL(9, 1);
391 AM = AM.Invert();
392 SOL = AM * BM;
393
394 for (int i = 0; i < 9; i++) {
395 params[i] = SOL[i][0];
396 }
397 }
398
399 void Spline2D::PrintParams(Double_t x, Double_t y) {
400 std::cout << "PARAMS" << x << " " << y << std::endl;
401 Int_t ix = fXaxis->FindFixBin(x);
402 Int_t iy = fYaxis->FindFixBin(y);
403 for (int i = 0; i < 9; i++) {
404 std::cout << fA[ix][iy][i] << std::endl;
405 }
406 }
407
408 void Spline2D::Extrapolate(TH2* h, Option_t* interpolation) const {
409 TString option = interpolation;
410 Int_t interpolation_opt = 0;
411 if (option.EqualTo("linear")) {
412 interpolation_opt = 1;
413 } else if (option.EqualTo("const")) {
414 interpolation_opt = 2;
415 } else {
416 return;
417 }
418 for (int i = 0; i <= fNbinsX + 1; i++) { // last row & first
419 h->SetBinContent(i, 0, Extrapolate(h, i, i, i, 1, 2, 0, interpolation_opt));
420 h->SetBinError(i, 0, h->GetBinError(i, 1));
421 h->SetBinContent(i, fNbinsY + 1, Extrapolate(h, i, i, i, fNbinsY - 1, fNbinsY, fNbinsY + 1, interpolation_opt));
422 h->SetBinError(i, fNbinsY + 1, h->GetBinError(i, fNbinsY));
423 }
424 for (int i = 0; i <= fNbinsY + 1; i++) { // last column & first
425 h->SetBinContent(0, i, Extrapolate(h, 1, 2, 0, i, i, i, interpolation_opt));
426 h->SetBinError(0, i, h->GetBinError(1, i));
427 h->SetBinContent(fNbinsX + 1, i, Extrapolate(h, fNbinsX - 1, fNbinsX, fNbinsX + 1, i, i, i, interpolation_opt));
428 h->SetBinError(fNbinsX + 1, i, h->GetBinError(fNbinsX, i));
429 }
430 h->SetBinContent(0, 0, Extrapolate(h, 1, 2, 0, 1, 2, 0, interpolation_opt)); //->00
431 h->SetBinError(0, 0, h->GetBinError(1, 1));
432 h->SetBinContent(fNbinsX + 1, 0, Extrapolate(h, fNbinsX - 1, fNbinsX, fNbinsX + 1, 1, 2, 0,
433 interpolation_opt)); //->L0
434 h->SetBinError(fNbinsX + 1, 0, h->GetBinError(fNbinsX, 1));
435 h->SetBinContent(0, fNbinsY + 1, Extrapolate(h, 1, 2, 0, fNbinsY - 1, fNbinsY, fNbinsY + 1,
436 interpolation_opt)); //->0L
437 h->SetBinError(0, fNbinsY + 1, h->GetBinError(1, fNbinsY));
438 h->SetBinContent(
439 fNbinsX + 1,
440 fNbinsY + 1,
441 Extrapolate(h, fNbinsX - 1, fNbinsX, fNbinsX + 1, fNbinsY - 1, fNbinsY, fNbinsY + 1, interpolation_opt)); //->LL
442 h->SetBinError(fNbinsX + 1, fNbinsY + 1, h->GetBinError(fNbinsX, fNbinsY));
443 }
444
445 Double_t Spline2D::GetExtrapolationParam(Double_t x, Double_t y, Int_t param) const {
446 Int_t ix = fXaxis->FindFixBin(x);
447 Int_t iy = fYaxis->FindFixBin(y);
448 return fA.Get(ix, iy, param);
449 }
450
451 Double_t Spline2D::Extrapolate(TH2* h, Int_t ix1, Int_t ix2, Int_t ix, Int_t iy1, Int_t iy2, Int_t iy, Int_t opt) const {
452 Double_t z1, z2, dS, dI, Z;
453 Double_t dX = h->GetXaxis()->GetBinCenter(ix2) - h->GetXaxis()->GetBinCenter(ix1);
454 Double_t dY = h->GetYaxis()->GetBinCenter(iy2) - h->GetYaxis()->GetBinCenter(iy1);
455 dS = TMath::Sqrt(dX * dX + dY * dY);
456 Int_t closest_x = ix1;
457 Int_t closest_y = iy1;
458 if (TMath::Abs(ix1 - ix) > TMath::Abs(ix2 - ix)) closest_x = ix2;
459 if (TMath::Abs(iy1 - iy) > TMath::Abs(iy2 - iy)) closest_y = iy2;
460 dX = h->GetXaxis()->GetBinCenter(closest_x) - h->GetXaxis()->GetBinCenter(ix);
461 dY = h->GetYaxis()->GetBinCenter(closest_y) - h->GetYaxis()->GetBinCenter(iy);
462 dI = TMath::Sqrt(dX * dX + dY * dY);
463 z1 = h->GetBinContent(ix1, iy1);
464 z2 = h->GetBinContent(ix2, iy2);
465 switch (opt) {
466 case 1: {
467 Double_t dZ = (z2 - z1) * dI / dS;
468 if (ix < ix1 || iy < iy1) {
469 Z = h->GetBinContent(ix1, iy1);
470 return Z - dZ;
471 } else {
472 Z = h->GetBinContent(ix2, iy2);
473 return Z + dZ;
474 }
475 } break;
476 case 2: return h->GetBinContent(closest_x, closest_y); break;
477 default: return 0; break;
478 }
479 }
480
481 Double_t Spline2D::GetError(Double_t x, Double_t y) const {
482 Int_t ix = fXaxis->FindFixBin(x);
483 Int_t iy = fYaxis->FindFixBin(y);
484 return fAe->Get(ix, iy);
485 }
486
487 Spline2D::~Spline2D() {}
488
489 void Spline1D::FastOverwrite(TH1D* h, Double_t begval, Double_t endval) {
490 SetFirstBin(h, begval, (Bool_t) begval);
491 SetLastBin(h, endval, (Bool_t) endval);
492 for (int i = 2; i < fXaxis->GetNbins(); i++) {
493 Double_t x1 = h->GetBinCenter(i - 1);
494 Double_t x2 = h->GetBinCenter(i);
495 Double_t x3 = h->GetBinCenter(i + 1);
496 Double_t y1 = h->GetBinContent(i - 1);
497 Double_t y2 = h->GetBinContent(i);
498 Double_t y3 = h->GetBinContent(i + 1);
499 Double_t A, B, C;
500 CalcParams(x1, y1, x2, y2, x3, y3, A, B, C);
501 fA[i] = A;
502 fB[i] = B;
503 fC[i] = C;
504 fAe[i] = TMath::Abs(A * h->GetBinContent(i) / h->GetBinError(i));
505 }
506 }
507
508 void Spline2D::SetPoint(Int_t xbin, Int_t ybin, Double_t val, Double_t error) {
509 if (xbin == 0 || ybin == 0 || xbin == fXaxis->GetNbins() || ybin == fYaxis->GetNbins()) { // this is edge reset parameters
510 fA[xbin][ybin][0] = val;
511 fAe[xbin][ybin] = error;
512 for (int i = 1; i < 9; i++) {
513 fA[xbin][ybin][i] = 0;
514 }
515 }
516 }
517
518 void Spline3D::CalcParams(Double_t x[2], Double_t y[2], Double_t z[2], Double_t v[2][2][2], Double_t params[8]) {
519 TMatrixD AM(8, 8);
520 TMatrixD BM(8, 1);
521 Int_t IX[8] = {0, 1, 0, 1, 0, 1, 0, 1};
522 Int_t IY[8] = {0, 0, 1, 1, 0, 0, 1, 1};
523 Int_t IZ[8] = {0, 0, 0, 0, 1, 1, 1, 1};
524 for (int i = 0; i < 8; i++) {
525 Int_t ix = IX[i];
526 Int_t iy = IY[i];
527 Int_t iz = IZ[i];
528 Double_t X = x[ix];
529 Double_t Y = y[iy];
530 Double_t Z = z[iz];
531 AM[i][0] = 1;
532 AM[i][1] = X;
533 AM[i][2] = Y;
534 AM[i][3] = Z;
535 AM[i][4] = X * Y;
536 AM[i][5] = X * Z;
537 AM[i][6] = Y * Z;
538 AM[i][7] = X * Y * Z;
539 BM[i][0] = v[ix][iy][iz];
540 }
541 TMatrixD SOL(8, 1);
542 AM = AM.Invert();
543 SOL = AM * BM;
544
545 for (int i = 0; i < 8; i++) {
546 params[i] = SOL[i][0];
547 }
548 }
549
550 Spline3D::Spline3D(TH3* h, Option_t* interpolation) {
551 Bool_t mid = kFALSE;
552 fXaxis = (TAxis*) h->GetXaxis()->Clone();
553 fYaxis = (TAxis*) h->GetYaxis()->Clone();
554 fZaxis = (TAxis*) h->GetZaxis()->Clone();
555 fNbinsX = fXaxis->GetNbins();
556 fNbinsY = fYaxis->GetNbins();
557 fNbinsZ = fZaxis->GetNbins();
558 fA.MakeBigger(fNbinsX + 2, fNbinsY + 2, fNbinsZ + 2, 8);
559 fAe.MakeBigger(fNbinsX + 2, fNbinsY + 2, fNbinsZ + 2);
560 TString opt = interpolation;
561 Int_t fix_method = 0;
562 if (Hal::Std::FindParam(opt, "yes") || Hal::Std::FindParam(opt, "fix1")) { fix_method = 1; }
563 if (Hal::Std::FindParam(opt, "nointerpolation")) {
564 for (int i = 0; i <= fNbinsX + 1; i++) {
565 for (int j = 0; j <= fNbinsY + 1; j++) {
566 for (int k = 0; k <= fNbinsZ + 1; k++) {
567 for (int a = 1; a < 8; a++) {
568 fA[i][j][k][a] = 0;
569 }
570 fA[i][j][k][0] = h->GetBinContent(i, j, k);
571 }
572 }
573 }
574 return;
575 }
576 if (Hal::Std::FindParam(opt, "mid")) { mid = kTRUE; }
577 Double_t x[2], y[2], z[2];
578 Double_t v[2][2][2];
579 Double_t params[8];
580 Int_t idx[2] = {-1, 1};
581 for (int i = 1; i <= fNbinsX; i++) {
582 x[0] = fXaxis->GetBinLowEdge(i);
583 x[1] = fXaxis->GetBinUpEdge(i);
584 for (int j = 1; j <= fNbinsY; j++) {
585 y[0] = fYaxis->GetBinLowEdge(j);
586 y[1] = fYaxis->GetBinUpEdge(j);
587 for (int k = 1; k <= fNbinsZ; k++) {
588 z[0] = fZaxis->GetBinLowEdge(k);
589 z[1] = fZaxis->GetBinUpEdge(k);
590 Double_t C = h->GetBinContent(i, j, k);
591 switch (fix_method) {
592 case 0: // do nothing
593 {
594 for (int a = 0; a < 2; a++) {
595 Int_t I = i + idx[a];
596 for (int b = 0; b < 2; b++) {
597 Int_t J = j + idx[b];
598 for (int c = 0; c < 2; c++) {
599 Int_t K = k + idx[c];
600 v[a][b][c] = (C + h->GetBinContent(I, J, K)) * 0.5;
601 }
602 }
603 }
604 } break;
605 case 1: // mean method
606 {
607 for (int a = 0; a < 2; a++) {
608 Int_t I = i + idx[a];
609 if (I == 0) I++;
610 if (I == fNbinsX + 1) I--;
611 for (int b = 0; b < 2; b++) {
612 Int_t J = j + idx[b];
613 if (J == 0) J++;
614 if (J == fNbinsY + 1) J--;
615 for (int c = 0; c < 2; c++) {
616 Int_t K = k + idx[c];
617 if (K == 0) K++;
618 if (K == fNbinsZ + 1) K--;
619 v[a][b][c] = 0.5 * (C + h->GetBinContent(I, J, K));
620 }
621 }
622 }
623 } break;
624 case 2: // copy method
625 {
626 // TODO implement
627 } break;
628 }
629
630 CalcParams(x, y, z, v, params);
631 for (int a = 0; a < 8; a++) {
632 fA[i][j][k][a] = params[a];
633 }
634 if (mid) {
635 Double_t eval = Eval(0.5 * (x[0] + x[1]), 0.5 * (y[0] + y[1]), 0.5 * (z[0] + z[1]));
636 Double_t val = h->GetBinContent(i, j, k);
637 fA[i][j][k][0] += val - eval;
638 }
639 fAe[i][j][k] = h->GetBinError(i, j, k);
640 }
641 }
642 }
643 }
644
645 Double_t Spline3D::Eval(Double_t x, Double_t y, Double_t z) const {
646 Int_t ix = fXaxis->FindFixBin(x);
647 Int_t iy = fYaxis->FindFixBin(y);
648 Int_t iz = fZaxis->FindFixBin(z);
649 return fA.Get(ix, iy, iz, 0) + fA.Get(ix, iy, iz, 1) * x + fA.Get(ix, iy, iz, 2) * y + fA.Get(ix, iy, iz, 3) * z
650 + fA.Get(ix, iy, iz, 4) * x * y + fA.Get(ix, iy, iz, 5) * x * z + fA.Get(ix, iy, iz, 6) * y * z
651 + fA.Get(ix, iy, iz, 7) * x * y * z;
652 }
653
654 Double_t Spline3D::GetError(Double_t x, Double_t y, Double_t z) const {
655 Int_t ix = fXaxis->FindFixBin(x);
656 Int_t iy = fYaxis->FindFixBin(y);
657 Int_t iz = fZaxis->FindFixBin(z);
658 return fAe.Get(ix, iy, iz);
659 }
660
661 void Spline3D::PrintParams(Double_t x, Double_t y, Double_t z) const {
662 Int_t ix = fXaxis->FindFixBin(x);
663 Int_t iy = fYaxis->FindFixBin(y);
664 Int_t iz = fZaxis->FindFixBin(z);
665 for (int i = 0; i < 8; i++)
666 std::cout << "Par " << i << " " << fA.Get(ix, iy, iz, i) << std::endl;
667 }
668
669 void Spline3D::Eval(Double_t x, Double_t y, Double_t z, Double_t& f, Double_t& error) const {
670 Int_t ix = fXaxis->FindFixBin(x);
671 Int_t iy = fYaxis->FindFixBin(y);
672 Int_t iz = fZaxis->FindFixBin(z);
673 f = fA.Get(ix, iy, iz, 0) + fA.Get(ix, iy, iz, 1) * x + fA.Get(ix, iy, iz, 2) * y + fA.Get(ix, iy, iz, 3) * z
674 + fA.Get(ix, iy, iz, 4) * x * y + fA.Get(ix, iy, iz, 5) * x * z + fA.Get(ix, iy, iz, 6) * y * z
675 + fA.Get(ix, iy, iz, 7) * x * y * z;
676 error = fAe.Get(ix, iy, iz);
677 }
678
679 Double_t Spline3D::EvalBin(Int_t ix, Int_t iy, Int_t iz) const {
680 Double_t x = fXaxis->GetBinCenter(ix);
681 Double_t y = fYaxis->GetBinCenter(iy);
682 Double_t z = fZaxis->GetBinCenter(iz);
683 return fA.Get(ix, iy, iz, 0) + fA.Get(ix, iy, iz, 1) * x + fA.Get(ix, iy, iz, 2) * y + fA.Get(ix, iy, iz, 3) * z
684 + fA.Get(ix, iy, iz, 4) * x * y + fA.Get(ix, iy, iz, 5) * x * z + fA.Get(ix, iy, iz, 6) * y * z
685 + fA.Get(ix, iy, iz, 7) * x * y * z;
686 }
687
688 Double_t Spline3D::ErrorBin(Int_t ix, Int_t iy, Int_t iz) const { return fAe.Get(ix, iy, iz); }
689
690 Spline3D::~Spline3D() {
691 if (fXaxis) {
692 delete fXaxis;
693 delete fYaxis;
694 delete fZaxis;
695 }
696 }
697
698 TH2D* Spline2D::GetTHD(TString name) const {
699 TH2D* res = new TH2D(name,
700 name,
701 fXaxis->GetNbins(),
702 fXaxis->GetXmin(),
703 fXaxis->GetXmax(),
704 fYaxis->GetNbins(),
705 fYaxis->GetXmin(),
706 fYaxis->GetXmax());
707 res->SetXTitle(fXaxis->GetTitle());
708 res->SetYTitle(fYaxis->GetTitle());
709
710 for (int i = 0; i <= fXaxis->GetNbins() + 1; i++) {
711 Double_t x = fXaxis->GetBinCenter(i);
712 for (int j = 0; j <= fYaxis->GetNbins() + 1; j++) {
713 Double_t y = fYaxis->GetBinCenter(j);
714 res->SetBinContent(i, j, Eval(x, y));
715 res->SetBinError(i, j, GetError(x, y));
716 }
717 }
718 return res;
719 }
720
721 TH3D* Spline3D::GetTHD(TString name) const {
722 TH3D* res = new TH3D(name,
723 name,
724 fXaxis->GetNbins(),
725 fXaxis->GetXmin(),
726 fXaxis->GetXmax(),
727 fYaxis->GetNbins(),
728 fYaxis->GetXmin(),
729 fYaxis->GetXmax(),
730 fZaxis->GetNbins(),
731 fZaxis->GetXmin(),
732 fZaxis->GetXmax());
733 res->SetXTitle(fXaxis->GetTitle());
734 res->SetYTitle(fYaxis->GetTitle());
735 res->SetZTitle(fZaxis->GetTitle());
736 for (int i = 0; i <= fXaxis->GetNbins() + 1; i++) {
737 Double_t x = fXaxis->GetBinCenter(i);
738 for (int j = 0; j <= fYaxis->GetNbins() + 1; j++) {
739 Double_t y = fYaxis->GetBinCenter(j);
740 for (int k = 0; k <= fZaxis->GetNbins() + 1; k++) {
741 Double_t z = fZaxis->GetBinCenter(k);
742 res->SetBinContent(i, j, Eval(x, y, z));
743 res->SetBinError(i, j, k, GetError(x, y, z));
744 }
745 }
746 }
747 return res;
748 }
749
750 void Spline3D::Refit() {
751 Array_4<Double_t> ANew;
752 ANew.MakeBigger(fNbinsX + 2, fNbinsY + 2, fNbinsZ + 2, 8);
753 Double_t x[2], y[2], z[2];
754 Double_t v[2][2][2];
755 Double_t params[8];
756
757 for (int i = 1; i <= fNbinsX; i++) {
758 x[0] = fXaxis->GetBinLowEdge(i);
759 x[1] = fXaxis->GetBinUpEdge(i);
760 Double_t dx = fXaxis->GetBinWidth(i) * 0.1;
761 for (int j = 1; j <= fNbinsY; j++) {
762 y[0] = fYaxis->GetBinLowEdge(j);
763 y[1] = fYaxis->GetBinUpEdge(j);
764 Double_t dy = fYaxis->GetBinWidth(j) * 0.1;
765 for (int k = 1; k <= fNbinsZ; k++) {
766 z[0] = fZaxis->GetBinLowEdge(k);
767 z[1] = fZaxis->GetBinUpEdge(k);
768 Double_t dz = fZaxis->GetBinWidth(k) * 0.1;
769 for (int a = 0; a < 2; a++) {
770 for (int b = 0; b < 2; b++) {
771 for (int c = 0; c < 2; c++) {
772 v[a][b][c] = 0;
773 Double_t sum = 0;
774
775 // calcuate "averaged corner"
776 Int_t D[2], E[2], F[2];
777 D[0] = fXaxis->FindBin(x[a] - dx);
778 D[1] = fXaxis->FindBin(x[a] + dx);
779 D[0] = TMath::Max(1, D[0]);
780 D[1] = TMath::Min(fNbinsX, D[1]);
781
782 E[0] = fYaxis->FindBin(y[b] - dy);
783 E[1] = fYaxis->FindBin(y[b] + dy);
784 E[0] = TMath::Max(1, E[0]);
785 E[1] = TMath::Min(fNbinsY, E[1]);
786
787 F[0] = fZaxis->FindBin(z[c] - dz);
788 F[1] = fZaxis->FindBin(z[c] + dz);
789 F[0] = TMath::Max(1, F[0]);
790 F[1] = TMath::Min(fNbinsZ, F[1]);
791
792 for (int d = 0; d < 2; d++) {
793 for (int e = 0; e < 2; e++) {
794 for (int f = 0; f < 2; f++) {
795 sum += EvalBin(D[d], E[e], F[f], x[a], y[b], z[c]);
796 }
797 }
798 }
799
800 v[a][b][c] = sum * 0.125;
801 }
802 }
803 }
804 CalcParams(x, y, z, v, params);
805 for (int a = 0; a < 8; a++) {
806 ANew[i][j][k][a] = params[a];
807 }
808 }
809 }
810 }
811 fA = ANew;
812 }
813
814 Double_t Spline3D::EvalBin(Int_t ix, Int_t iy, Int_t iz, Double_t x, Double_t y, Double_t z) const {
815 return fA.Get(ix, iy, iz, 0) + fA.Get(ix, iy, iz, 1) * x + fA.Get(ix, iy, iz, 2) * y + fA.Get(ix, iy, iz, 3) * z
816 + fA.Get(ix, iy, iz, 4) * x * y + fA.Get(ix, iy, iz, 5) * x * z + fA.Get(ix, iy, iz, 6) * y * z
817 + fA.Get(ix, iy, iz, 7) * x * y * z;
818 }
819
820 void Spline3D::FastOverwrite(TH3* h, Option_t* interpolation) {
821 Bool_t mid = kFALSE;
822 TString opt = interpolation;
823 Int_t fix_method = 0;
824 if (Hal::Std::FindParam(opt, "yes") || Hal::Std::FindParam(opt, "fix1")) { fix_method = 1; }
825 if (Hal::Std::FindParam(opt, "nointerpolation")) {
826 for (int i = 0; i <= fNbinsX + 1; i++) {
827 for (int j = 0; j <= fNbinsY + 1; j++) {
828 for (int k = 0; k <= fNbinsZ + 1; k++) {
829 for (int a = 1; a < 8; a++) {
830 fA[i][j][k][a] = 0;
831 }
832 fA[i][j][k][0] = h->GetBinContent(i, j, k);
833 }
834 }
835 }
836 return;
837 }
838 if (Hal::Std::FindParam(opt, "mid")) { mid = kTRUE; }
839 Double_t x[2], y[2], z[2];
840 Double_t v[2][2][2];
841 Double_t params[8];
842 Int_t idx[2] = {-1, 1};
843 for (int i = 1; i <= fNbinsX; i++) {
844 x[0] = fXaxis->GetBinLowEdge(i);
845 x[1] = fXaxis->GetBinUpEdge(i);
846 for (int j = 1; j <= fNbinsY; j++) {
847 y[0] = fYaxis->GetBinLowEdge(j);
848 y[1] = fYaxis->GetBinUpEdge(j);
849 for (int k = 1; k <= fNbinsZ; k++) {
850 z[0] = fZaxis->GetBinLowEdge(k);
851 z[1] = fZaxis->GetBinUpEdge(k);
852 Double_t C = h->GetBinContent(i, j, k);
853 switch (fix_method) {
854 case 0: // do nothing
855 {
856 for (int a = 0; a < 2; a++) {
857 Int_t I = i + idx[a];
858 for (int b = 0; b < 2; b++) {
859 Int_t J = j + idx[b];
860 for (int c = 0; c < 2; c++) {
861 Int_t K = k + idx[c];
862 v[a][b][c] = (C + h->GetBinContent(I, J, K)) * 0.5;
863 }
864 }
865 }
866 } break;
867 case 1: // mean method
868 {
869 for (int a = 0; a < 2; a++) {
870 Int_t I = i + idx[a];
871 if (I == 0) I++;
872 if (I == fNbinsX + 1) I--;
873 for (int b = 0; b < 2; b++) {
874 Int_t J = j + idx[b];
875 if (J == 0) J++;
876 if (J == fNbinsY + 1) J--;
877 for (int c = 0; c < 2; c++) {
878 Int_t K = k + idx[c];
879 if (K == 0) K++;
880 if (K == fNbinsZ + 1) K--;
881 v[a][b][c] = 0.5 * (C + h->GetBinContent(I, J, K));
882 }
883 }
884 }
885 } break;
886 case 2: // copy method
887 {
888 // TODO implement
889 } break;
890 }
891
892 CalcParams(x, y, z, v, params);
893 for (int a = 0; a < 8; a++) {
894 fA[i][j][k][a] = params[a];
895 }
896 if (mid) {
897 Double_t eval = Eval(0.5 * (x[0] + x[1]), 0.5 * (y[0] + y[1]), 0.5 * (z[0] + z[1]));
898 Double_t val = h->GetBinContent(i, j, k);
899 fA[i][j][k][0] += val - eval;
900 }
901 fAe[i][j][k] = h->GetBinError(i, j, k);
902 }
903 }
904 }
905 }
906} // namespace Hal
Hal::Array_2
Definition Array.h:127
Hal::Array_2::MakeBigger
void MakeBigger(Int_t sizeA, Int_t sizeB)
Definition Array.cxx:197
Hal::Array_2::Get
T Get(Int_t i, Int_t j) const
Definition Array.h:174
Hal::Array_3
Definition Array.h:199
Hal::Array_3::MakeBigger
void MakeBigger(Int_t sizeA, Int_t sizeB, Int_t sizeC)
Definition Array.cxx:263
Hal::Array_3::Get
T Get(Int_t A, Int_t B, Int_t C) const
Definition Array.h:248
Hal::Array_4
Definition Array.h:276
Hal::Array_4::Get
T Get(Int_t A, Int_t B, Int_t C, Int_t D) const
Definition Array.h:327
Hal::Array_4::MakeBigger
void MakeBigger(Int_t sizeA, Int_t sizeB, Int_t sizeC, Int_t sizeD)
Definition Array.cxx:326
Hal::Cout::Text
static void Text(TString text, TString option="L", Color_t color=-1)
Definition Cout.cxx:92
Hal::Spline1D::FastOverwrite
void FastOverwrite(TH1D *h, Double_t begval=0, Double_t endval=0)
Definition Splines.cxx:489
Hal::Spline2D::PrintParams
void PrintParams(Double_t x, Double_t y)
Definition Splines.cxx:399
Hal::Spline2D::GetTHD
virtual TH2D * GetTHD(TString name="spline_thd") const
Definition Splines.cxx:698
Hal::Spline2D::Eval
Double_t Eval(Double_t x, Double_t y) const
Definition Splines.cxx:327
Hal::Spline2D::ErrorBin
Double_t ErrorBin(Int_t bin_x, Int_t bin_y) const
Definition Splines.cxx:356
Hal::Spline2D::EvalBin
Double_t EvalBin(Int_t bin_x, Int_t bin_y) const
Definition Splines.cxx:347
Hal::Spline2D::SetPoint
void SetPoint(Int_t xbin, Int_t ybin, Double_t val, Double_t err=0)
Definition Splines.cxx:508
Hal::Spline2D::GetError
Double_t GetError(Double_t x, Double_t y) const
Definition Splines.cxx:481
Hal::Spline2D::GetExtrapolationParam
Double_t GetExtrapolationParam(Double_t x, Double_t y, Int_t param) const
Definition Splines.cxx:445
Hal::Spline3D::PrintParams
void PrintParams(Double_t x, Double_t y, Double_t z) const
Definition Splines.cxx:661
Hal::Spline3D::Spline3D
Spline3D(TH3 *h=NULL, Option_t *interpolation="")
Definition Splines.cxx:550
Hal::Spline3D::GetTHD
virtual TH3D * GetTHD(TString name="spline_thd") const
Definition Splines.cxx:721
Hal::Spline3D::ErrorBin
Double_t ErrorBin(Int_t x, Int_t y, Int_t z) const
Definition Splines.cxx:688
Hal::Spline3D::Eval
Double_t Eval(Double_t x, Double_t y, Double_t z) const
Definition Splines.cxx:645
Hal::Spline3D::GetError
Double_t GetError(Double_t x, Double_t y, Double_t z) const
Definition Splines.cxx:654
Hal::Spline3D::FastOverwrite
void FastOverwrite(TH3 *h, Option_t *interpolation="")
Definition Splines.cxx:820