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 Double_t
 133 BGBlastWave_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
 154 TF1 *
 155 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)
 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
 168 TF1 * 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
 185 static TF1 *fTsallisBlastWave_Integrand_r = NULL;
 186 Double_t
 187 TsallisBlastWave_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
 227 static TF1 *fTsallisBlastWave_Integrand_phi = NULL;
 228 Double_t
 229 TsallisBlastWave_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
 241 static TF1 *fTsallisBlastWave_Integrand_y = NULL;
 242 Double_t
 243 TsallisBlastWave_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
 255 Double_t
 256 TsallisBlastWave_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
 281 TF1 *
 282 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)
 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
 301 TF1 *
 302 BGBlastWave_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
 313 Int_t nBW;
 314 TF1 *fBGBW[1000];
 315 TGraphErrors *gBW[1000];
 316
 317 TObjArray *
 318 BGBlastWave_GlobalFit(TObjArray *data, Double_t *mass, Double_t profile = .9, 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 + 0, "<beta>", 0.55, 0.01, 0., 1., ierflg);
 357   minuit->mnparm(nBW + 1, "T", 0.13, 0.01, 0., 1., ierflg);
 358   minuit->mnparm(nBW + 2, "n", profile, 0.1, 0., 10., ierflg);
 359   if (fixProfile) minuit->FixParameter(nBW + 2);
 360
 361   /* set strategy */
 362   arglist[0] = 1;
 363   minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
 364
 365   /* start MIGRAD minimization */
 366   arglist[0] = 500000;
 367   arglist[1] = 1.;
 368   minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
 369
 370   /* set strategy */
 371   arglist[0] = 2;
 372   minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
 373
 374   /* start MIGRAD minimization */
 375   arglist[0] = 500000;
 376   arglist[1] = 1.;
 377   minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
 378
 379   /* start IMPROVE minimization */
 380   arglist[0] = 500000;
 381   minuit->mnexcm("IMPROVE", arglist, 1, ierflg);
 382
 383   /* start MINOS */
 384   arglist[0] = 500000;
 385   arglist[1] = nBW + 1;
 386   arglist[2] = nBW + 2;
 387   arglist[3] = nBW + 3;
 388   minuit->mnexcm("MINOS", arglist, 4, ierflg);
 389
 390   /* print results */
 391   Double_t amin,edm,errdef;
 392   Int_t nvpar,nparx,icstat;
 393   minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
 394   minuit->mnprin(4, amin);
 395
 396   /* get parameters */
 397   Double_t beta, betae, betaeplus, betaeminus, betagcc, temp, tempe, tempeplus, tempeminus, tempgcc, prof, profe, profeplus, profeminus, profgcc;
 398   minuit->GetParameter(nBW + 0, beta, betae);
 399   minuit->mnerrs(nBW + 0, betaeplus, betaeminus, betae, betagcc);
 400   minuit->GetParameter(nBW + 1, temp, tempe);
 401   minuit->mnerrs(nBW + 1, tempeplus, tempeminus, tempe, tempgcc);
 402   minuit->GetParameter(nBW + 2, prof, profe);
 403   minuit->mnerrs(nBW + 2, profeplus, profeminus, profe, profgcc);
 404   Double_t beta_max = 0.5 * (2. + prof) * beta;
 405   Double_t norm[1000], norme[1000];
 406   for (Int_t idata = 0; idata < nBW; idata++)
 407     minuit->GetParameter(idata, norm[idata], norme[idata]);
 408
 409   /* printout */
 410   printf("*********************************\n");
 411   printf("beta_max = %f\n", beta_max);
 412   printf("<beta>   = %f +- %f (e+ = %f, e- = %f)\n", beta, betae, betaeplus, betaeminus);
 413   printf("T        = %f +- %f (e+ = %f, e- = %f)\n", temp, tempe, tempeplus, tempeminus);
 414   printf("n        = %f +- %f (e+ = %f, e- = %f)\n", prof, profe, profeplus, profeminus);
 415
 416   /* 1-sigma contour */
 417   minuit->SetErrorDef(1);
 418   TGraph *gCont1 = NULL;
 419   //  gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
 420   if (gCont1) gCont1->SetName("gCont1");
 421
 422   /* 2-sigma contour */
 423   minuit->SetErrorDef(4);
 424   TGraph *gCont2 = NULL;
 425   //  gCont2 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
 426   if (gCont2) gCont2->SetName("gCont2");
 427
 428   /* display fits */
 429   for (Int_t idata = 0; idata < nBW; idata++) {
 430     cBW->cd(idata + 1);
 431     fBGBW[idata]->SetParameter(4, norm[idata]);
 432     fBGBW[idata]->SetParameter(1, beta_max);
 433     fBGBW[idata]->SetParameter(2, temp);
 434     fBGBW[idata]->SetParameter(3, prof);
 435     fBGBW[idata]->Draw("same");
 436   }
 437   cBW->Update();
 438
 439   /* histo params */
 440   TH1D *hBW = new TH1D("hBW", "", 3, 0., 3.);
 441   hBW->SetBinContent(1, beta);
 442   hBW->SetBinError(1, betae);
 443   hBW->SetBinContent(2, temp);
 444   hBW->SetBinError(2, tempe);
 445   hBW->SetBinContent(3, prof);
 446   hBW->SetBinError(3, profe);
 447
 448   /* BW graph */
 449   TGraphAsymmErrors *gBetaT = new TGraphAsymmErrors();
 450   gBetaT->SetName("gBetaT");
 451   gBetaT->SetPoint(0, beta, temp);
 452   gBetaT->SetPointEXlow(0, TMath::Abs(betaeminus));
 453   gBetaT->SetPointEXhigh(0, TMath::Abs(betaeplus));
 454   gBetaT->SetPointEYlow(0, TMath::Abs(tempeminus));
 455   gBetaT->SetPointEYhigh(0, TMath::Abs(tempeplus));
 456
 457   /* prepare output array */
 458   TObjArray *outoa = new TObjArray();
 459   for (Int_t idata = 0; idata < nBW; idata++) {
 460     outoa->Add(gBW[idata]);
 461     outoa->Add(fBGBW[idata]);
 462   }
 463   outoa->Add(cBW);
 464   outoa->Add(hBW);
 465   outoa->Add(gBetaT);
 466   if (gCont1) outoa->Add(gCont1);
 467   if (gCont2) outoa->Add(gCont2);
 468
 469   return outoa;
 470
 471 }
 472
 473 void
 474 BGBlastWave_FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
 475 {
 476
 477   /* beta -> beta_max */
 478   Double_t beta = par[nBW+0];
 479   Double_t T = par[nBW+1];
 480   Double_t n = par[nBW+2];
 481   Double_t beta_max = 0.5 * (2. + n) * beta;
 482 #if 0
 483   /* check beta_max */
 484   if (beta_max >= 1. || beta_max <= 0.) {
 485     f = kMaxInt;
 486     return;
 487   }
 488   /* check T */
 489   if (T <= 0.) {
 490     f = kMaxInt;
 491     return;
 492   }
 493 #endif
 494
 495   Double_t pt, pte, val, vale, func, pull, chi = 0;
 496   /* loop over all the data */
 497   for (Int_t iBW = 0; iBW < nBW; iBW++) {
 498     /* set BGBW parameters */
 499     fBGBW[iBW]->SetParameter(4, par[iBW]);
 500     fBGBW[iBW]->SetParameter(1, beta_max);
 501     fBGBW[iBW]->SetParameter(2, T);
 502     fBGBW[iBW]->SetParameter(3, n);
 503     /* loop over all the points */
 504     for (Int_t ipt = 0; ipt < gBW[iBW]->GetN(); ipt++) {
 505       pt = gBW[iBW]->GetX()[ipt];
 506       pte = gBW[iBW]->GetEX()[ipt];
 507       val = gBW[iBW]->GetY()[ipt];
 508       vale = gBW[iBW]->GetEY()[ipt];
 509       func = fBGBW[iBW]->Eval(pt);
 510       //      func = fBGBW[iBW]->Integral(pt - pte, pt + pte);
 511       pull = (val - func) / vale;
 512       chi += pull * pull;
 513     }
 514   }
 515
 516   f = chi;
 517 }
 518
 519 /*****************************************************************/
 520
 521 IntegratedProduction(TH1 *h, Double_t mass, Option_t *opt = "")
 522 {
 523
 524   Double_t yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, partyield[3], partyielderr[3], partyielderrcorr[3];
 525   TF1 *f = BGBlastWave_SingleFit(h, mass, opt);
 526   GetYieldAndMean(h, f, yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, 0., 10., partyield, partyielderr, partyielderrcorr);
 527
 528   //  Double_t fint = f->Integral(0.,10.);
 529   //  Double_t finte = f->IntegralError(0.,10.);
 530   //  Double_t fmean = f->Mean(0., 10.);
 531
 532   printf("dN/dy        = %f +- %f (%f)\n", yield, yielderr, yielderrcorr);
 533   printf("<pt>         = %f +- %f (%f)\n", mean, meanerr, meanerrcorr);
 534   printf("dN/dy (data) = %f +- %f (%f)\n", partyield[0], partyielderr[0], partyielderrcorr[0]);
 535   printf("dN/dy (low)  = %f +- %f (%f)\n", partyield[1], partyielderr[1], partyielderrcorr[1]);
 536   printf("dN/dy (high) = %f +- %f (%f)\n", partyield[2], partyielderr[2], partyielderrcorr[2]);
 537   //  printf("----\n");
 538   //  printf("dN/dy (func) = %f +- %f\n", fint, finte);
 539   //  printf("<pT> (func)  = %f +- %f\n", fmean, 0.);
 540
 541   //  TH1 *hr = (TH1 *)h->Clone("hr");
 542   //  hr->Divide(f);
 543   //  new TCanvas;
 544   //  hr->Draw();
 545
 546   //  TProfile *p = new TProfile("p", "", 100, 0., 10.);
 547   //  gROOT->LoadMacro("HistoUtils.C");
 548   //  HistoUtils_Function2Profile(f, p);
 549   //  p->Draw();
 550 }
 551
 552 GetYieldAndMean(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)
 553 {
 554
 555   /* find lowest edge in histo */
 556   Int_t binlo;
 557   Double_t lo;
 558   for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
 559     if (h->GetBinContent(ibin) != 0.) {
 560       binlo = ibin;
 561       lo = h->GetBinLowEdge(ibin);
 562       break;
 563     }
 564   }
 565
 566   /* find highest edge in histo */
 567   Int_t binhi;
 568   Double_t hi;
 569   for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
 570     if (h->GetBinContent(ibin) != 0.) {
 571       binhi = ibin + 1;
 572       hi = h->GetBinLowEdge(ibin + 1);
 573       break;
 574     }
 575   }
 576
 577   /* integrate the data */
 578   Double_t cont, err, width, cent, integral_data = 0., integralerr_data = 0., integralerrcorr_data = 0., meanintegral_data = 0., meanintegralerr_data = 0., meanintegralerrcorr_data = 0.;
 579   for (Int_t ibin = binlo; ibin < binhi; ibin++) {
 580     cent = h->GetBinCenter(ibin);
 581     width = h->GetBinWidth(ibin);
 582     cont = h->GetBinContent(ibin);
 583     err = h->GetBinError(ibin);
 584     /* check we didn't get an empty bin in between */
 585     if (cont != 0. && err != 0.) {
 586       /* all right, use data */
 587       integral_data += cont * width;
 588       integralerr_data += err * err * width * width;
 589       integralerrcorr_data += err * width;
 590       meanintegral_data += cont * width * cent;
 591       meanintegralerr_data += err * err * width * width * cent * cent;
 592       meanintegralerrcorr_data += err * width * cent;
 593     }
 594     else {
 595       /* missing data-point, complain and use function */
 596       printf("WARNING: missing data-point at %f\n", cent);
 597       printf("         using function as a patch\n");
 598       integral_data += f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
 599       integralerr_data += f->IntegralError(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1), 0, 0, 1.e-6);
 600       meanintegral_data += f->Mean(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1)) * f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
 601       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);
 602     }
 603   }
 604   integralerr_data = TMath::Sqrt(integralerr_data);
 605   meanintegralerr_data = TMath::Sqrt(meanintegralerr_data);
 606
 607   /* integrate below the data */
 608   Double_t integral_lo = min < lo ? f->Integral(min, lo) : 0.;
 609   Double_t integralerr_lo = min < lo ? f->IntegralError(min, lo, 0, 0, 1.e-6) : 0.;
 610   Double_t meanintegral_lo = min < lo ? f->Mean(min, lo) * integral_lo : 0.;
 611   Double_t meanintegralerr_lo = min < lo ? f->Mean(min, lo) * integralerr_lo : 0.;
 612
 613   /* integrate above the data */
 614   Double_t integral_hi = max > hi ? f->Integral(hi, max) : 0.;
 615   Double_t integralerr_hi = max > hi ? f->IntegralError(hi, max, 0, 0, 1.e-6) : 0.;
 616   Double_t meanintegral_hi = max > hi ? f->Mean(hi, max) * integral_hi : 0.;
 617   Double_t meanintegralerr_hi = max > hi ? f->Mean(hi, max) * integralerr_hi : 0.;
 618
 619   /* compute integrated yield */
 620   yield = integral_data + integral_lo + integral_hi;
 621   yielderr = TMath::Sqrt(integralerr_data * integralerr_data +
 622                          integralerr_lo * integralerr_lo +
 623                          integralerr_hi * integralerr_hi);
 624   yielderrcorr = TMath::Sqrt(integralerrcorr_data * integralerrcorr_data +
 625                              integralerr_lo * integralerr_lo +
 626                              integralerr_hi * integralerr_hi);
 627
 628   /* compute integrated mean */
 629   mean = (meanintegral_data + meanintegral_lo + meanintegral_hi) / yield;
 630   meanerr = TMath::Sqrt(meanintegralerr_data * meanintegralerr_data +
 631                         meanintegralerr_lo * meanintegralerr_lo +
 632                         meanintegralerr_hi * meanintegralerr_hi) / yield;
 633   meanerrcorr = TMath::Sqrt(meanintegralerrcorr_data * meanintegralerrcorr_data +
 634                             meanintegralerr_lo * meanintegralerr_lo +
 635                             meanintegralerr_hi * meanintegralerr_hi) / yield;
 636
 637   /* set partial yields */
 638   partyield[0] = integral_data;
 639   partyielderr[0] = integralerr_data;
 640   partyielderrcorr[0] = integralerrcorr_data;
 641   partyield[1] = integral_lo;
 642   partyielderr[1] = integralerr_lo;
 643   partyielderrcorr[1] = integralerr_lo;
 644   partyield[2] = integral_hi;
 645   partyielderr[2] = integralerr_hi;
 646   partyielderrcorr[2] = integralerr_hi;
 647
 648 }
 649
 650 /*****************************************************************/
 651
 652 Double_t
 653 y2eta(Double_t pt, Double_t mass, Double_t y){
 654   Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
 655   return TMath::ASinH(mt / pt * TMath::SinH(y));
 656 }
 657 Double_t
 658 eta2y(Double_t pt, Double_t mass, Double_t eta){
 659   Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
 660   return TMath::ASinH(pt / mt * TMath::SinH(eta));
 661 }
 662
 663 TH1 *
 664 Convert_dNdy_1over2pipt_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
 665 {
 666
 667   TH1 *hout = hin->Clone("hout");
 668   hout->Reset();
 669   Double_t pt, mt, conv, val, vale;
 670   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
 671     pt = hin->GetBinCenter(ibin + 1);
 672     conv = eta2y(pt, mass, eta) / eta;
 673     val = hin->GetBinContent(ibin + 1);
 674     vale = hin->GetBinError(ibin + 1);
 675     val /= (2. * TMath::Pi() * pt);
 676     vale /= (2. * TMath::Pi() * pt);
 677     val *= conv;
 678     vale *= conv;
 679     hout->SetBinContent(ibin + 1, val);
 680     hout->SetBinError(ibin + 1, vale);
 681   }
 682   return hout;
 683 }
 684
 685 TH1 *
 686 Convert_dNdy_1over2pipt_dNdy(TH1 *hin)
 687 {
 688
 689   TH1 *hout = hin->Clone("hout");
 690   hout->Reset();
 691   Double_t pt, mt, conv, val, vale;
 692   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
 693     pt = hin->GetBinCenter(ibin + 1);
 694     val = hin->GetBinContent(ibin + 1);
 695     vale = hin->GetBinError(ibin + 1);
 696     val /= (2. * TMath::Pi() * pt);
 697     vale /= (2. * TMath::Pi() * pt);
 698     hout->SetBinContent(ibin + 1, val);
 699     hout->SetBinError(ibin + 1, vale);
 700   }
 701   return hout;
 702 }
 703
 704 TH1 *
 705 Convert_dNdy_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
 706 {
 707
 708   TH1 *hout = hin->Clone("hout");
 709   hout->Reset();
 710   Double_t pt, mt, conv, val, vale;
 711   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
 712     pt = hin->GetBinCenter(ibin + 1);
 713     conv = eta2y(pt, mass, eta) / eta;
 714     val = hin->GetBinContent(ibin + 1);
 715     vale = hin->GetBinError(ibin + 1);
 716     val *= conv;
 717     vale *= conv;
 718     hout->SetBinContent(ibin + 1, val);
 719     hout->SetBinError(ibin + 1, vale);
 720   }
 721   return hout;
 722 }
 723
 724 TGraph *
 725 Convert_dNdy_dNdeta(TGraph *hin, Double_t mass, Double_t eta = 0.8)
 726 {
 727
 728   TGraph *hout = hin->Clone("hout");
 729   //  hout->Reset();
 730   Double_t pt, mt, conv, val, valelo, valehi;
 731   for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
 732     pt = hin->GetX()[ibin];
 733     conv = eta2y(pt, mass, eta) / eta;
 734     val = hin->GetY()[ibin];
 735     valelo = hin->GetEYlow()[ibin];
 736     valehi = hin->GetEYhigh()[ibin];
 737     val *= conv;
 738     valelo *= conv;
 739     valehi *= conv;
 740     hout->GetX()[ibin] = pt;
 741     hout->GetY()[ibin] = val;
 742     hout->GetEYlow()[ibin] = valelo;
 743     hout->GetEYhigh()[ibin] = valehi;
 744   }
 745   return hout;
 746 }
 747
 748 TH1 *
 749 SummedId_1over2pipt_dNdeta(const Char_t *filename, Int_t icent)
 750 {
 751
 752   const Char_t *chargeName[2] = {
 753     "plus", "minus"
 754   };
 755
 756   TFile *filein = TFile::Open(filename);
 757   TH1 *hy[AliPID::kSPECIES][2];
 758   TH1 *heta[AliPID::kSPECIES][2];
 759   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
 760     for (Int_t icharge = 0; icharge < 2; icharge++) {
 761       hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
 762       if (!hy[ipart][icharge]) {
 763         printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
 764         return NULL;
 765       }
 766       heta[ipart][icharge] = Convert_dNdy_1over2pipt_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
 767     }
 768
 769   /* sum */
 770   TH1D *hsum = heta[2][0]->Clone("hsum");
 771   hsum->Reset();
 772   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
 773     for (Int_t icharge = 0; icharge < 2; icharge++)
 774       hsum->Add(heta[ipart][icharge]);
 775
 776   return hsum;
 777 }
 778
 779 TH1 *
 780 SummedId_dNdeta(const Char_t *filename, Int_t icent)
 781 {
 782
 783   const Char_t *chargeName[2] = {
 784     "plus", "minus"
 785   };
 786
 787   TFile *filein = TFile::Open(filename);
 788   TH1 *hy[AliPID::kSPECIES][2];
 789   TH1 *heta[AliPID::kSPECIES][2];
 790   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
 791     for (Int_t icharge = 0; icharge < 2; icharge++) {
 792       hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
 793       if (!hy[ipart][icharge]) {
 794         printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
 795         return NULL;
 796       }
 797       heta[ipart][icharge] = Convert_dNdy_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
 798     }
 799
 800   /* sum */
 801   TH1D *hsum = heta[2][0]->Clone("hsum");
 802   hsum->Reset();
 803   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
 804     for (Int_t icharge = 0; icharge < 2; icharge++)
 805       hsum->Add(heta[ipart][icharge]);
 806
 807   return hsum;
 808 }
 809
 810 /*****************************************************************/
 811