]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/TOF/pPb502/macros/SpectraUtils.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TOF / pPb502 / macros / SpectraUtils.C
CommitLineData
59e49925 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
132Double_t
133BGBlastWave_Func_OneOverPt(const Double_t *x, const Double_t *p)
134{
135 /* 1/pt dN/dpt */
136
137 Double_t pt = x[0];
138 Double_t mass = p[0];
139 Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
140 Double_t beta_max = p[1];
141 Double_t temp = p[2];
142 Double_t n = p[3];
143 Double_t norm = p[4];
144
145 if (!fBGBlastWave_Integrand)
146 fBGBlastWave_Integrand = new TF1("fBGBlastWave_Integrand", BGBlastWave_Integrand, 0., 1., 5);
147 fBGBlastWave_Integrand->SetParameters(mt, pt, beta_max, temp, n);
148 Double_t integral = fBGBlastWave_Integrand->Integral(0., 1.);
149
150 return norm * integral;
151}
152
153
154TF1 *
155BGBlastWave(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)
156{
157
158 TF1 *fBGBlastWave = new TF1(name, BGBlastWave_Func, 0., 10., 5);
159 fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm);
160 fBGBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
161 fBGBlastWave->FixParameter(0, mass);
162 fBGBlastWave->SetParLimits(1, 0.01, 0.99);
163 fBGBlastWave->SetParLimits(2, 0.01, 1.);
164 fBGBlastWave->SetParLimits(3, 0.01, 10.);
165 return fBGBlastWave;
166}
167
168TF1 * BGBlastWave_OneOverPT(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)
169{
170
171 TF1 *fBGBlastWave = new TF1(name, BGBlastWave_Func_OneOverPt, 0., 10., 5);
172 fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm);
173 fBGBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
174 fBGBlastWave->FixParameter(0, mass);
175 fBGBlastWave->SetParLimits(1, 0.01, 0.99);
176 fBGBlastWave->SetParLimits(2, 0.01, 1.);
177 fBGBlastWave->SetParLimits(3, 0.01, 10.);
178 return fBGBlastWave;
179}
180
181/*****************************************************************/
182/* TSALLIS BLAST-WAVE */
183/*****************************************************************/
184
185static TF1 *fTsallisBlastWave_Integrand_r = NULL;
186Double_t
187TsallisBlastWave_Integrand_r(const Double_t *x, const Double_t *p)
188{
189 /*
190 x[0] -> r (radius)
191 p[0] -> mT (transverse mass)
192 p[1] -> pT (transverse momentum)
193 p[2] -> beta_max (surface velocity)
194 p[3] -> T (freezout temperature)
195 p[4] -> n (velocity profile)
196 p[5] -> q
197 p[6] -> y (rapidity)
198 p[7] -> phi (azimuthal angle)
199 */
200
201 Double_t r = x[0];
202 Double_t mt = p[0];
203 Double_t pt = p[1];
204 Double_t beta_max = p[2];
205 Double_t temp_1 = 1. / p[3];
206 Double_t n = p[4];
207 Double_t q = p[5];
208 Double_t y = p[6];
209 Double_t phi = p[7];
210
211 if (q <= 1.) return r;
212
213 Double_t beta = beta_max * TMath::Power(r, n);
214 Double_t rho = TMath::ATanH(beta);
215
216 Double_t part1 = mt * TMath::CosH(y) * TMath::CosH(rho);
217 Double_t part2 = pt * TMath::SinH(rho) * TMath::Cos(phi);
218 Double_t part3 = part1 - part2;
219 Double_t part4 = 1 + (q - 1.) * temp_1 * part3;
220 Double_t expo = -1. / (q - 1.);
221 // printf("part1=%f, part2=%f, part3=%f, part4=%f, expo=%f\n", part1, part2, part3, part4, expo);
222 Double_t part5 = TMath::Power(part4, expo);
223
224 return r * part5;
225}
226
227static TF1 *fTsallisBlastWave_Integrand_phi = NULL;
228Double_t
229TsallisBlastWave_Integrand_phi(const Double_t *x, const Double_t *p)
230{
231 /*
232 x[0] -> phi (azimuthal angle)
233 */
234
235 Double_t phi = x[0];
236 fTsallisBlastWave_Integrand_r->SetParameter(7, phi);
237 Double_t integral = fTsallisBlastWave_Integrand_r->Integral(0., 1.);
238 return integral;
239}
240
241static TF1 *fTsallisBlastWave_Integrand_y = NULL;
242Double_t
243TsallisBlastWave_Integrand_y(const Double_t *x, const Double_t *p)
244{
245 /*
246 x[0] -> y (rapidity)
247 */
248
249 Double_t y = x[0];
250 fTsallisBlastWave_Integrand_r->SetParameter(6, y);
251 Double_t integral = fTsallisBlastWave_Integrand_phi->Integral(-TMath::Pi(), TMath::Pi());
252 return TMath::CosH(y) * integral;
253}
254
255Double_t
256TsallisBlastWave_Func(const Double_t *x, const Double_t *p)
257{
258 /* dN/dpt */
259
260 Double_t pt = x[0];
261 Double_t mass = p[0];
262 Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
263 Double_t beta_max = p[1];
264 Double_t temp = p[2];
265 Double_t n = p[3];
266 Double_t q = p[4];
267 Double_t norm = p[5];
268
269 if (!fTsallisBlastWave_Integrand_r)
270 fTsallisBlastWave_Integrand_r = new TF1("fTsallisBlastWave_Integrand_r", TsallisBlastWave_Integrand_r, 0., 1., 8);
271 if (!fTsallisBlastWave_Integrand_phi)
272 fTsallisBlastWave_Integrand_phi = new TF1("fTsallisBlastWave_Integrand_phi", TsallisBlastWave_Integrand_phi, -TMath::Pi(), TMath::Pi(), 0);
273 if (!fTsallisBlastWave_Integrand_y)
274 fTsallisBlastWave_Integrand_y = new TF1("fTsallisBlastWave_Integrand_y", TsallisBlastWave_Integrand_y, -0.5, 0.5, 0);
275
276 fTsallisBlastWave_Integrand_r->SetParameters(mt, pt, beta_max, temp, n, q, 0., 0.);
277 Double_t integral = fTsallisBlastWave_Integrand_y->Integral(-0.5, 0.5);
278 return norm * pt * integral;
279}
280
281TF1 *
282TsallisBlastWave(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)
283{
284
285 TF1 *fTsallisBlastWave = new TF1(name, TsallisBlastWave_Func, 0., 10., 6);
286 fTsallisBlastWave->SetParameters(mass, beta_max, temp, n, q, norm);
287 fTsallisBlastWave->SetParNames("mass", "beta_max", "T", "n", "q", "norm");
288 fTsallisBlastWave->FixParameter(0, mass);
289 fTsallisBlastWave->SetParLimits(1, 0.01, 0.99);
290 fTsallisBlastWave->SetParLimits(2, 0.01, 1.);
291 fTsallisBlastWave->SetParLimits(3, 0.1, 10.);
292 fTsallisBlastWave->SetParLimits(4, 1., 10.);
293 return fTsallisBlastWave;
294}
295
296/*****************************************************************/
297/*****************************************************************/
298/*****************************************************************/
299
300
301TF1 *
302BGBlastWave_SingleFit(TH1 *h, Double_t mass, Option_t *opt = "")
303{
304
305 TF1 *f = BGBlastWave(Form("fBGBW_%s", h->GetName()), mass);
306 h->Fit(f);
307 h->Fit(f);
308 h->Fit(f, opt);
309 return f;
310
311}
312
313Int_t nBW;
314TF1 *fBGBW[1000];
315TGraphErrors *gBW[1000];
316
317TObjArray *
318BGBlastWave_GlobalFit(TObjArray *data, Double_t *mass, Double_t profile = .7, Bool_t fixProfile = kFALSE)
319{
320
321 /* get data */
322 nBW = data->GetEntries();
323 for (Int_t idata = 0; idata < nBW; idata++) {
324 gBW[idata] = (TGraphErrors *)data->At(idata);
325 gBW[idata]->SetName(Form("gBW%d", idata));
326 }
327
328 /* init BG blast-wave functions */
329 for (Int_t idata = 0; idata < nBW; idata++) {
330 printf("init BG-BlastWave function #%d: mass = %f\n", idata, mass[idata]);
331 fBGBW[idata] = BGBlastWave(Form("fBGBW%d", idata), mass[idata]);
332 }
333
334 /* display data */
335 TCanvas *cBW = new TCanvas("cBW");
336 cBW->Divide(nBW, 1);
337 for (Int_t idata = 0; idata < nBW; idata++) {
338 cBW->cd(idata + 1);
339 gBW[idata]->Draw("ap*");
340 }
341 cBW->Update();
342
343 /* init minuit: nBW normalizations + 3 (beta, T, n) BG-BlastWave params */
344 const Int_t nbwpars = 3;
345 const Int_t nfitpars = nBW + nbwpars;
346 TMinuit *minuit = new TMinuit(nfitpars);
347 minuit->SetFCN(BGBlastWave_FCN);
348 Double_t arglist[10];
349 Int_t ierflg = 0;
350 arglist[0] = 1;
351 minuit->mnexcm("SET ERR", arglist, 1, ierflg);
352 for (Int_t idata = 0; idata < nBW; idata++)
353 minuit->mnparm(idata, Form("norm%d", idata), 1.e6, 1., 0., 0., ierflg);
354 minuit->mnparm(nBW + 0, "<beta>", 0.65, 0.01, 0., 1., ierflg);
355 minuit->mnparm(nBW + 1, "T", 0.1, 0.01, 0., 1., ierflg);
356 minuit->mnparm(nBW + 2, "n", profile, 0.1, 0., 10., ierflg);
357 if (fixProfile) minuit->FixParameter(nBW + 2);
358
359 /* set strategy */
360 arglist[0] = 1;
361 minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
362
363 /* start MIGRAD minimization */
364 arglist[0] = 500000;
365 arglist[1] = 1.;
366 minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
367
368 /* set strategy */
369 arglist[0] = 2;
370 minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
371
372 /* start MIGRAD minimization */
373 arglist[0] = 500000;
374 arglist[1] = 1.;
375 minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
376
377 /* start IMPROVE minimization */
378 arglist[0] = 500000;
379 minuit->mnexcm("IMPROVE", arglist, 1, ierflg);
380
381 /* start MINOS */
382 arglist[0] = 500000;
383 arglist[1] = nBW + 1;
384 arglist[2] = nBW + 2;
385 arglist[3] = nBW + 3;
386 minuit->mnexcm("MINOS", arglist, 4, ierflg);
387
388 /* print results */
389 Double_t amin,edm,errdef;
390 Int_t nvpar,nparx,icstat;
391 minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
392 minuit->mnprin(4, amin);
393
394 /* get parameters */
395 Double_t beta, betae, betaeplus, betaeminus, betagcc, temp, tempe, tempeplus, tempeminus, tempgcc, prof, profe, profeplus, profeminus, profgcc;
396 minuit->GetParameter(nBW + 0, beta, betae);
397 minuit->mnerrs(nBW + 0, betaeplus, betaeminus, betae, betagcc);
398 minuit->GetParameter(nBW + 1, temp, tempe);
399 minuit->mnerrs(nBW + 1, tempeplus, tempeminus, tempe, tempgcc);
400 minuit->GetParameter(nBW + 2, prof, profe);
401 minuit->mnerrs(nBW + 2, profeplus, profeminus, profe, profgcc);
402 Double_t beta_max = 0.5 * (2. + prof) * beta;
403 Double_t norm[1000], norme[1000];
404 for (Int_t idata = 0; idata < nBW; idata++)
405 minuit->GetParameter(idata, norm[idata], norme[idata]);
406
407 /* printout */
408 printf("*********************************\n");
409 printf("beta_max = %f\n", beta_max);
410 printf("<beta> = %f +- %f (e+ = %f, e- = %f)\n", beta, betae, betaeplus, betaeminus);
411 printf("T = %f +- %f (e+ = %f, e- = %f)\n", temp, tempe, tempeplus, tempeminus);
412 printf("n = %f +- %f (e+ = %f, e- = %f)\n", prof, profe, profeplus, profeminus);
413
414 /* 1-sigma contour */
415 minuit->SetErrorDef(1);
416 TGraph *gCont1 = NULL;
417 gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
418 if (gCont1) gCont1->SetName("gCont1");
419
420 /* 2-sigma contour */
421 minuit->SetErrorDef(4);
422 TGraph *gCont2 = NULL;
423 // gCont2 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
424 if (gCont2) gCont2->SetName("gCont2");
425
426 /* display fits */
427 for (Int_t idata = 0; idata < nBW; idata++) {
428 cBW->cd(idata + 1);
429 fBGBW[idata]->SetParameter(4, norm[idata]);
430 fBGBW[idata]->SetParameter(1, beta_max);
431 fBGBW[idata]->SetParameter(2, temp);
432 fBGBW[idata]->SetParameter(3, prof);
433 fBGBW[idata]->Draw("same");
434 }
435 cBW->Update();
436
437 /* histo params */
438 TH1D *hBW = new TH1D("hBW", "", 3, 0., 3.);
439 hBW->SetBinContent(1, beta);
440 hBW->SetBinError(1, betae);
441 hBW->SetBinContent(2, temp);
442 hBW->SetBinError(2, tempe);
443 hBW->SetBinContent(3, prof);
444 hBW->SetBinError(3, profe);
445
446 /* BW graph */
447 TGraphAsymmErrors *gBetaT = new TGraphAsymmErrors();
448 gBetaT->SetName("gBetaT");
449 gBetaT->SetPoint(0, beta, temp);
450 gBetaT->SetPointEXlow(0, TMath::Abs(betaeminus));
451 gBetaT->SetPointEXhigh(0, TMath::Abs(betaeplus));
452 gBetaT->SetPointEYlow(0, TMath::Abs(tempeminus));
453 gBetaT->SetPointEYhigh(0, TMath::Abs(tempeplus));
454
455 /* prepare output array */
456 TObjArray *outoa = new TObjArray();
457 for (Int_t idata = 0; idata < nBW; idata++) {
458 outoa->Add(gBW[idata]);
459 outoa->Add(fBGBW[idata]);
460 }
461 outoa->Add(cBW);
462 outoa->Add(hBW);
463 outoa->Add(gBetaT);
464 if (gCont1) outoa->Add(gCont1);
465 if (gCont2) outoa->Add(gCont2);
466
467 return outoa;
468
469}
470
471void
472BGBlastWave_FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
473{
474
475 /* beta -> beta_max */
476 Double_t beta = par[nBW+0];
477 Double_t T = par[nBW+1];
478 Double_t n = par[nBW+2];
479 Double_t beta_max = 0.5 * (2. + n) * beta;
480#if 0
481 /* check beta_max */
482 if (beta_max >= 1. || beta_max <= 0.) {
483 f = kMaxInt;
484 return;
485 }
486 /* check T */
487 if (T <= 0.) {
488 f = kMaxInt;
489 return;
490 }
491#endif
492
493 Double_t pt, pte, val, vale, func, pull, chi = 0;
494 /* loop over all the data */
495 for (Int_t iBW = 0; iBW < nBW; iBW++) {
496 /* set BGBW parameters */
497 fBGBW[iBW]->SetParameter(4, par[iBW]);
498 fBGBW[iBW]->SetParameter(1, beta_max);
499 fBGBW[iBW]->SetParameter(2, T);
500 fBGBW[iBW]->SetParameter(3, n);
501 /* loop over all the points */
502 for (Int_t ipt = 0; ipt < gBW[iBW]->GetN(); ipt++) {
503 pt = gBW[iBW]->GetX()[ipt];
504 pte = gBW[iBW]->GetEX()[ipt];
505 val = gBW[iBW]->GetY()[ipt];
506 vale = gBW[iBW]->GetEY()[ipt];
507 func = fBGBW[iBW]->Eval(pt);
508 // func = fBGBW[iBW]->Integral(pt - pte, pt + pte);
509 pull = (val - func) / vale;
510 chi += pull * pull;
511 }
512 }
513
514 f = chi;
515}
516
517/*****************************************************************/
518
519IntegratedProduction(TH1 *h, Double_t mass, Option_t *opt = "")
520{
521
522 Double_t yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, partyield[3], partyielderr[3], partyielderrcorr[3];
523 TF1 *f = BGBlastWave_SingleFit(h, mass, opt);
524 GetYieldAndMean(h, f, yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, 0., 10., partyield, partyielderr, partyielderrcorr);
525
526 // Double_t fint = f->Integral(0.,10.);
527 // Double_t finte = f->IntegralError(0.,10.);
528 // Double_t fmean = f->Mean(0., 10.);
529
530 printf("dN/dy = %f +- %f (%f)\n", yield, yielderr, yielderrcorr);
531 printf("<pt> = %f +- %f (%f)\n", mean, meanerr, meanerrcorr);
532 printf("dN/dy (data) = %f +- %f (%f)\n", partyield[0], partyielderr[0], partyielderrcorr[0]);
533 printf("dN/dy (low) = %f +- %f (%f)\n", partyield[1], partyielderr[1], partyielderrcorr[1]);
534 printf("dN/dy (high) = %f +- %f (%f)\n", partyield[2], partyielderr[2], partyielderrcorr[2]);
535 // printf("----\n");
536 // printf("dN/dy (func) = %f +- %f\n", fint, finte);
537 // printf("<pT> (func) = %f +- %f\n", fmean, 0.);
538
539 // TH1 *hr = (TH1 *)h->Clone("hr");
540 // hr->Divide(f);
541 // new TCanvas;
542 // hr->Draw();
543
544 // TProfile *p = new TProfile("p", "", 100, 0., 10.);
545 // gROOT->LoadMacro("HistoUtils.C");
546 // HistoUtils_Function2Profile(f, p);
547 // p->Draw();
548}
549
550GetYieldAndMean(TH1 *h, TF1 *f, Double_t &yield, Double_t &yielderr, Double_t &yielderrcorr, Double_t &mean, Double_t &meanerr, Double_t &meanerrcorr, Double_t min, Double_t max, Double_t *partyield, Double_t *partyielderr, Double_t *partyielderrcorr)
551{
552
553 /* find lowest edge in histo */
554 Int_t binlo;
555 Double_t lo;
556 for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
557 if (h->GetBinContent(ibin) != 0.) {
558 binlo = ibin;
559 lo = h->GetBinLowEdge(ibin);
560 break;
561 }
562 }
563
564 /* find highest edge in histo */
565 Int_t binhi;
566 Double_t hi;
567 for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
568 if (h->GetBinContent(ibin) != 0.) {
569 binhi = ibin + 1;
570 hi = h->GetBinLowEdge(ibin + 1);
571 break;
572 }
573 }
574
575 /* integrate the data */
576 Double_t cont, err, width, cent, integral_data = 0., integralerr_data = 0., integralerrcorr_data = 0., meanintegral_data = 0., meanintegralerr_data = 0., meanintegralerrcorr_data = 0.;
577 for (Int_t ibin = binlo; ibin < binhi; ibin++) {
578 cent = h->GetBinCenter(ibin);
579 width = h->GetBinWidth(ibin);
580 cont = h->GetBinContent(ibin);
581 err = h->GetBinError(ibin);
582 /* check we didn't get an empty bin in between */
583 if (cont != 0. && err != 0.) {
584 /* all right, use data */
585 integral_data += cont * width;
586 integralerr_data += err * err * width * width;
587 integralerrcorr_data += err * width;
588 meanintegral_data += cont * width * cent;
589 meanintegralerr_data += err * err * width * width * cent * cent;
590 meanintegralerrcorr_data += err * width * cent;
591 }
592 else {
593 /* missing data-point, complain and use function */
594 printf("WARNING: missing data-point at %f\n", cent);
595 printf(" using function as a patch\n");
596 integral_data += f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
597 integralerr_data += f->IntegralError(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1), 0, 0, 1.e-6);
598 meanintegral_data += f->Mean(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1)) * f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
599 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);
600 }
601 }
602 integralerr_data = TMath::Sqrt(integralerr_data);
603 meanintegralerr_data = TMath::Sqrt(meanintegralerr_data);
604
605 /* integrate below the data */
606 Double_t integral_lo = min < lo ? f->Integral(min, lo) : 0.;
607 Double_t integralerr_lo = min < lo ? f->IntegralError(min, lo, 0, 0, 1.e-6) : 0.;
608 Double_t meanintegral_lo = min < lo ? f->Mean(min, lo) * integral_lo : 0.;
609 Double_t meanintegralerr_lo = min < lo ? f->Mean(min, lo) * integralerr_lo : 0.;
610
611 /* integrate above the data */
612 Double_t integral_hi = max > hi ? f->Integral(hi, max) : 0.;
613 Double_t integralerr_hi = max > hi ? f->IntegralError(hi, max, 0, 0, 1.e-6) : 0.;
614 Double_t meanintegral_hi = max > hi ? f->Mean(hi, max) * integral_hi : 0.;
615 Double_t meanintegralerr_hi = max > hi ? f->Mean(hi, max) * integralerr_hi : 0.;
616
617 /* compute integrated yield */
618 yield = integral_data + integral_lo + integral_hi;
619 yielderr = TMath::Sqrt(integralerr_data * integralerr_data +
620 integralerr_lo * integralerr_lo +
621 integralerr_hi * integralerr_hi);
622 yielderrcorr = TMath::Sqrt(integralerrcorr_data * integralerrcorr_data +
623 integralerr_lo * integralerr_lo +
624 integralerr_hi * integralerr_hi);
625
626 /* compute integrated mean */
627 mean = (meanintegral_data + meanintegral_lo + meanintegral_hi) / yield;
628 meanerr = TMath::Sqrt(meanintegralerr_data * meanintegralerr_data +
629 meanintegralerr_lo * meanintegralerr_lo +
630 meanintegralerr_hi * meanintegralerr_hi) / yield;
631 meanerrcorr = TMath::Sqrt(meanintegralerrcorr_data * meanintegralerrcorr_data +
632 meanintegralerr_lo * meanintegralerr_lo +
633 meanintegralerr_hi * meanintegralerr_hi) / yield;
634
635 /* set partial yields */
636 partyield[0] = integral_data;
637 partyielderr[0] = integralerr_data;
638 partyielderrcorr[0] = integralerrcorr_data;
639 partyield[1] = integral_lo;
640 partyielderr[1] = integralerr_lo;
641 partyielderrcorr[1] = integralerr_lo;
642 partyield[2] = integral_hi;
643 partyielderr[2] = integralerr_hi;
644 partyielderrcorr[2] = integralerr_hi;
645
646}
647
648/*****************************************************************/
649
650Double_t
651y2eta(Double_t pt, Double_t mass, Double_t y){
652 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
653 return TMath::ASinH(mt / pt * TMath::SinH(y));
654}
655Double_t
656eta2y(Double_t pt, Double_t mass, Double_t eta){
657 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
658 return TMath::ASinH(pt / mt * TMath::SinH(eta));
659}
660
661TH1 *
662Convert_dNdy_1over2pipt_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
663{
664
665 TH1 *hout = hin->Clone("hout");
666 hout->Reset();
667 Double_t pt, mt, conv, val, vale;
668 for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
669 pt = hin->GetBinCenter(ibin + 1);
670 conv = eta2y(pt, mass, eta) / eta;
671 val = hin->GetBinContent(ibin + 1);
672 vale = hin->GetBinError(ibin + 1);
673 val /= (2. * TMath::Pi() * pt);
674 vale /= (2. * TMath::Pi() * pt);
675 val *= conv;
676 vale *= conv;
677 hout->SetBinContent(ibin + 1, val);
678 hout->SetBinError(ibin + 1, vale);
679 }
680 return hout;
681}
682
683TH1 *
684Convert_dNdy_1over2pipt_dNdy(TH1 *hin)
685{
686
687 TH1 *hout = hin->Clone("hout");
688 hout->Reset();
689 Double_t pt, mt, conv, val, vale;
690 for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
691 pt = hin->GetBinCenter(ibin + 1);
692 val = hin->GetBinContent(ibin + 1);
693 vale = hin->GetBinError(ibin + 1);
694 val /= (2. * TMath::Pi() * pt);
695 vale /= (2. * TMath::Pi() * pt);
696 hout->SetBinContent(ibin + 1, val);
697 hout->SetBinError(ibin + 1, vale);
698 }
699 return hout;
700}
701
702TH1 *
703Convert_dNdy_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
704{
705
706 TH1 *hout = hin->Clone("hout");
707 hout->Reset();
708 Double_t pt, mt, conv, val, vale;
709 for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
710 pt = hin->GetBinCenter(ibin + 1);
711 conv = eta2y(pt, mass, eta) / eta;
712 val = hin->GetBinContent(ibin + 1);
713 vale = hin->GetBinError(ibin + 1);
714 val *= conv;
715 vale *= conv;
716 hout->SetBinContent(ibin + 1, val);
717 hout->SetBinError(ibin + 1, vale);
718 }
719 return hout;
720}
721
722TGraph *
723Convert_dNdy_dNdeta(TGraph *hin, Double_t mass, Double_t eta = 0.8)
724{
725
726 TGraph *hout = hin->Clone("hout");
727 // hout->Reset();
728 Double_t pt, mt, conv, val, valelo, valehi;
729 for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
730 pt = hin->GetX()[ibin];
731 conv = eta2y(pt, mass, eta) / eta;
732 val = hin->GetY()[ibin];
733 valelo = hin->GetEYlow()[ibin];
734 valehi = hin->GetEYhigh()[ibin];
735 val *= conv;
736 valelo *= conv;
737 valehi *= conv;
738 hout->GetX()[ibin] = pt;
739 hout->GetY()[ibin] = val;
740 hout->GetEYlow()[ibin] = valelo;
741 hout->GetEYhigh()[ibin] = valehi;
742 }
743 return hout;
744}
745
746TH1 *
747SummedId_1over2pipt_dNdeta(const Char_t *filename, Int_t icent)
748{
749
750 const Char_t *chargeName[2] = {
751 "plus", "minus"
752 };
753
754 TFile *filein = TFile::Open(filename);
755 TH1 *hy[AliPID::kSPECIES][2];
756 TH1 *heta[AliPID::kSPECIES][2];
757 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
758 for (Int_t icharge = 0; icharge < 2; icharge++) {
759 hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
760 if (!hy[ipart][icharge]) {
761 printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
762 return NULL;
763 }
764 heta[ipart][icharge] = Convert_dNdy_1over2pipt_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
765 }
766
767 /* sum */
768 TH1D *hsum = heta[2][0]->Clone("hsum");
769 hsum->Reset();
770 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
771 for (Int_t icharge = 0; icharge < 2; icharge++)
772 hsum->Add(heta[ipart][icharge]);
773
774 return hsum;
775}
776
777TH1 *
778SummedId_dNdeta(const Char_t *filename, Int_t icent)
779{
780
781 const Char_t *chargeName[2] = {
782 "plus", "minus"
783 };
784
785 TFile *filein = TFile::Open(filename);
786 TH1 *hy[AliPID::kSPECIES][2];
787 TH1 *heta[AliPID::kSPECIES][2];
788 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
789 for (Int_t icharge = 0; icharge < 2; icharge++) {
790 hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
791 if (!hy[ipart][icharge]) {
792 printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
793 return NULL;
794 }
795 heta[ipart][icharge] = Convert_dNdy_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
796 }
797
798 /* sum */
799 TH1D *hsum = heta[2][0]->Clone("hsum");
800 hsum->Reset();
801 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
802 for (Int_t icharge = 0; icharge < 2; icharge++)
803 hsum->Add(heta[ipart][icharge]);
804
805 return hsum;
806}
807
808/*****************************************************************/
809