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_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
 625 TH1 *
 626 Convert_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
 645 TGraph *
 646 Convert_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
 669 TH1 *
 670 SummedId_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
 700 TH1 *
 701 SummedId_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