]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGLF/SPECTRA/PiKaPr/COMBINED/SpectraUtils.C
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
spectra-utils update
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / COMBINED / SpectraUtils.C
... / ...
CommitLineData
1/*
2 several function used for PbPb combined spectra
3 Blast Wave is also implemented here
4 further documentation will come
5
6 author: Roberto Preghenella
7 email : preghenella@bo.infn.it
8*/
9
10
11/*****************************************************************/
12/* BOLTZMANN
13/*****************************************************************/
14
15Double_t
16Boltzmann_Func(const Double_t *x, const Double_t *p)
17{
18 /* dN/dpt */
19
20 Double_t pt = x[0];
21 Double_t mass = p[0];
22 Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
23 Double_t T = p[1];
24 Double_t norm = p[2];
25
26 return pt * norm * mt * TMath::Exp(-mt / T);
27}
28
29TF1 *
30Boltzmann(const Char_t *name, Double_t mass, Double_t T = 0.1, Double_t norm = 1.)
31{
32
33 TF1 *fBoltzmann = new TF1(name, Boltzmann_Func, 0., 10., 3);
34 fBoltzmann->SetParameters(mass, T, norm);
35 fBoltzmann->SetParNames("mass", "T", "norm");
36 fBoltzmann->FixParameter(0, mass);
37 return fBoltzmann;
38}
39
40/*****************************************************************/
41/* LEVY-TSALLIS */
42/*****************************************************************/
43
44Double_t
45LevyTsallis_Func(const Double_t *x, const Double_t *p)
46{
47 /* dN/dpt */
48
49 Double_t pt = x[0];
50 Double_t mass = p[0];
51 Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
52 Double_t n = p[1];
53 Double_t C = p[2];
54 Double_t norm = p[3];
55
56 Double_t part1 = (n - 1.) * (n - 2);
57 Double_t part2 = n * C * (n * C + mass * (n - 2.));
58 Double_t part3 = part1 / part2;
59 Double_t part4 = 1. + (mt - mass) / n / C;
60 Double_t part5 = TMath::Power(part4, -n);
61 return pt * norm * part3 * part5;
62}
63
64TF1 *
65LevyTsallis(const Char_t *name, Double_t mass, Double_t n = 5., Double_t C = 0.1, Double_t norm = 1.)
66{
67
68 TF1 *fLevyTsallis = new TF1(name, LevyTsallis_Func, 0., 10., 4);
69 fLevyTsallis->SetParameters(mass, n, C, norm);
70 fLevyTsallis->SetParNames("mass", "n", "C", "norm");
71 fLevyTsallis->FixParameter(0, mass);
72 return fLevyTsallis;
73}
74
75/*****************************************************************/
76/* BOLTZMANN-GIBBS BLAST-WAVE */
77/*****************************************************************/
78
79static TF1 *fBGBlastWave_Integrand = NULL;
80Double_t
81BGBlastWave_Integrand(const Double_t *x, const Double_t *p)
82{
83
84 /*
85 x[0] -> r (radius)
86 p[0] -> mT (transverse mass)
87 p[1] -> pT (transverse momentum)
88 p[2] -> beta_max (surface velocity)
89 p[3] -> T (freezout temperature)
90 p[4] -> n (velocity profile)
91 */
92
93 Double_t r = x[0];
94 Double_t mt = p[0];
95 Double_t pt = p[1];
96 Double_t beta_max = p[2];
97 Double_t temp_1 = 1. / p[3];
98 Double_t n = p[4];
99
100 Double_t beta = beta_max * TMath::Power(r, n);
101 if (beta > 0.9999999999999999) beta = 0.9999999999999999;
102 Double_t rho = TMath::ATanH(beta);
103 Double_t argI0 = pt * TMath::SinH(rho) * temp_1;
104 if (argI0 > 700.) argI0 = 700.;
105 Double_t argK1 = mt * TMath::CosH(rho) * temp_1;
106 // if (argI0 > 100 || argI0 < -100)
107 // printf("r=%f, pt=%f, beta_max=%f, temp=%f, n=%f, mt=%f, beta=%f, rho=%f, argI0=%f, argK1=%f\n", r, pt, beta_max, 1. / temp_1, n, mt, beta, rho, argI0, argK1);
108 return r * mt * TMath::BesselI0(argI0) * TMath::BesselK1(argK1);
109
110}
111
112Double_t
113BGBlastWave_Func(const Double_t *x, const Double_t *p)
114{
115 /* dN/dpt */
116
117 Double_t pt = x[0];
118 Double_t mass = p[0];
119 Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
120 Double_t beta_max = p[1];
121 Double_t temp = p[2];
122 Double_t n = p[3];
123 Double_t norm = p[4];
124
125 if (!fBGBlastWave_Integrand)
126 fBGBlastWave_Integrand = new TF1("fBGBlastWave_Integrand", BGBlastWave_Integrand, 0., 1., 5);
127 fBGBlastWave_Integrand->SetParameters(mt, pt, beta_max, temp, n);
128 Double_t integral = fBGBlastWave_Integrand->Integral(0., 1.);
129 return norm * pt * integral;
130}
131
132TF1 *
133BGBlastWave(const Char_t *name, Double_t mass, Double_t beta_max = 0.9, Double_t temp = 0.1, Double_t n = 1., Double_t norm = 1.e6)
134{
135
136 TF1 *fBGBlastWave = new TF1(name, BGBlastWave_Func, 0., 10., 5);
137 fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm);
138 fBGBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
139 fBGBlastWave->FixParameter(0, mass);
140 fBGBlastWave->SetParLimits(1, 0.01, 0.99);
141 fBGBlastWave->SetParLimits(2, 0.01, 1.);
142 fBGBlastWave->SetParLimits(3, 0.1, 10.);
143 return fBGBlastWave;
144}
145
146/*****************************************************************/
147/* TSALLIS BLAST-WAVE */
148/*****************************************************************/
149
150static TF1 *fTsallisBlastWave_Integrand_r = NULL;
151Double_t
152TsallisBlastWave_Integrand_r(const Double_t *x, const Double_t *p)
153{
154 /*
155 x[0] -> r (radius)
156 p[0] -> mT (transverse mass)
157 p[1] -> pT (transverse momentum)
158 p[2] -> beta_max (surface velocity)
159 p[3] -> T (freezout temperature)
160 p[4] -> n (velocity profile)
161 p[5] -> q
162 p[6] -> y (rapidity)
163 p[7] -> phi (azimuthal angle)
164 */
165
166 Double_t r = x[0];
167 Double_t mt = p[0];
168 Double_t pt = p[1];
169 Double_t beta_max = p[2];
170 Double_t temp_1 = 1. / p[3];
171 Double_t n = p[4];
172 Double_t q = p[5];
173 Double_t y = p[6];
174 Double_t phi = p[7];
175
176 if (q <= 1.) return r;
177
178 Double_t beta = beta_max * TMath::Power(r, n);
179 Double_t rho = TMath::ATanH(beta);
180
181 Double_t part1 = mt * TMath::CosH(y) * TMath::CosH(rho);
182 Double_t part2 = pt * TMath::SinH(rho) * TMath::Cos(phi);
183 Double_t part3 = part1 - part2;
184 Double_t part4 = 1 + (q - 1.) * temp_1 * part3;
185 Double_t expo = -1. / (q - 1.);
186 // printf("part1=%f, part2=%f, part3=%f, part4=%f, expo=%f\n", part1, part2, part3, part4, expo);
187 Double_t part5 = TMath::Power(part4, expo);
188
189 return r * part5;
190}
191
192static TF1 *fTsallisBlastWave_Integrand_phi = NULL;
193Double_t
194TsallisBlastWave_Integrand_phi(const Double_t *x, const Double_t *p)
195{
196 /*
197 x[0] -> phi (azimuthal angle)
198 */
199
200 Double_t phi = x[0];
201 fTsallisBlastWave_Integrand_r->SetParameter(7, phi);
202 Double_t integral = fTsallisBlastWave_Integrand_r->Integral(0., 1.);
203 return integral;
204}
205
206static TF1 *fTsallisBlastWave_Integrand_y = NULL;
207Double_t
208TsallisBlastWave_Integrand_y(const Double_t *x, const Double_t *p)
209{
210 /*
211 x[0] -> y (rapidity)
212 */
213
214 Double_t y = x[0];
215 fTsallisBlastWave_Integrand_r->SetParameter(6, y);
216 Double_t integral = fTsallisBlastWave_Integrand_phi->Integral(-TMath::Pi(), TMath::Pi());
217 return TMath::CosH(y) * integral;
218}
219
220Double_t
221TsallisBlastWave_Func(const Double_t *x, const Double_t *p)
222{
223 /* dN/dpt */
224
225 Double_t pt = x[0];
226 Double_t mass = p[0];
227 Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
228 Double_t beta_max = p[1];
229 Double_t temp = p[2];
230 Double_t n = p[3];
231 Double_t q = p[4];
232 Double_t norm = p[5];
233
234 if (!fTsallisBlastWave_Integrand_r)
235 fTsallisBlastWave_Integrand_r = new TF1("fTsallisBlastWave_Integrand_r", TsallisBlastWave_Integrand_r, 0., 1., 8);
236 if (!fTsallisBlastWave_Integrand_phi)
237 fTsallisBlastWave_Integrand_phi = new TF1("fTsallisBlastWave_Integrand_phi", TsallisBlastWave_Integrand_phi, -TMath::Pi(), TMath::Pi(), 0);
238 if (!fTsallisBlastWave_Integrand_y)
239 fTsallisBlastWave_Integrand_y = new TF1("fTsallisBlastWave_Integrand_y", TsallisBlastWave_Integrand_y, -0.5, 0.5, 0);
240
241 fTsallisBlastWave_Integrand_r->SetParameters(mt, pt, beta_max, temp, n, q, 0., 0.);
242 Double_t integral = fTsallisBlastWave_Integrand_y->Integral(-0.5, 0.5);
243 return norm * pt * integral;
244}
245
246TF1 *
247TsallisBlastWave(const Char_t *name, Double_t mass, Double_t beta_max = 0.9, Double_t temp = 0.1, Double_t n = 1., Double_t q = 2., Double_t norm = 1.e6)
248{
249
250 TF1 *fTsallisBlastWave = new TF1(name, TsallisBlastWave_Func, 0., 10., 6);
251 fTsallisBlastWave->SetParameters(mass, beta_max, temp, n, q, norm);
252 fTsallisBlastWave->SetParNames("mass", "beta_max", "T", "n", "q", "norm");
253 fTsallisBlastWave->FixParameter(0, mass);
254 fTsallisBlastWave->SetParLimits(1, 0.01, 0.99);
255 fTsallisBlastWave->SetParLimits(2, 0.01, 1.);
256 fTsallisBlastWave->SetParLimits(3, 0.1, 10.);
257 fTsallisBlastWave->SetParLimits(4, 1., 10.);
258 return fTsallisBlastWave;
259}
260
261/*****************************************************************/
262/*****************************************************************/
263/*****************************************************************/
264
265
266TF1 *
267BGBlastWave_SingleFit(TH1 *h, Double_t mass, Option_t *opt = "")
268{
269
270 TF1 *f = BGBlastWave(Form("fBGBW_%s", h->GetName()), mass);
271 h->Fit(f);
272 h->Fit(f);
273 h->Fit(f, opt);
274 return f;
275
276}
277
278Int_t nBW;
279TF1 *fBGBW[1000];
280TGraphErrors *gBW[1000];
281
282TObjArray *
283BGBlastWave_GlobalFit(TObjArray *data, Double_t *mass, Double_t profile = 1., Bool_t fixProfile = kFALSE)
284{
285
286 /* get data */
287 nBW = data->GetEntries();
288 for (Int_t idata = 0; idata < nBW; idata++) {
289 gBW[idata] = (TGraphErrors *)data->At(idata);
290 gBW[idata]->SetName(Form("gBW%d", idata));
291 }
292
293 /* init BG blast-wave functions */
294 for (Int_t idata = 0; idata < nBW; idata++) {
295 printf("init BG-BlastWave function #%d: mass = %f\n", idata, mass[idata]);
296 fBGBW[idata] = BGBlastWave(Form("fBGBW%d", idata), mass[idata]);
297 }
298
299 /* display data */
300 TCanvas *cBW = new TCanvas("cBW");
301 cBW->Divide(nBW, 1);
302 for (Int_t idata = 0; idata < nBW; idata++) {
303 cBW->cd(idata + 1);
304 gBW[idata]->Draw("ap*");
305 }
306 cBW->Update();
307
308 /* init minuit: nBW normalizations + 3 (beta, T, n) BG-BlastWave params */
309 const Int_t nbwpars = 3;
310 const Int_t nfitpars = nBW + nbwpars;
311 TMinuit *minuit = new TMinuit(nfitpars);
312 minuit->SetFCN(BGBlastWave_FCN);
313 Double_t arglist[10];
314 Int_t ierflg = 0;
315 arglist[0] = 1;
316 minuit->mnexcm("SET ERR", arglist, 1, ierflg);
317 for (Int_t idata = 0; idata < nBW; idata++)
318 minuit->mnparm(idata, Form("norm%d", idata), 1.e6, 1., 0., 0., ierflg);
319 minuit->mnparm(nBW + 0, "<beta>", 0.5, 0.1, 0., 1., ierflg);
320 minuit->mnparm(nBW + 1, "T", 0.2, 0.1, 0., 1., ierflg);
321 minuit->mnparm(nBW + 2, "n", profile, 0.1, 0., 10., ierflg);
322 if (fixProfile) minuit->FixParameter(nBW + 2);
323
324 /* set strategy */
325 arglist[0] = 1;
326 minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
327
328 /* start MIGRAD minimization */
329 arglist[0] = 500000;
330 arglist[1] = 1.;
331 minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
332
333 /* set strategy */
334 arglist[0] = 2;
335 minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
336
337 /* start MIGRAD minimization */
338 arglist[0] = 500000;
339 arglist[1] = 1.;
340 minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
341
342 /* start IMPROVE minimization */
343 arglist[0] = 500000;
344 minuit->mnexcm("IMPROVE", arglist, 1, ierflg);
345
346 /* start MINOS */
347 arglist[0] = 500000;
348 arglist[1] = nBW + 1;
349 arglist[2] = nBW + 2;
350 arglist[3] = nBW + 3;
351 minuit->mnexcm("MINOS", arglist, 4, ierflg);
352
353 /* print results */
354 Double_t amin,edm,errdef;
355 Int_t nvpar,nparx,icstat;
356 minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
357 minuit->mnprin(4, amin);
358
359 /* get parameters */
360 Double_t beta, betae, betaeplus, betaeminus, betagcc, temp, tempe, tempeplus, tempeminus, tempgcc, prof, profe, profeplus, profeminus, profgcc;
361 minuit->GetParameter(nBW + 0, beta, betae);
362 minuit->mnerrs(nBW + 0, betaeplus, betaeminus, betae, betagcc);
363 minuit->GetParameter(nBW + 1, temp, tempe);
364 minuit->mnerrs(nBW + 1, tempeplus, tempeminus, tempe, tempgcc);
365 minuit->GetParameter(nBW + 2, prof, profe);
366 minuit->mnerrs(nBW + 2, profeplus, profeminus, profe, profgcc);
367 Double_t beta_max = 0.5 * (2. + prof) * beta;
368 Double_t norm[1000], norme[1000];
369 for (Int_t idata = 0; idata < nBW; idata++)
370 minuit->GetParameter(idata, norm[idata], norme[idata]);
371
372 /* printout */
373 printf("*********************************\n");
374 printf("beta_max = %f\n", beta_max);
375 printf("<beta> = %f +- %f (e+ = %f, e- = %f)\n", beta, betae, betaeplus, betaeminus);
376 printf("T = %f +- %f (e+ = %f, e- = %f)\n", temp, tempe, tempeplus, tempeminus);
377 printf("n = %f +- %f (e+ = %f, e- = %f)\n", prof, profe, profeplus, profeminus);
378
379 /* 1-sigma contour */
380 minuit->SetErrorDef(1);
381 TGraph *gCont1 = NULL;
382 gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
383 if (gCont1) gCont1->SetName("gCont1");
384
385 /* 2-sigma contour */
386 minuit->SetErrorDef(4);
387 TGraph *gCont2 = NULL;
388 // gCont2 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
389 if (gCont2) gCont2->SetName("gCont2");
390
391 /* display fits */
392 for (Int_t idata = 0; idata < nBW; idata++) {
393 cBW->cd(idata + 1);
394 fBGBW[idata]->SetParameter(4, norm[idata]);
395 fBGBW[idata]->SetParameter(1, beta_max);
396 fBGBW[idata]->SetParameter(2, temp);
397 fBGBW[idata]->SetParameter(3, prof);
398 fBGBW[idata]->Draw("same");
399 }
400 cBW->Update();
401
402 /* histo params */
403 TH1D *hBW = new TH1D("hBW", "", 3, 0., 3.);
404 hBW->SetBinContent(1, beta);
405 hBW->SetBinError(1, betae);
406 hBW->SetBinContent(2, temp);
407 hBW->SetBinError(2, tempe);
408 hBW->SetBinContent(3, prof);
409 hBW->SetBinError(3, profe);
410
411 /* BW graph */
412 TGraphAsymmErrors *gBetaT = new TGraphAsymmErrors();
413 gBetaT->SetName("gBetaT");
414 gBetaT->SetPoint(0, beta, temp);
415 gBetaT->SetPointEXlow(0, TMath::Abs(betaeminus));
416 gBetaT->SetPointEXhigh(0, TMath::Abs(betaeplus));
417 gBetaT->SetPointEYlow(0, TMath::Abs(tempeminus));
418 gBetaT->SetPointEYhigh(0, TMath::Abs(tempeplus));
419
420 /* prepare output array */
421 TObjArray *outoa = new TObjArray();
422 for (Int_t idata = 0; idata < nBW; idata++) {
423 outoa->Add(gBW[idata]);
424 outoa->Add(fBGBW[idata]);
425 }
426 outoa->Add(cBW);
427 outoa->Add(hBW);
428 outoa->Add(gBetaT);
429 if (gCont1) outoa->Add(gCont1);
430 if (gCont2) outoa->Add(gCont2);
431
432 return outoa;
433
434}
435
436void
437BGBlastWave_FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
438{
439
440 /* beta -> beta_max */
441 Double_t beta = par[nBW+0];
442 Double_t T = par[nBW+1];
443 Double_t n = par[nBW+2];
444 Double_t beta_max = 0.5 * (2. + n) * beta;
445#if 0
446 /* check beta_max */
447 if (beta_max >= 1. || beta_max <= 0.) {
448 f = kMaxInt;
449 return;
450 }
451 /* check T */
452 if (T <= 0.) {
453 f = kMaxInt;
454 return;
455 }
456#endif
457
458 Double_t pt, pte, val, vale, func, pull, chi = 0;
459 /* loop over all the data */
460 for (Int_t iBW = 0; iBW < nBW; iBW++) {
461 /* set BGBW parameters */
462 fBGBW[iBW]->SetParameter(4, par[iBW]);
463 fBGBW[iBW]->SetParameter(1, beta_max);
464 fBGBW[iBW]->SetParameter(2, T);
465 fBGBW[iBW]->SetParameter(3, n);
466 /* loop over all the points */
467 for (Int_t ipt = 0; ipt < gBW[iBW]->GetN(); ipt++) {
468 pt = gBW[iBW]->GetX()[ipt];
469 pte = gBW[iBW]->GetEX()[ipt];
470 val = gBW[iBW]->GetY()[ipt];
471 vale = gBW[iBW]->GetEY()[ipt];
472 func = fBGBW[iBW]->Eval(pt);
473 // func = fBGBW[iBW]->Integral(pt - pte, pt + pte);
474 pull = (val - func) / vale;
475 chi += pull * pull;
476 }
477 }
478
479 f = chi;
480}
481
482/*****************************************************************/
483
484GetYieldAndMean(TH1 *h, TF1 *f, Double_t &yield, Double_t &yielderr, Double_t &mean, Double_t &meanerr, Double_t min, Double_t max, Double_t *partyield, Double_t *partyielderr)
485{
486
487 /* find lowest edge in histo */
488 Int_t binlo;
489 Double_t lo;
490 for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
491 if (h->GetBinContent(ibin) != 0.) {
492 binlo = ibin;
493 lo = h->GetBinLowEdge(ibin);
494 break;
495 }
496 }
497
498 /* find highest edge in histo */
499 Int_t binhi;
500 Double_t hi;
501 for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
502 if (h->GetBinContent(ibin) != 0.) {
503 binhi = ibin + 1;
504 hi = h->GetBinLowEdge(ibin + 1);
505 break;
506 }
507 }
508
509 /* integrate the data */
510 Double_t cont, err, width, cent, integral_data = 0., integralerr_data = 0., meanintegral_data = 0., meanintegralerr_data = 0.;
511 for (Int_t ibin = binlo; ibin < binhi; ibin++) {
512 cent = h->GetBinCenter(ibin);
513 width = h->GetBinWidth(ibin);
514 cont = h->GetBinContent(ibin);
515 err = h->GetBinError(ibin);
516 /* check we didn't get an empty bin in between */
517 if (cont != 0. && err != 0.) {
518 /* all right, use data */
519 integral_data += cont * width;
520 integralerr_data += err * err * width * width;
521 meanintegral_data += cont * width * cent;
522 meanintegralerr_data += err * err * width * width * cent * cent;
523 }
524 else {
525 /* missing data-point, complain and use function */
526 printf("WARNING: missing data-point at %f\n", cent);
527 printf(" using function as a patch\n");
528 integral_data += f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
529 integralerr_data += f->IntegralError(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1), 0, 0, 1.e-6);
530 meanintegral_data += f->Mean(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1)) * f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
531 meanintegralerr_data += f->Mean(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1)) * f->IntegralError(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1), 0, 0, 1.e-6);
532 }
533 }
534 integralerr_data = TMath::Sqrt(integralerr_data);
535 meanintegralerr_data = TMath::Sqrt(meanintegralerr_data);
536
537 /* integrate below the data */
538 Double_t integral_lo = min < lo ? f->Integral(min, lo) : 0.;
539 Double_t integralerr_lo = min < lo ? f->IntegralError(min, lo, 0, 0, 1.e-6) : 0.;
540 Double_t meanintegral_lo = min < lo ? f->Mean(min, lo) * integral_lo : 0.;
541 Double_t meanintegralerr_lo = min < lo ? f->Mean(min, lo) * integralerr_lo : 0.;
542
543 /* integrate above the data */
544 Double_t integral_hi = max > hi ? f->Integral(hi, max) : 0.;
545 Double_t integralerr_hi = max > hi ? f->IntegralError(hi, max, 0, 0, 1.e-6) : 0.;
546 Double_t meanintegral_hi = max > hi ? f->Mean(hi, max) * integral_hi : 0.;
547 Double_t meanintegralerr_hi = max > hi ? f->Mean(hi, max) * integralerr_hi : 0.;
548
549 /* compute integrated yield */
550 yield = integral_data + integral_lo + integral_hi;
551 yielderr = TMath::Sqrt(integralerr_data * integralerr_data +
552 integralerr_lo * integralerr_lo +
553 integralerr_hi * integralerr_hi);
554
555 /* compute integrated mean */
556 mean = (meanintegral_data + meanintegral_lo + meanintegral_hi) / yield;
557 meanerr = TMath::Sqrt(meanintegralerr_data * meanintegralerr_data +
558 meanintegralerr_lo * meanintegralerr_lo +
559 meanintegralerr_hi * meanintegralerr_hi) / yield;
560
561 /* set partial yields */
562 partyield[0] = integral_data;
563 partyielderr[0] = integralerr_data;
564 partyield[1] = integral_lo;
565 partyielderr[1] = integralerr_lo;
566 partyield[2] = integral_hi;
567 partyielderr[2] = integralerr_hi;
568
569}
570
571/*****************************************************************/
572
573Double_t
574y2eta(Double_t pt, Double_t mass, Double_t y){
575 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
576 return TMath::ASinH(mt / pt * TMath::SinH(y));
577}
578Double_t
579eta2y(Double_t pt, Double_t mass, Double_t eta){
580 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
581 return TMath::ASinH(pt / mt * TMath::SinH(eta));
582}
583
584TH1 *
585Convert_dNdy_1over2pipt_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
586{
587
588 TH1 *hout = hin->Clone("hout");
589 hout->Reset();
590 Double_t pt, mt, conv, val, vale;
591 for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
592 pt = hin->GetBinCenter(ibin + 1);
593 conv = eta2y(pt, mass, eta) / eta;
594 val = hin->GetBinContent(ibin + 1);
595 vale = hin->GetBinError(ibin + 1);
596 val /= (2. * TMath::Pi() * pt);
597 vale /= (2. * TMath::Pi() * pt);
598 val *= conv;
599 vale *= conv;
600 hout->SetBinContent(ibin + 1, val);
601 hout->SetBinError(ibin + 1, vale);
602 }
603 return hout;
604}
605
606TH1 *
607Convert_dNdy_1over2pipt_dNdy(TH1 *hin)
608{
609
610 TH1 *hout = hin->Clone("hout");
611 hout->Reset();
612 Double_t pt, mt, conv, val, vale;
613 for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
614 pt = hin->GetBinCenter(ibin + 1);
615 val = hin->GetBinContent(ibin + 1);
616 vale = hin->GetBinError(ibin + 1);
617 val /= (2. * TMath::Pi() * pt);
618 vale /= (2. * TMath::Pi() * pt);
619 hout->SetBinContent(ibin + 1, val);
620 hout->SetBinError(ibin + 1, vale);
621 }
622 return hout;
623}
624
625TH1 *
626Convert_dNdy_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
627{
628
629 TH1 *hout = hin->Clone("hout");
630 hout->Reset();
631 Double_t pt, mt, conv, val, vale;
632 for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
633 pt = hin->GetBinCenter(ibin + 1);
634 conv = eta2y(pt, mass, eta) / eta;
635 val = hin->GetBinContent(ibin + 1);
636 vale = hin->GetBinError(ibin + 1);
637 val *= conv;
638 vale *= conv;
639 hout->SetBinContent(ibin + 1, val);
640 hout->SetBinError(ibin + 1, vale);
641 }
642 return hout;
643}
644
645TGraph *
646Convert_dNdy_dNdeta(TGraph *hin, Double_t mass, Double_t eta = 0.8)
647{
648
649 TGraph *hout = hin->Clone("hout");
650 // hout->Reset();
651 Double_t pt, mt, conv, val, valelo, valehi;
652 for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
653 pt = hin->GetX()[ibin];
654 conv = eta2y(pt, mass, eta) / eta;
655 val = hin->GetY()[ibin];
656 valelo = hin->GetEYlow()[ibin];
657 valehi = hin->GetEYhigh()[ibin];
658 val *= conv;
659 valelo *= conv;
660 valehi *= conv;
661 hout->GetX()[ibin] = pt;
662 hout->GetY()[ibin] = val;
663 hout->GetEYlow()[ibin] = valelo;
664 hout->GetEYhigh()[ibin] = valehi;
665 }
666 return hout;
667}
668
669TH1 *
670SummedId_1over2pipt_dNdeta(const Char_t *filename, Int_t icent)
671{
672
673 const Char_t *chargeName[2] = {
674 "plus", "minus"
675 };
676
677 TFile *filein = TFile::Open(filename);
678 TH1 *hy[AliPID::kSPECIES][2];
679 TH1 *heta[AliPID::kSPECIES][2];
680 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
681 for (Int_t icharge = 0; icharge < 2; icharge++) {
682 hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
683 if (!hy[ipart][icharge]) {
684 printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
685 return NULL;
686 }
687 heta[ipart][icharge] = Convert_dNdy_1over2pipt_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
688 }
689
690 /* sum */
691 TH1D *hsum = heta[2][0]->Clone("hsum");
692 hsum->Reset();
693 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
694 for (Int_t icharge = 0; icharge < 2; icharge++)
695 hsum->Add(heta[ipart][icharge]);
696
697 return hsum;
698}
699
700TH1 *
701SummedId_dNdeta(const Char_t *filename, Int_t icent)
702{
703
704 const Char_t *chargeName[2] = {
705 "plus", "minus"
706 };
707
708 TFile *filein = TFile::Open(filename);
709 TH1 *hy[AliPID::kSPECIES][2];
710 TH1 *heta[AliPID::kSPECIES][2];
711 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
712 for (Int_t icharge = 0; icharge < 2; icharge++) {
713 hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
714 if (!hy[ipart][icharge]) {
715 printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
716 return NULL;
717 }
718 heta[ipart][icharge] = Convert_dNdy_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
719 }
720
721 /* sum */
722 TH1D *hsum = heta[2][0]->Clone("hsum");
723 hsum->Reset();
724 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
725 for (Int_t icharge = 0; icharge < 2; icharge++)
726 hsum->Add(heta[ipart][icharge]);
727
728 return hsum;
729}
730
731/*****************************************************************/
732