added dN/dy -> dN/deta conversion functions
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / COMBINED / SpectraUtils.C
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
15 Double_t
16 Boltzmann_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
29 TF1 *
30 Boltzmann(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
44 Double_t
45 LevyTsallis_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
64 TF1 *
65 LevyTsallis(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
79 static TF1 *fBGBlastWave_Integrand = NULL;
80 Double_t
81 BGBlastWave_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
112 Double_t
113 BGBlastWave_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
132 TF1 *
133 BGBlastWave(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
150 static TF1 *fTsallisBlastWave_Integrand_r = NULL;
151 Double_t
152 TsallisBlastWave_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
192 static TF1 *fTsallisBlastWave_Integrand_phi = NULL;
193 Double_t
194 TsallisBlastWave_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
206 static TF1 *fTsallisBlastWave_Integrand_y = NULL;
207 Double_t
208 TsallisBlastWave_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
220 Double_t
221 TsallisBlastWave_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
246 TF1 *
247 TsallisBlastWave(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
266 TF1 *
267 BGBlastWave_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
278 Int_t nBW;
279 TF1 *fBGBW[1000];
280 TGraphErrors *gBW[1000];
281
282 TObjArray *
283 BGBlastWave_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
436 void 
437 BGBlastWave_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
484 GetYieldAndMean(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
573 Double_t 
574 y2eta(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 }
578 Double_t 
579 eta2y(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
584 TH1 *
585 Convert_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
606 TH1 *
607 Convert_dNdy_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
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     conv = eta2y(pt, mass, eta) / eta;
616     val = hin->GetBinContent(ibin + 1);
617     vale = hin->GetBinError(ibin + 1);
618     val *= conv;
619     vale *= conv;
620     hout->SetBinContent(ibin + 1, val);
621     hout->SetBinError(ibin + 1, vale);
622   }
623   return hout;
624 }
625
626 TH1 *
627 SummedId_1over2pipt_dNdeta(const Char_t *filename, Int_t icent)
628 {
629
630   const Char_t *chargeName[2] = {
631     "plus", "minus"
632   };
633
634   TFile *filein = TFile::Open(filename);
635   TH1 *hy[AliPID::kSPECIES][2];
636   TH1 *heta[AliPID::kSPECIES][2];
637   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
638     for (Int_t icharge = 0; icharge < 2; icharge++) {
639       hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
640       if (!hy[ipart][icharge]) {
641         printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
642         return NULL;
643       }
644       heta[ipart][icharge] = Convert_dNdy_1over2pipt_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
645     }
646
647   /* sum */
648   TH1D *hsum = heta[2][0]->Clone("hsum");
649   hsum->Reset();
650   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
651     for (Int_t icharge = 0; icharge < 2; icharge++)
652       hsum->Add(heta[ipart][icharge]);
653
654   return hsum;
655 }
656
657 TH1 *
658 SummedId_dNdeta(const Char_t *filename, Int_t icent)
659 {
660
661   const Char_t *chargeName[2] = {
662     "plus", "minus"
663   };
664
665   TFile *filein = TFile::Open(filename);
666   TH1 *hy[AliPID::kSPECIES][2];
667   TH1 *heta[AliPID::kSPECIES][2];
668   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
669     for (Int_t icharge = 0; icharge < 2; icharge++) {
670       hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
671       if (!hy[ipart][icharge]) {
672         printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
673         return NULL;
674       }
675       heta[ipart][icharge] = Convert_dNdy_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
676     }
677
678   /* sum */
679   TH1D *hsum = heta[2][0]->Clone("hsum");
680   hsum->Reset();
681   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
682     for (Int_t icharge = 0; icharge < 2; icharge++)
683       hsum->Add(heta[ipart][icharge]);
684
685   return hsum;
686 }
687
688 /*****************************************************************/
689