Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
CorrFit3DCF.cxx
1/*
2 * CorrFit3DCF.cxx
3 *
4 * Created on: 08-04-2015
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9
10#include "CorrFit.h"
11
12#include <TAttLine.h>
13#include <TAttMarker.h>
14#include <TAxis.h>
15#include <TCanvas.h>
16#include <TF1.h>
17#include <TF3.h>
18#include <TH1.h>
19#include <TH3.h>
20#include <TLegend.h>
21#include <TMath.h>
22#include <TNamed.h>
23#include <TString.h>
24#include <iostream>
25
26#include "Array.h"
27#include "CorrFit3DCF.h"
28#include "CorrFit3DCFPainter.h"
29#include "CorrFitHDFunc3D.h"
30#include "CorrFitMask.h"
31#include "CorrFitMask3D.h"
32#include "Cout.h"
33#include "DividedHisto.h"
34#include "Femto3DCF.h"
35#include "Splines.h"
36#include "Std.h"
37#include "StdHist.h"
38#include "StdString.h"
39
40namespace Hal {
41 CorrFit3DCF::CorrFit3DCF(e3DMode mode, Int_t parameters) : CorrFitFunc3D(mode, parameters, 3) {
42 for (int i = 0; i < 3; i++) {
43 fRange[2 * i] = 0;
44 fRange[2 * i + 1] = 1.0;
45 }
46 if (parameters < 4) {
47 Cout::PrintInfo(Form("%s must have at least 3 parameters", this->ClassName()), EInfo::kWarning);
48 return;
49 }
50 fXbins.Resize(1);
51 fYbins.Resize(1);
52 fZbins.Resize(1);
53 SetParameterName(NormID(), "N");
54 }
55
56 CorrFit3DCF::CorrFit3DCF(Int_t parameters) : CorrFit3DCF(e3DMode::kNormal3R, parameters) {
57 SetParameterName(RoutID(), "R_{out}");
58 SetParameterName(RsideID(), "R_{side}");
59 SetParameterName(RlongID(), "R_{long}");
60 }
61
62 Double_t CorrFit3DCF::EvalDenominator(Double_t x, Double_t y, Double_t z) const {
63 Int_t binx = fDenominatorHistogram->GetXaxis()->FindBin(x);
64 Int_t biny = fDenominatorHistogram->GetYaxis()->FindBin(y);
65 Int_t binz = fDenominatorHistogram->GetZaxis()->FindBin(z);
66 return fDenominatorHistogram->GetBinContent(binx, biny, binz);
67 }
68
69 Double_t CorrFit3DCF::GetFunX(Double_t* x, Double_t* params) const {
70 // co z mapowanymcorrfitem?
71 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
72 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(fYbins.Get(0));
73 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(fZbins.Get(0));
74 if (fBinCalc == kExtrapolated) { // calculate "step function"
75 Double_t X[3] = {
76 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
77 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
78 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
79 };
80 return GetScaledValue(EvalCF(X, params), params);
81 }
82
83 // simple make simple average operation
84 Double_t sum = 0;
85 Double_t ctn = fYbins.GetSize() + fZbins.GetSize();
86 Int_t binY = (fYbins.GetSize() - 1) / 2;
87 Int_t binZ = (fZbins.GetSize() - 1) / 2;
88 for (int i = 0; i < fYbins.GetSize(); i++) {
89 Double_t tempX[3] = {x[0], fYbins.Get(i), fZbins.Get(binZ)};
90 sum += CalculateCF(tempX, params);
91 }
92 for (int j = 0; j < fZbins.GetSize(); j++) {
93 Double_t tempX[3] = {x[0], fYbins.Get(binY), fZbins.Get(j)};
94 sum += CalculateCF(tempX, params);
95 }
96 return GetScaledValue(sum, params) / ctn;
97 }
98
99 Double_t CorrFit3DCF::GetFunY(Double_t* x, Double_t* params) const {
100 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
101 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(fXbins.Get(0));
102 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(fZbins.Get(0));
103 if (fBinCalc == kExtrapolated) { // calculate "step function"
104 Double_t X[3] = {
105 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
106 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
107 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
108 };
109 return GetScaledValue(EvalCF(X, params), params);
110 }
111 Double_t sum = 0;
112 Double_t ctn = fXbins.GetSize() + fZbins.GetSize();
113 Int_t binZ = (fZbins.GetSize() - 1) / 2;
114 Int_t binX = (fXbins.GetSize() - 1) / 2;
115 for (int i = 0; i < fXbins.GetSize(); i++) {
116 Double_t tempX[3] = {fXbins.Get(i), x[0], fZbins.Get(binZ)};
117 sum += GetScaledValue(CalculateCF(tempX, params), params);
118 }
119 for (int j = 0; j < fZbins.GetSize(); j++) {
120 Double_t tempX[3] = {fXbins.Get(binX), x[0], fZbins.Get(j)};
121 sum += GetScaledValue(CalculateCF(tempX, params), params);
122 }
123 return sum / ctn;
124 }
125
126 Double_t CorrFit3DCF::GetFunZ(Double_t* x, Double_t* params) const {
127 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[0]);
128 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(fXbins.Get(0));
129 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(fYbins.Get(0));
130 if (fBinCalc == kExtrapolated) { // calculate "step function"
131 Double_t X[3] = {
132 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
133 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
134 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
135 };
136 return GetScaledValue(EvalCF(X, params), params);
137 }
138 Double_t sum = 0;
139 Double_t ctn = fXbins.GetSize() + fYbins.GetSize();
140 Int_t binX = (fXbins.GetSize() - 1) / 2;
141 Int_t binY = (fYbins.GetSize() - 1) / 2;
142 for (int i = 0; i < fXbins.GetSize(); i++) {
143 Double_t tempX[3] = {fXbins.Get(i), fYbins.Get(binY), x[0]};
144 sum += GetScaledValue(CalculateCF(tempX, params), params);
145 }
146 for (int j = 0; j < fYbins.GetSize(); j++) {
147 Double_t tempX[3] = {fXbins.Get(binX), fYbins.Get(j), x[0]};
148 sum += GetScaledValue(CalculateCF(tempX, params), params);
149 }
150 return sum / ctn;
151 }
152
153 void CorrFit3DCF::SetFuncRange(Double_t x_min, Double_t x_max, Double_t y_min, Double_t y_max, Double_t z_min, Double_t z_max) {
154 fRange[0] = x_min;
155 fRange[1] = x_max;
156 fRange[2] = y_min;
157 fRange[3] = y_max;
158 fRange[4] = z_min;
159 fRange[5] = z_max;
160 }
161
162 void CorrFit3DCF::SetParametersToTF1(TF1* f) const {
163 for (int i = 0; i < GetParametersNo(); i++) {
164 f->FixParameter(i, GetParameter(i));
165 }
166 }
167
168 void CorrFit3DCF::SetErrors(TH1* num, const TH1* den) const {
169 for (int i = 0; i <= num->GetXaxis()->GetNbins(); i++) {
170 for (int j = 0; j <= num->GetYaxis()->GetNbins(); j++) {
171 for (int k = 0; k <= num->GetZaxis()->GetNbins(); k++) {
172 Double_t ea = num->GetBinError(i, j, k);
173 Double_t eb = den->GetBinError(i, j, k);
174 num->SetBinError(i, j, k, TMath::Sqrt(ea * ea + eb * eb));
175 }
176 }
177 }
178 }
179
181 Prepare();
184 fNumeratorHistogram = ((DividedHisto1D*) fCF)->GetNum();
185 fCorrelationFunctionHistogram = ((DividedHisto1D*) fCF)->GetHist(kFALSE);
186 fCorrelationFunctionHistogram->SetDirectory(0);
187 fKinematics = static_cast<Femto3DCF*>(fCF)->GetFrame();
188 if (fHDMaps == nullptr) fHDMaps = new CorrFitHDFunc3D();
189 }
190
191 double CorrFit3DCF::GetChiTFD(const double* /*par*/) const {
192 Double_t f = 0.0;
193 Double_t A, B, C;
194 Double_t e, e1, e2, chi; /* FIXME */
195 Bool_t useHD = kFALSE;
196 if (fBinCalc == kExtrapolated) useHD = kTRUE;
197 CorrFitHDFunc3D* cf = static_cast<CorrFitHDFunc3D*>(fHDMaps);
198 for (int i = 0; i < cf->GetNBins(); i++) {
199 fBinX = cf->GetXBin(i);
200 fBinY = cf->GetYBin(i);
201 fBinZ = cf->GetZBin(i);
202
203 A = fNumeratorHistogram->GetBinContent(fBinX, fBinY, fBinZ);
204 B = fDenominatorHistogram->GetBinContent(fBinX, fBinY, fBinZ);
205 e1 = fNumeratorHistogram->GetBinError(fBinX, fBinY, fBinZ);
206 e = e1 * e1;
207 C = cf->GetBinCFVal(fBinX, fBinY, fBinZ, useHD);
208 chi = A - C * B;
209 e2 = B * GetNumericalError(fBinX, fBinY, fBinZ);
210 e += e2 * e2;
211 if (e == 0) continue;
212 f += chi * chi / e;
213 }
214#ifdef CF_FIT_TRACKING
215 for (int i = 0; i < GetParametersNo(); i++) {
216 std::cout << Form("%4.3f\t", par[i]);
217 }
218 std::cout << "->" << f / z << std::endl;
219#endif
220 return f;
221 }
222
223 double CorrFit3DCF::GetChiTF(const double* /*par*/) const {
224 Double_t f = 0.0;
225 Double_t A, B, C;
226 Double_t e, e2, chi;
227 Bool_t useHD = kFALSE;
228 if (fBinCalc == kExtrapolated) useHD = kTRUE;
229 CorrFitHDFunc3D* cf = static_cast<CorrFitHDFunc3D*>(fHDMaps);
230 for (int i = 0; i < cf->GetNBins(); i++) {
231 fBinX = cf->GetXBin(i);
232 fBinY = cf->GetYBin(i);
233 fBinZ = cf->GetZBin(i);
234 A = fNumeratorHistogram->GetBinContent(fBinX, fBinY, fBinZ);
235 B = fDenominatorHistogram->GetBinContent(fBinX, fBinY, fBinZ);
236 e = fCorrelationFunctionHistogram->GetBinError(fBinX, fBinY, fBinZ);
237 e *= e;
238 C = cf->GetBinCFVal(fBinX, fBinY, fBinZ, useHD);
239 chi = A / B - C;
240 e2 = GetNumericalError(fBinX, fBinY, fBinZ);
241 e += e2 * e2;
242 if (e == 0) continue;
243 f += chi * chi / e;
244 }
245#ifdef CF_FIT_TRACKING
246 for (int i = 0; i < GetParametersNo(); i++) {
247 std::cout << Form("%4.3f\t", par[i]);
248 }
249 std::cout << "->" << f << std::endl;
250#endif
251 return f;
252 }
253
254 double CorrFit3DCF::GetLogTFD(const double* /*par*/) const {
255 Double_t f = 0.0;
256 Double_t A, B, C;
257 Bool_t useHD = kFALSE;
258 if (fBinCalc == kExtrapolated) useHD = kTRUE;
259 CorrFitHDFunc3D* cf = static_cast<CorrFitHDFunc3D*>(fHDMaps);
260 for (int i = 0; i < cf->GetNBins(); i++) {
261 fBinX = cf->GetXBin(i);
262 fBinY = cf->GetYBin(i);
263 fBinZ = cf->GetZBin(i);
264 A = fNumeratorHistogram->GetBinContent(fBinX, fBinY, fBinZ);
265 B = fDenominatorHistogram->GetBinContent(fBinX, fBinY, fBinZ);
266 C = cf->GetBinCFVal(fBinX, fBinY, fBinZ, useHD);
267 Double_t logA = (C * (A + B)) / (A * (C + 1.0));
268 Double_t logB = (A + B) / (B * (C + 1.0));
269 Double_t step = -(A * TMath::Log(logA) + B * TMath::Log(logB));
270 // Double_t step = -A*(TMath::Log(logA*TMath::Power(logB,B/A)));
271 if (TMath::IsNaN(step)) continue;
272 f += step;
273 }
274#ifdef CF_FIT_TRACKING
275 for (int i = 0; i < GetParametersNo(); i++) {
276 std::cout << Form("%4.3f\t", par[i]);
277 }
278 std::cout << "->" << f << std::endl;
279#endif
280 return f;
281 }
282
284 fActiveBins = 0;
285 if (fMask == nullptr) {
286 fMask = new CorrFitMask3D(fDenominatorHistogram->GetNbinsX(),
287 fDenominatorHistogram->GetXaxis()->GetXmin(),
288 fDenominatorHistogram->GetXaxis()->GetXmax(),
289 fDenominatorHistogram->GetNbinsY(),
290 fDenominatorHistogram->GetYaxis()->GetXmin(),
291 fDenominatorHistogram->GetYaxis()->GetXmax(),
292 fDenominatorHistogram->GetNbinsZ(),
293 fDenominatorHistogram->GetZaxis()->GetXmin(),
294 fDenominatorHistogram->GetZaxis()->GetXmax());
295 fOwnRangeMap = kFALSE;
296 }
297 if (fOwnRangeMap) { // use own map
298 } else { // get own mask
299 GetMask()->Reset(true);
300 GetMask()->ApplyThreshold(*fNumeratorHistogram, fThreshold);
301 GetMask()->ApplyThreshold(*fDenominatorHistogram, fThreshold);
302 GetMask()->ApplyRange(fRange[0], fRange[1], fRange[2], fRange[3], fRange[4], fRange[5], kFALSE);
303 }
304 GetMask()->ApplyThreshold(*fNumeratorHistogram, 0);
305 GetMask()->ApplyThreshold(*fDenominatorHistogram, 0);
306 fMask->Init();
307 fActiveBins = fMask->GetActiveBins();
308 Bool_t useHD = kFALSE;
309 if (fBinCalc == kExtrapolated) useHD = kTRUE;
310
312 Double_t free_parameters = 0;
313 for (int i = 0; i < GetParametersNo(); i++) {
314 if (!fParameters[i].IsFixed()) free_parameters++;
315 }
316 fNDF = fActiveBins - free_parameters;
317 }
318
319 void CorrFit3DCF::SetRadiusLimits(Double_t min, Double_t max) {
320 SetRoutLimits(min, max);
321 SetRsideLimits(min, max);
322 SetRlongLimits(min, max);
323 }
324
325 CorrFit3DCF::~CorrFit3DCF() {}
326
327 Double_t CorrFit3DCF::Eval(Double_t x, Double_t y, Double_t z) {
328 Double_t q[3] = {x, y, z};
329 return CalculateCF(q, fTempParamsEval);
330 }
331
332 Double_t CorrFit3DCF::EvalCF(const Double_t* xx, const Double_t* params) const {
333 if (fBinCalc == kExtrapolated) {
334 Double_t x[3];
335 Double_t y[3];
336 Double_t z[3];
337 Double_t val = 0;
338 x[1] = fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX);
339 x[0] = x[1] - fXBinf;
340 x[2] = x[1] + fXBinf;
341
342 y[1] = fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY);
343 y[0] = y[1] - fYBinf;
344 y[2] = y[1] + fYBinf;
345
346 z[1] = fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ);
347 z[0] = z[1] - fZBinf;
348 z[2] = z[1] + fZBinf;
349 Int_t xbin0 = fBinX * 2 - 2;
350 Int_t ybin0 = fBinY * 2 - 2;
351 Int_t zbin0 = fBinZ * 2 - 2;
352
353 CorrFitHDFunc3D* cf = static_cast<CorrFitHDFunc3D*>(fHDMaps);
354 Double_t X[3];
355 for (int i = 0; i < 3; i++) {
356 X[0] = x[i];
357 for (int j = 0; j < 3; j++) {
358 X[1] = y[j];
359 for (int k = 0; k < 3; k++) {
360 X[2] = z[k];
361 Double_t n_temp = cf->GetDenominatorHD()[xbin0 + i][ybin0 + j][zbin0 + k];
362 if (n_temp < 0) continue;
363 Double_t Cf = CalculateCF(X, params);
364 val += n_temp * Cf;
365 }
366 }
367 }
368 return val * cf->GetDenominatorSum()[fBinX][fBinY][fBinZ];
369 }
370 return CalculateCF(xx, params);
371 }
372
373 Double_t CorrFit3DCF::GetFunXYpp(Double_t* x, Double_t* params) const {
374 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
375 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
376 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(fZbins.Get(0));
377 if (fBinCalc == kExtrapolated) { // calculate "step function"
378 Double_t X[3] = {
379 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
380 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
381 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
382 };
383 return GetScaledValue(EvalCF(X, params), params);
384 }
385
386 // simple make simple average operation
387 // TODO take width into account
388 Double_t tempX[3] = {x[0], x[0], fZbins.Get(0)};
389 return GetScaledValue(CalculateCF(tempX, params), params);
390 }
391
392 Double_t CorrFit3DCF::GetFunXYpm(Double_t* x, Double_t* params) const {
393 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
394 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
395 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(fZbins.Get(0));
396 fBinY = fDenominatorHistogram->GetNbinsY() - fBinY + 1;
397 if (fBinCalc == kExtrapolated) { // calculate "step function"
398 Double_t X[3] = {
399 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
400 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
401 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
402 };
403 return GetScaledValue(EvalCF(X, params), params);
404 }
405
406 // simple make simple average operation
407 // TODO take width into account
408 Double_t sY = 2.0 * fYAxisf - x[0];
409 Double_t tempX[3] = {x[0], sY, fZbins.Get(0)};
410 return GetScaledValue(CalculateCF(tempX, params), params);
411 }
412
413 Double_t CorrFit3DCF::GetFunXZpp(Double_t* x, Double_t* params) const {
414 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
415 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(fYbins.Get(0));
416 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[0]);
417 if (fBinCalc == kExtrapolated) { // calculate "step function"
418 Double_t X[3] = {
419 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
420 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
421 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
422 };
423 return GetScaledValue(EvalCF(X, params), params);
424 }
425
426 // simple make simple average operation
427 // TODO take width into account
428 Double_t tempX[3] = {x[0], fYbins.Get(0), x[0]};
429 return GetScaledValue(CalculateCF(tempX, params), params);
430 }
431
432 Double_t CorrFit3DCF::GetFunXZpm(Double_t* x, Double_t* params) const {
433 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
434 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(fYbins.Get(0));
435 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[0]);
436 fBinZ = fDenominatorHistogram->GetNbinsZ() - fBinZ + 1;
437 if (fBinCalc == kExtrapolated) { // calculate "step function"
438 Double_t X[3] = {
439 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
440 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
441 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
442 };
443 return GetScaledValue(EvalCF(X, params), params);
444 }
445
446 // simple make simple average operation
447 // TODO take width into account
448 Double_t sZ = 2.0 * fZAxisf - x[0];
449 Double_t tempX[3] = {x[0], fYbins.Get(0), sZ};
450 return GetScaledValue(CalculateCF(tempX, params), params);
451 }
452
453 Double_t CorrFit3DCF::GetFunYZpp(Double_t* x, Double_t* params) const {
454 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(fXbins.Get(0));
455 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
456 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[0]);
457 if (fBinCalc == kExtrapolated) { // calculate "step function"
458 Double_t X[3] = {
459 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
460 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
461 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
462 };
463 return GetScaledValue(EvalCF(X, params), params);
464 }
465
466 // simple make simple average operation
467 // TODO take width into account
468 Double_t tempX[3] = {fXbins.Get(0), x[0], x[0]};
469 return GetScaledValue(CalculateCF(tempX, params), params);
470 }
471
472 Double_t CorrFit3DCF::GetFunYZpm(Double_t* x, Double_t* params) const {
473 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(fXbins.Get(0));
474 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
475 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[0]);
476 fBinZ = fDenominatorHistogram->GetNbinsZ() - fBinZ + 1;
477 if (fBinCalc == kExtrapolated) { // calculate "step function"
478 Double_t X[3] = {
479 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
480 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
481 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
482 };
483 return GetScaledValue(EvalCF(X, params), params);
484 }
485
486 // simple make simple average operation
487 // TODO take width into account
488 Double_t sZ = 2.0 * fZAxisf - x[0];
489 Double_t tempX[3] = {fXbins.Get(0), x[0], sZ};
490 return GetScaledValue(CalculateCF(tempX, params), params);
491 }
492
493 Double_t CorrFit3DCF::GetFunXYZppp(Double_t* x, Double_t* params) const {
494 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
495 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
496 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[0]);
497 if (fBinCalc == kExtrapolated) { // calculate "step function"
498 Double_t X[3] = {
499 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
500 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
501 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
502 };
503 return GetScaledValue(EvalCF(X, params), params);
504 }
505
506 // simple make simple average operation
507 // TODO take width into account
508 Double_t tempX[3] = {x[0], x[0], x[0]};
509 return GetScaledValue(CalculateCF(tempX, params), params);
510 }
511
512 Double_t CorrFit3DCF::GetFunXYZpmp(Double_t* x, Double_t* params) const {
513 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
514 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
515 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[0]);
516 fBinY = fDenominatorHistogram->GetNbinsY() - fBinY + 1;
517 if (fBinCalc == kExtrapolated) { // calculate "step function"
518 Double_t X[3] = {
519 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
520 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
521 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
522 };
523 return GetScaledValue(EvalCF(X, params), params);
524 }
525
526 // simple make simple average operation
527 // TODO take width into account
528 Double_t sY = 2.0 * fYAxisf - x[0];
529 Double_t tempX[3] = {x[0], sY, x[0]};
530 return GetScaledValue(CalculateCF(tempX, params), params);
531 }
532
533 Double_t CorrFit3DCF::GetFunXYZpmm(Double_t* x, Double_t* params) const {
534 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
535 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
536 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[0]);
537 fBinY = fDenominatorHistogram->GetNbinsY() - fBinY + 1;
538 fBinZ = fDenominatorHistogram->GetNbinsZ() - fBinZ + 1;
539 if (fBinCalc == kExtrapolated) { // calculate "step function"
540 Double_t X[3] = {
541 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
542 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
543 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
544 };
545 return GetScaledValue(EvalCF(X, params), params);
546 }
547
548 // simple make simple average operation
549 // TODO take width into account
550 Double_t sY = 2.0 * fYAxisf - x[0];
551 Double_t sZ = 2.0 * fZAxisf - x[0];
552 Double_t tempX[3] = {x[0], sY, sZ};
553 return GetScaledValue(CalculateCF(tempX, params), params);
554 }
555
556 Double_t CorrFit3DCF::GetFunXYZppm(Double_t* x, Double_t* params) const {
557 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
558 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
559 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[0]);
560 fBinZ = fDenominatorHistogram->GetNbinsZ() - fBinZ + 1;
561 if (fBinCalc == kExtrapolated) { // calculate "step function"
562 Double_t X[3] = {
563 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
564 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
565 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
566 };
567 return GetScaledValue(EvalCF(X, params), params);
568 }
569
570 // simple make simple average operation
571 // TODO take width into account
572 Double_t sZ = 2.0 * fZAxisf - x[0];
573 Double_t tempX[3] = {x[0], x[0], sZ};
574 return GetScaledValue(CalculateCF(tempX, params), params);
575 }
576
577 Double_t CorrFit3DCF::GetFunXY2d(Double_t* x, Double_t* params) const {
578 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
579 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[1]);
580 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(fZbins.Get(0));
581 if (fBinCalc == kExtrapolated) { // calculate "step function"
582 Double_t X[3] = {
583 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
584 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
585 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
586 };
587 return GetScaledValue(EvalCF(X, params), params);
588 }
589
590 // simple make simple average operation
591 // TODO take width into account
592 Double_t tempX[3] = {x[1], x[0], fZbins.Get(0)};
593 return GetScaledValue(CalculateCF(tempX, params), params);
594 }
595
596 Double_t CorrFit3DCF::GetFunXZ2d(Double_t* x, Double_t* params) const {
597 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(x[0]);
598 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(fYbins.Get(0));
599 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[1]);
600 if (fBinCalc == kExtrapolated) { // calculate "step function"
601 Double_t X[3] = {
602 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
603 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
604 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
605 };
606 return GetScaledValue(EvalCF(X, params), params);
607 }
608
609 // simple make simple average operation
610 // TODO take width into account
611 Double_t tempX[3] = {x[1], fYbins.Get(0), x[0]};
612 return GetScaledValue(CalculateCF(tempX, params), params);
613 }
614
615 Double_t CorrFit3DCF::GetFunYZ2d(Double_t* x, Double_t* params) const {
616 fBinX = fDenominatorHistogram->GetXaxis()->FindBin(fXbins.Get(0));
617 fBinY = fDenominatorHistogram->GetYaxis()->FindBin(x[0]);
618 fBinZ = fDenominatorHistogram->GetZaxis()->FindBin(x[1]);
619 if (fBinCalc == kExtrapolated) { // calculate "step function"
620 Double_t X[3] = {
621 fDenominatorHistogram->GetXaxis()->GetBinCenter(fBinX),
622 fDenominatorHistogram->GetYaxis()->GetBinCenter(fBinY),
623 fDenominatorHistogram->GetZaxis()->GetBinCenter(fBinZ),
624 };
625 return GetScaledValue(EvalCF(X, params), params);
626 }
627
628 // simple make simple average operation
629 // TODO take width into account
630 Double_t tempX[3] = {fXbins.Get(0), x[1], x[0]};
631 return GetScaledValue(CalculateCF(tempX, params), params);
632 }
633
634 Double_t CorrFit3DCF::GetScaledValue(Double_t x, Double_t* params) const {
635 double y = x;
636 if (params[GetParametersNo()]) return x / params[fNormParIndex];
637 return x;
638 }
639
640 void CorrFit3DCF::Calculatef(Double_t width) {
641 TH3* h = (TH3*) ((DividedHisto1D*) fCF)->GetNum();
642
643 Int_t nbins;
644 Double_t min, max;
645 Hal::Std::GetAxisPar(*h, nbins, min, max, "x");
646 fXAxisf = 0.5 * (min + max);
647 Hal::Std::GetAxisPar(*h, nbins, min, max, "y");
648 fYAxisf = 0.5 * (min + max);
649 Hal::Std::GetAxisPar(*h, nbins, min, max, "z");
650 fZAxisf = 0.5 * (min + max);
651 const Femto3DCF* CF = static_cast<Femto3DCF*>(fCF);
652 Int_t middle_x = CF->GetNum()->GetXaxis()->FindBin(0.0);
653 Int_t middle_y = CF->GetNum()->GetYaxis()->FindBin(0.0);
654 Int_t middle_z = CF->GetNum()->GetZaxis()->FindBin(0.0);
655 fXbins.Resize(1);
656 fYbins.Resize(1);
657 fZbins.Resize(1);
658 fXbins[0] = CF->GetNum()->GetXaxis()->GetBinCenter(middle_x);
659 fYbins[0] = CF->GetNum()->GetYaxis()->GetBinCenter(middle_y);
660 fZbins[0] = CF->GetNum()->GetZaxis()->GetBinCenter(middle_z);
661 if (width == -1) { // center+edges
662 fXbins.Resize(3);
663 fYbins.Resize(3);
664 fZbins.Resize(3);
665 fXbins[0] = CF->GetNum()->GetXaxis()->GetBinLowEdge(middle_x);
666 fYbins[0] = CF->GetNum()->GetYaxis()->GetBinLowEdge(middle_y);
667 fZbins[0] = CF->GetNum()->GetZaxis()->GetBinLowEdge(middle_z);
668 fXbins[1] = CF->GetNum()->GetXaxis()->GetBinCenter(middle_x);
669 fYbins[1] = CF->GetNum()->GetYaxis()->GetBinCenter(middle_y);
670 fZbins[1] = CF->GetNum()->GetZaxis()->GetBinCenter(middle_z);
671 fXbins[2] = CF->GetNum()->GetXaxis()->GetBinUpEdge(middle_x);
672 fYbins[2] = CF->GetNum()->GetYaxis()->GetBinUpEdge(middle_y);
673 fZbins[2] = CF->GetNum()->GetZaxis()->GetBinUpEdge(middle_z);
674 width = 0;
675 } else if (width != 0) {
676 fXbins.Resize(2 * width + 1);
677 fYbins.Resize(2 * width + 1);
678 fZbins.Resize(2 * width + 1);
679 for (int i = 0; i <= width * 2; i++) {
680 fXbins[i] = CF->GetNum()->GetXaxis()->GetBinCenter(middle_x - width + i);
681 fYbins[i] = CF->GetNum()->GetYaxis()->GetBinCenter(middle_y - width + i);
682 fZbins[i] = CF->GetNum()->GetZaxis()->GetBinCenter(middle_z - width + i);
683 }
684 }
685 }
686
688 CorrFitHDFunc3D* cf = static_cast<CorrFitHDFunc3D*>(fHDMaps);
689 Double_t X[3];
690 for (int a = 0; a < cf->GetBinsHDX().GetSize(); a++) {
691 Int_t i = cf->GetBinsHDX()[a];
692 Int_t j = cf->GetBinsHDY()[a];
693 Int_t k = cf->GetBinsHDZ()[a];
694 fBinX = cf->HDBinToBin(i);
695 fBinY = cf->HDBinToBin(j);
696 fBinZ = cf->HDBinToBin(k);
697 X[0] = cf->EvalHDX(i);
698 X[1] = cf->EvalHDY(j);
699 X[2] = cf->EvalHDZ(k);
700 Double_t CF = CalculateCF(X, fTempParamsEval);
701 cf->CFMapHD()[i][j][k] = CF;
702 }
703 }
704
706 const CorrFitMask3D* mask = dynamic_cast<const CorrFitMask3D*>(&map);
707 if (mask == nullptr) return;
708 if (!mask->AreCompatible(fCF)) return;
709 if (fMask) delete fMask;
710 fMask = (CorrFitMask*) mask->Clone();
711 fOwnRangeMap = kTRUE;
712 }
713
714 TF1* CorrFit3DCF::GetDrawableFunc(TString option) {
715 const Int_t nPar = GetParametersNo() + 1;
716 TString className = this->ClassName();
717 TF1* func = nullptr;
718 if (option.EqualTo("x"))
719 func = new TF1("funcX", this, &CorrFit3DCF::GetFunX, fRange[0], fRange[1], nPar, className, "GetFunX");
720 if (option.EqualTo("y"))
721 func = new TF1("funcY", this, &CorrFit3DCF::GetFunY, fRange[2], fRange[3], nPar, className, "GetFunY");
722 if (option.EqualTo("z"))
723 func = new TF1("funcZ", this, &CorrFit3DCF::GetFunZ, fRange[4], fRange[5], nPar, className, "GetFunZ");
724
725 if (option.EqualTo("xy++"))
726 func = new TF1("funcXY++", this, &CorrFit3DCF::GetFunXYpp, fRange[0], fRange[1], nPar, className, "GetFunXYpp");
727 if (option.EqualTo("xy+-"))
728 func = new TF1("funcXY+-", this, &CorrFit3DCF::GetFunXYpm, fRange[0], fRange[1], nPar, className, "GetFunXYpm");
729 if (option.EqualTo("yz++"))
730 func = new TF1("funcYZ++", this, &CorrFit3DCF::GetFunYZpp, fRange[2], fRange[3], nPar, className, "GetFunXZpp");
731 if (option.EqualTo("yz+-"))
732 func = new TF1("funcYZ+-", this, &CorrFit3DCF::GetFunYZpm, fRange[2], fRange[3], nPar, className, "GetFunYZpm");
733 if (option.EqualTo("xz++"))
734 func = new TF1("funcXZ++", this, &CorrFit3DCF::GetFunXZpp, fRange[4], fRange[5], nPar, className, "GetFunYZpp");
735 if (option.EqualTo("xz+-"))
736 func = new TF1("funcXZ+-", this, &CorrFit3DCF::GetFunXZpm, fRange[2], fRange[3], nPar, className, "GetFunXZpm");
737
738
739 if (option.EqualTo("xyz+++"))
740 func = new TF1("funcXYZ+++", this, &CorrFit3DCF::GetFunXYZppp, fRange[0], fRange[1], nPar, className, "GetFunXYZppp");
741 if (option.EqualTo("xyz+-+"))
742 func = new TF1("funcXYZ+-+", this, &CorrFit3DCF::GetFunXYZpmp, fRange[0], fRange[1], nPar, className, "GetFunXYZpmp");
743 if (option.EqualTo("xyz+--"))
744 func = new TF1("funcXYZ+--", this, &CorrFit3DCF::GetFunXYZpmm, fRange[0], fRange[1], nPar, className, "GetFunXYZpmm");
745 if (option.EqualTo("xyz++-"))
746 func = new TF1("funcXYZ++-", this, &CorrFit3DCF::GetFunXYZppm, fRange[0], fRange[1], nPar, className, "GetFunXYZppm");
747
748 if (option.EqualTo("xy"))
749 func = new TF2(
750 "fun2dxy", this, &CorrFit3DCF::GetFunXY2d, fRange[0], fRange[1], fRange[2], fRange[3], nPar, className, "GetFunXY2d");
751 if (option.EqualTo("xz"))
752 func = new TF2(
753 "fun2dxz", this, &CorrFit3DCF::GetFunXZ2d, fRange[0], fRange[1], fRange[4], fRange[5], nPar, className, "GetFunXZ2d");
754 if (option.EqualTo("yz"))
755 func = new TF2(
756 "fun2dyz", this, &CorrFit3DCF::GetFunYZ2d, fRange[2], fRange[3], fRange[4], fRange[5], nPar, className, "GetFunYZ2d");
757
758 if (!func) { std::cout << __LINE__ << __FILE__ << " wrong option " << option << std::endl; }
759 SetParametersToTF1(func);
760 func->FixParameter(GetParametersNo(), 0);
761 func->SetLineColor(GetLineColor());
762 func->SetLineStyle(GetLineStyle());
763 func->SetLineWidth(GetLineWidth());
764 return func;
765 }
766
767 void CorrFit3DCF::MakePainter(TString options) {
769 fPainter->SetOption(options);
770 }
771
772} // namespace Hal
void Resize(Int_t new_dim)
Definition Array.cxx:21
T Get(Int_t i) const
Definition Array.h:97
Int_t GetSize() const
Definition Array.h:50
void SetRadiusLimits(Double_t min, Double_t max)
virtual Double_t GetNumericalError(Int_t, Int_t, Int_t) const
virtual Double_t CalculateCF(const Double_t *x, const Double_t *params) const =0
Double_t EvalDenominator(Double_t x, Double_t y, Double_t z) const
virtual void MakePainter(TString options)
virtual Double_t EvalCF(const Double_t *x, const Double_t *params) const
virtual void RecalculateSmoothFunction() const
virtual void EstimateActiveBins()
Double_t Eval(Double_t x, Double_t y, Double_t z)
void SetErrors(TH1 *num, const TH1 *den) const
void SetFittingMask(const CorrFitMask &map)
void SetFuncRange(Double_t x_min, Double_t x_max, Double_t y_min, Double_t y_max, Double_t z_min, Double_t z_max)
void SetRoutLimits(Double_t min, Double_t max)
Int_t RsideID() const
void SetRsideLimits(Double_t min, Double_t max)
Int_t RoutID() const
Int_t RlongID() const
void SetRlongLimits(Double_t min, Double_t max)
TH1 * fDenominatorHistogram
CorrFitPainter * fPainter
TH1 * fCorrelationFunctionHistogram
CorrFitHDFunc * fHDMaps
Femto::EKinematics fKinematics
Definition CorrFitFunc.h:69
virtual void Prepare()
Double_t fActiveBins
Definition CorrFitFunc.h:93
CorrFitMask * fMask
TH1 * fNumeratorHistogram
Array_1< Double_t > fRange
Definition CorrFitFunc.h:97
Double_t GetBinCFVal(Int_t binX, Int_t binY, Int_t binZ, Bool_t extrapolated) const
Int_t GetYBin(Int_t i) const
Int_t GetXBin(Int_t i) const
Int_t GetZBin(Int_t i) const
virtual void SetMask(const CorrFitMask &mask, TH1 *denominator, Bool_t hd)=0
Int_t HDBinToBin(Int_t hd_bin) const
Width_t GetLineWidth() const
Definition CorrFit.h:137
Int_t fNDF
Definition CorrFit.h:95
Int_t GetParametersNo() const
Definition CorrFit.h:269
Double_t fThreshold
Definition CorrFit.h:99
ECalcOption fBinCalc
Definition CorrFit.h:87
Double_t GetParameter(Int_t par) const
Definition CorrFit.cxx:141
Double_t * fTempParamsEval
Definition CorrFit.h:109
Style_t GetLineStyle() const
Definition CorrFit.h:132
Color_t GetLineColor() const
Definition CorrFit.h:127
void SetParameterName(Int_t par, TString name)
Definition CorrFit.cxx:110
std::vector< FitParam > fParameters
Definition CorrFit.h:103
@ kExtrapolated
Definition CorrFit.h:79
virtual void SetOption(TString option)
Definition Painter.cxx:124