PWGLF/SPECTRA/UTILS/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   fLevyTsallis->SetParLimits(1, 1.e-3, 1.e3);
  73   fLevyTsallis->SetParLimits(2, 1.e-3, 1.e3);
  74   fLevyTsallis->SetParLimits(3, 1.e-6, 1.e6);
  75   return fLevyTsallis;
  76 }
  77
  78 /*****************************************************************/
  79 /* BOLTZMANN-GIBBS BLAST-WAVE */
  80 /*****************************************************************/
  81
  82 static TF1 *fBGBlastWave_Integrand = NULL;
  83 static TF1 *fBGBlastWave_Integrand_num = NULL;
  84 static TF1 *fBGBlastWave_Integrand_den = NULL;
  85 Double_t
  86 BGBlastWave_Integrand(const Double_t *x, const Double_t *p)
  87 {
  88
  89   /*
  90      x[0] -> r (radius)
  91      p[0] -> mT (transverse mass)
  92      p[1] -> pT (transverse momentum)
  93      p[2] -> beta_max (surface velocity)
  94      p[3] -> T (freezout temperature)
  95      p[4] -> n (velocity profile)
  96   */
  97
  98   Double_t r = x[0];
  99   Double_t mt = p[0];
 100   Double_t pt = p[1];
 101   Double_t beta_max = p[2];
 102   Double_t temp_1 = 1. / p[3];
 103   Double_t n = p[4];
 104
 105   Double_t beta = beta_max * TMath::Power(r, n);
 106   if (beta > 0.9999999999999999) beta = 0.9999999999999999;
 107   Double_t rho = TMath::ATanH(beta);
 108   Double_t argI0 = pt * TMath::SinH(rho) * temp_1;
 109   if (argI0 > 700.) argI0 = 700.;
 110   Double_t argK1 = mt * TMath::CosH(rho) * temp_1;
 111   //  if (argI0 > 100 || argI0 < -100)
 112   //    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);
 113   return r * mt * TMath::BesselI0(argI0) * TMath::BesselK1(argK1);
 114
 115 }
 116
 117 Double_t
 118 BGBlastWave_Func(const Double_t *x, const Double_t *p)
 119 {
 120   /* dN/dpt */
 121
 122   Double_t pt = x[0];
 123   Double_t mass = p[0];
 124   Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
 125   Double_t beta_max = p[1];
 126   Double_t temp = p[2];
 127   Double_t n = p[3];
 128   Double_t norm = p[4];
 129
 130   if (!fBGBlastWave_Integrand)
 131     fBGBlastWave_Integrand = new TF1("fBGBlastWave_Integrand", BGBlastWave_Integrand, 0., 1., 5);
 132   fBGBlastWave_Integrand->SetParameters(mt, pt, beta_max, temp, n);
 133   Double_t integral = fBGBlastWave_Integrand->Integral(0., 1., (Double_t *)0, 1.e-6);
 134   return norm * pt * integral;
 135 }
 136
 137 Double_t
 138 BGBlastWaveRatio_Func(const Double_t *x, const Double_t *p)
 139 {
 140   /* dN/dpt */
 141
 142   Double_t pt = x[0];
 143   Double_t mass = p[0];
 144   Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
 145   Double_t beta_max_num = p[1];
 146   Double_t temp_num = p[2];
 147   Double_t n_num = p[3];
 148   Double_t norm_num = p[4];
 149   Double_t beta_max_den = p[5];
 150   Double_t temp_den = p[6];
 151   Double_t n_den = p[7];
 152   Double_t norm_den = p[8];
 153
 154   if (!fBGBlastWave_Integrand_num)
 155     fBGBlastWave_Integrand_num = new TF1("fBGBlastWave_Integrand_num", BGBlastWave_Integrand, 0., 1., 5);
 156   fBGBlastWave_Integrand_num->SetParameters(mt, pt, beta_max_num, temp_num, n_num);
 157   Double_t integral_num = fBGBlastWave_Integrand_num->Integral(0., 1.);
 158
 159   if (!fBGBlastWave_Integrand_den)
 160     fBGBlastWave_Integrand_den = new TF1("fBGBlastWave_Integrand_den", BGBlastWave_Integrand, 0., 1., 5);
 161   fBGBlastWave_Integrand_den->SetParameters(mt, pt, beta_max_den, temp_den, n_den);
 162   Double_t integral_den = fBGBlastWave_Integrand_den->Integral(0., 1.);
 163
 164   return (norm_num / norm_den) * (integral_num / integral_den);
 165 }
 166
 167 Double_t
 168 BGBlastWaveParticleRatio_Func(const Double_t *x, const Double_t *p)
 169 {
 170   /* dN/dpt */
 171
 172   Double_t pt = x[0];
 173   Double_t mass_num = p[0];
 174   Double_t mass_den = p[1];
 175   Double_t mt_num = TMath::Sqrt(pt * pt + mass_num * mass_num);
 176   Double_t mt_den = TMath::Sqrt(pt * pt + mass_den * mass_den);
 177   Double_t beta_max = p[2];
 178   Double_t temp = p[3];
 179   Double_t n = p[4];
 180   Double_t norm_num = p[5];
 181   Double_t norm_den = p[6];
 182
 183   if (!fBGBlastWave_Integrand_num)
 184     fBGBlastWave_Integrand_num = new TF1("fBGBlastWave_Integrand_num", BGBlastWave_Integrand, 0., 1., 5);
 185   fBGBlastWave_Integrand_num->SetParameters(mt_num, pt, beta_max, temp, n);
 186   Double_t integral_num = fBGBlastWave_Integrand_num->Integral(0., 1.);
 187
 188   if (!fBGBlastWave_Integrand_den)
 189     fBGBlastWave_Integrand_den = new TF1("fBGBlastWave_Integrand_den", BGBlastWave_Integrand, 0., 1., 5);
 190   fBGBlastWave_Integrand_den->SetParameters(mt_den, pt, beta_max, temp, n);
 191   Double_t integral_den = fBGBlastWave_Integrand_den->Integral(0., 1.);
 192
 193   return (norm_num / norm_den) * (integral_num / integral_den);
 194 }
 195
 196 Double_t
 197 BGBlastWave_Func_OneOverPt(const Double_t *x, const Double_t *p)
 198 {
 199   /* 1/pt dN/dpt */
 200
 201   Double_t pt = x[0];
 202   Double_t mass = p[0];
 203   Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
 204   Double_t beta_max = p[1];
 205   Double_t temp = p[2];
 206   Double_t n = p[3];
 207   Double_t norm = p[4];
 208
 209   if (!fBGBlastWave_Integrand)
 210     fBGBlastWave_Integrand = new TF1("fBGBlastWave_Integrand", BGBlastWave_Integrand, 0., 1., 5);
 211   fBGBlastWave_Integrand->SetParameters(mt, pt, beta_max, temp, n);
 212   Double_t integral = fBGBlastWave_Integrand->Integral(0., 1., (Double_t *)0, 1.e-3);
 213
 214   return norm * integral;
 215 }
 216
 217 TF1 *
 218 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)
 219 {
 220
 221   TF1 *fBGBlastWave = new TF1(name, BGBlastWave_Func, 0., 10., 5);
 222   fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm);
 223   fBGBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
 224   fBGBlastWave->FixParameter(0, mass);
 225   fBGBlastWave->SetParLimits(1, 0.01, 0.99);
 226   fBGBlastWave->SetParLimits(2, 0.01, 1.);
 227   fBGBlastWave->SetParLimits(3, 0.01, 50.);
 228   return fBGBlastWave;
 229 }
 230
 231 TF1 *
 232 BGBlastWaveRatio(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)
 233 {
 234
 235   TF1 *fBGBlastWave = new TF1(name, BGBlastWaveRatio_Func, 0., 10., 9);
 236   fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm, beta_max, temp, n, norm);
 237   fBGBlastWave->SetParNames("mass", "beta_max_num", "T_num", "n_num", "norm_num", "beta_max_den", "T_den", "n_den", "norm_den");
 238   fBGBlastWave->FixParameter(0, mass);
 239   fBGBlastWave->SetParLimits(1, 0.01, 0.99);
 240   fBGBlastWave->SetParLimits(2, 0.01, 1.);
 241   fBGBlastWave->SetParLimits(3, 0.01, 10.);
 242   fBGBlastWave->SetParLimits(5, 0.01, 0.99);
 243   fBGBlastWave->SetParLimits(6, 0.01, 1.);
 244   fBGBlastWave->SetParLimits(7, 0.01, 10.);
 245   return fBGBlastWave;
 246 }
 247
 248 TF1 *
 249 BGBlastWaveParticleRatio(const Char_t *name, Double_t mass_num, Double_t mass_den, Double_t beta_max = 0.9, Double_t temp = 0.1, Double_t n = 1., Double_t norm_num = 1.e6, Double_t norm_den = 1.e6)
 250 {
 251
 252   TF1 *fBGBlastWave = new TF1(name, BGBlastWaveParticleRatio_Func, 0., 10., 7);
 253   fBGBlastWave->SetParameters(mass_num, mass_den, beta_max, temp, n, norm_num, norm_den);
 254   fBGBlastWave->SetParNames("mass_num", "mass_den", "beta_max", "T", "n", "norm_num", "norm_den");
 255   fBGBlastWave->FixParameter(0, mass_num);
 256   fBGBlastWave->FixParameter(1, mass_den);
 257   fBGBlastWave->SetParLimits(2, 0.01, 0.99);
 258   fBGBlastWave->SetParLimits(3, 0.01, 1.);
 259   fBGBlastWave->SetParLimits(4, 0.01, 10.);
 260   return fBGBlastWave;
 261 }
 262
 263 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)
 264 {
 265
 266   TF1 *fBGBlastWave = new TF1(name, BGBlastWave_Func_OneOverPt, 0., 10., 5);
 267   fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm);
 268   fBGBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
 269   fBGBlastWave->FixParameter(0, mass);
 270   fBGBlastWave->SetParLimits(1, 0.01, 0.99);
 271   fBGBlastWave->SetParLimits(2, 0.01, 1.);
 272   fBGBlastWave->SetParLimits(3, 0.01, 50.);
 273   return fBGBlastWave;
 274 }
 275
 276 /*****************************************************************/
 277 /* TSALLIS BLAST-WAVE */
 278 /*****************************************************************/
 279
 280 static TF1 *fTsallisBlastWave_Integrand_r = NULL;
 281 Double_t
 282 TsallisBlastWave_Integrand_r(const Double_t *x, const Double_t *p)
 283 {
 284   /*
 285      x[0] -> r (radius)
 286      p[0] -> mT (transverse mass)
 287      p[1] -> pT (transverse momentum)
 288      p[2] -> beta_max (surface velocity)
 289      p[3] -> T (freezout temperature)
 290      p[4] -> n (velocity profile)
 291      p[5] -> q
 292      p[6] -> y (rapidity)
 293      p[7] -> phi (azimuthal angle)
 294   */
 295
 296   Double_t r = x[0];
 297   Double_t mt = p[0];
 298   Double_t pt = p[1];
 299   Double_t beta_max = p[2];
 300   Double_t temp_1 = 1. / p[3];
 301   Double_t n = p[4];
 302   Double_t q = p[5];
 303   Double_t y = p[6];
 304   Double_t phi = p[7];
 305
 306   if (q <= 1.) return r;
 307
 308   Double_t beta = beta_max * TMath::Power(r, n);
 309   Double_t rho = TMath::ATanH(beta);
 310
 311   Double_t part1 = mt * TMath::CosH(y) * TMath::CosH(rho);
 312   Double_t part2 = pt * TMath::SinH(rho) * TMath::Cos(phi);
 313   Double_t part3 = part1 - part2;
 314   Double_t part4 = 1 + (q - 1.) * temp_1 * part3;
 315   Double_t expo = -1. / (q - 1.);
 316   //  printf("part1=%f, part2=%f, part3=%f, part4=%f, expo=%f\n", part1, part2, part3, part4, expo);
 317   Double_t part5 = TMath::Power(part4, expo);
 318
 319   return r * part5;
 320 }
 321
 322 static TF1 *fTsallisBlastWave_Integrand_phi = NULL;
 323 Double_t
 324 TsallisBlastWave_Integrand_phi(const Double_t *x, const Double_t *p)
 325 {
 326   /*
 327      x[0] -> phi (azimuthal angle)
 328   */
 329
 330   Double_t phi = x[0];
 331   fTsallisBlastWave_Integrand_r->SetParameter(7, phi);
 332   Double_t integral = fTsallisBlastWave_Integrand_r->Integral(0., 1.);
 333   return integral;
 334 }
 335
 336 static TF1 *fTsallisBlastWave_Integrand_y = NULL;
 337 Double_t
 338 TsallisBlastWave_Integrand_y(const Double_t *x, const Double_t *p)
 339 {
 340   /*
 341      x[0] -> y (rapidity)
 342   */
 343
 344   Double_t y = x[0];
 345   fTsallisBlastWave_Integrand_r->SetParameter(6, y);
 346   Double_t integral = fTsallisBlastWave_Integrand_phi->Integral(-TMath::Pi(), TMath::Pi());
 347   return TMath::CosH(y) * integral;
 348 }
 349
 350 Double_t
 351 TsallisBlastWave_Func(const Double_t *x, const Double_t *p)
 352 {
 353   /* dN/dpt */
 354
 355   Double_t pt = x[0];
 356   Double_t mass = p[0];
 357   Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
 358   Double_t beta_max = p[1];
 359   Double_t temp = p[2];
 360   Double_t n = p[3];
 361   Double_t q = p[4];
 362   Double_t norm = p[5];
 363
 364   if (!fTsallisBlastWave_Integrand_r)
 365     fTsallisBlastWave_Integrand_r = new TF1("fTsallisBlastWave_Integrand_r", TsallisBlastWave_Integrand_r, 0., 1., 8);
 366   if (!fTsallisBlastWave_Integrand_phi)
 367     fTsallisBlastWave_Integrand_phi = new TF1("fTsallisBlastWave_Integrand_phi", TsallisBlastWave_Integrand_phi, -TMath::Pi(), TMath::Pi(), 0);
 368   if (!fTsallisBlastWave_Integrand_y)
 369     fTsallisBlastWave_Integrand_y = new TF1("fTsallisBlastWave_Integrand_y", TsallisBlastWave_Integrand_y, -0.5, 0.5, 0);
 370
 371   fTsallisBlastWave_Integrand_r->SetParameters(mt, pt, beta_max, temp, n, q, 0., 0.);
 372   Double_t integral = fTsallisBlastWave_Integrand_y->Integral(-0.5, 0.5);
 373   return norm * pt * integral;
 374 }
 375
 376 TF1 *
 377 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)
 378 {
 379
 380   TF1 *fTsallisBlastWave = new TF1(name, TsallisBlastWave_Func, 0., 10., 6);
 381   fTsallisBlastWave->SetParameters(mass, beta_max, temp, n, q, norm);
 382   fTsallisBlastWave->SetParNames("mass", "beta_max", "T", "n", "q", "norm");
 383   fTsallisBlastWave->FixParameter(0, mass);
 384   fTsallisBlastWave->SetParLimits(1, 0.01, 0.99);
 385   fTsallisBlastWave->SetParLimits(2, 0.01, 1.);
 386   fTsallisBlastWave->SetParLimits(3, 0.1, 10.);
 387   fTsallisBlastWave->SetParLimits(4, 1., 10.);
 388   return fTsallisBlastWave;
 389 }
 390
 391 /*****************************************************************/
 392 /*****************************************************************/
 393 /*****************************************************************/
 394
 395
 396 TF1 *
 397 BGBlastWave_SingleFit(TH1 *h, Double_t mass, Option_t *opt = "")
 398 {
 399
 400   TF1 *f = BGBlastWave(Form("fBGBW_%s", h->GetName()), mass);
 401   h->Fit(f);
 402   h->Fit(f);
 403   h->Fit(f, opt);
 404   return f;
 405
 406 }
 407
 408 Int_t nBW;
 409 TF1 *fBGBW[1000];
 410 TF1 *fBGBWratio[1000];
 411 TGraphErrors *gBW[1000];
 412
 413 TObjArray *
 414 BGBlastWave_GlobalFit(TObjArray *data, Double_t *mass, Double_t profile = 0.5, Bool_t computeCont = kFALSE, Bool_t fixProfile = kFALSE)
 415 {
 416
 417   /* get data */
 418   Int_t ndf = 0;
 419   nBW = data->GetEntries();
 420   for (Int_t idata = 0; idata < nBW; idata++) {
 421     gBW[idata] = (TGraphErrors *)data->At(idata);
 422     gBW[idata]->SetName(Form("gBW%d", idata));
 423     ndf += gBW[idata]->GetN();
 424   }
 425
 426   /* init BG blast-wave functions */
 427   for (Int_t idata = 0; idata < nBW; idata++) {
 428     printf("init BG-BlastWave function #%d: mass = %f\n", idata, mass[idata]);
 429     fBGBW[idata] = BGBlastWave(Form("fBGBW%d", idata), mass[idata]);
 430   }
 431
 432   if (computeCont)
 433     printf("-> compute contours requested\n");
 434
 435   /* display data */
 436   TCanvas *cBW = new TCanvas("cBW");
 437   cBW->Divide(nBW, 1);
 438   for (Int_t idata = 0; idata < nBW; idata++) {
 439     cBW->cd(idata + 1);
 440     gBW[idata]->Draw("ap*");
 441   }
 442   cBW->Update();
 443
 444   /* init minuit: nBW normalizations + 3 (beta, T, n) BG-BlastWave params */
 445   const Int_t nbwpars = 3;
 446   const Int_t nfitpars = nBW + nbwpars;
 447   TMinuit *minuit = new TMinuit(nfitpars);
 448   minuit->SetFCN(BGBlastWave_FCN);
 449   Double_t arglist[10];
 450   Int_t ierflg = 0;
 451   arglist[0] = 1;
 452   minuit->mnexcm("SET ERR", arglist, 1, ierflg);
 453   for (Int_t idata = 0; idata < nBW; idata++) {
 454     minuit->mnparm(idata, Form("norm%d", idata), 1.e6, 1., 0., 0., ierflg);
 455     ndf--;
 456   }
 457   //  minuit->mnparm(nBW + 0, "<beta>", 0.55, 0.01, 0., 1., ierflg);
 458   //  minuit->mnparm(nBW + 1, "T", 0.14, 0.01, 0., 1., ierflg);
 459   //  minuit->mnparm(nBW + 2, "n", profile, 0.1, 0., 10., ierflg);
 460
 461   minuit->mnparm(nBW + 0, "<beta>", 0.7, 0.01, 0.2, 0.7, ierflg);
 462   minuit->mnparm(nBW + 1, "T", 0.07, 0.001, 0.07, 0.2, ierflg);
 463   minuit->mnparm(nBW + 2, "n", profile, 0.1, 0.6, 5., ierflg);
 464
 465   ndf -= 3;
 466   if (fixProfile) {
 467     minuit->FixParameter(nBW + 2);
 468     ndf++;
 469   }
 470
 471   /* set strategy */
 472   arglist[0] = 1;
 473   minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
 474
 475   /* start MIGRAD minimization */
 476   arglist[0] = 500000;
 477   arglist[1] = 1.;
 478   minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
 479
 480   /* set strategy */
 481   arglist[0] = 2;
 482   minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
 483
 484   /* start MIGRAD minimization */
 485   arglist[0] = 500000;
 486   arglist[1] = 1.;
 487   minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
 488
 489   /* start IMPROVE minimization */
 490   arglist[0] = 500000;
 491   minuit->mnexcm("IMPROVE", arglist, 1, ierflg);
 492
 493   /* start MINOS */
 494   arglist[0] = 500000;
 495   arglist[1] = nBW + 1;
 496   arglist[2] = nBW + 2;
 497   arglist[3] = nBW + 3;
 498   minuit->mnexcm("MINOS", arglist, 4, ierflg);
 499
 500   /* print results */
 501   Double_t amin,edm,errdef;
 502   Int_t nvpar,nparx,icstat;
 503   minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
 504   minuit->mnprin(4, amin);
 505
 506   /* get parameters */
 507   Double_t beta, betae, betaeplus, betaeminus, betagcc, temp, tempe, tempeplus, tempeminus, tempgcc, prof, profe, profeplus, profeminus, profgcc;
 508   minuit->GetParameter(nBW + 0, beta, betae);
 509   minuit->mnerrs(nBW + 0, betaeplus, betaeminus, betae, betagcc);
 510   minuit->GetParameter(nBW + 1, temp, tempe);
 511   minuit->mnerrs(nBW + 1, tempeplus, tempeminus, tempe, tempgcc);
 512   minuit->GetParameter(nBW + 2, prof, profe);
 513   minuit->mnerrs(nBW + 2, profeplus, profeminus, profe, profgcc);
 514   Double_t beta_max = 0.5 * (2. + prof) * beta;
 515   Double_t norm[1000], norme[1000];
 516   for (Int_t idata = 0; idata < nBW; idata++)
 517     minuit->GetParameter(idata, norm[idata], norme[idata]);
 518
 519   /* printout */
 520   printf("[x] *********************************\n");
 521   printf("[x] beta_max = %f\n", beta_max);
 522   printf("[x] <beta>   = %f +- %f (e+ = %f, e- = %f)\n", beta, betae, betaeplus, betaeminus);
 523   printf("[x] T        = %f +- %f (e+ = %f, e- = %f)\n", temp, tempe, tempeplus, tempeminus);
 524   printf("[x] n        = %f +- %f (e+ = %f, e- = %f)\n", prof, profe, profeplus, profeminus);
 525   printf("[x] chi2     = %f\n", amin);
 526   printf("[x] ndf      = %f\n", ndf);
 527
 528   /* 1-sigma contour */
 529   minuit->SetErrorDef(1);
 530   TGraph *gCont1 = NULL;
 531   if (computeCont) gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
 532   if (gCont1) gCont1->SetName("gCont1");
 533
 534   /* 2-sigma contour */
 535   minuit->SetErrorDef(4);
 536   TGraph *gCont2 = NULL;
 537   //  if (computeCont) gCont2 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
 538   if (gCont2) gCont2->SetName("gCont2");
 539
 540   /* display fits */
 541   for (Int_t idata = 0; idata < nBW; idata++) {
 542     cBW->cd(idata + 1);
 543     fBGBW[idata]->SetParameter(4, norm[idata]);
 544     fBGBW[idata]->SetParameter(1, beta_max);
 545     fBGBW[idata]->SetParameter(2, temp);
 546     fBGBW[idata]->SetParameter(3, prof);
 547     fBGBW[idata]->Draw("same");
 548   }
 549   cBW->Update();
 550
 551   /* histo params */
 552   TH1D *hBW = new TH1D("hBW", "", 4, 0., 4.);
 553   hBW->SetBinContent(1, beta);
 554   hBW->SetBinError(1, betae);
 555   hBW->SetBinContent(2, temp);
 556   hBW->SetBinError(2, tempe);
 557   hBW->SetBinContent(3, prof);
 558   hBW->SetBinError(3, profe);
 559   hBW->SetBinContent(4, amin/ndf);
 560
 561   /* BW graph */
 562   TGraphAsymmErrors *gBetaT = new TGraphAsymmErrors();
 563   gBetaT->SetName("gBetaT");
 564   gBetaT->SetPoint(0, beta, temp);
 565   gBetaT->SetPointEXlow(0, TMath::Abs(betaeminus));
 566   gBetaT->SetPointEXhigh(0, TMath::Abs(betaeplus));
 567   gBetaT->SetPointEYlow(0, TMath::Abs(tempeminus));
 568   gBetaT->SetPointEYhigh(0, TMath::Abs(tempeplus));
 569
 570   /* prepare output array */
 571   TObjArray *outoa = new TObjArray();
 572   for (Int_t idata = 0; idata < nBW; idata++) {
 573     outoa->Add(gBW[idata]);
 574     outoa->Add(fBGBW[idata]);
 575   }
 576   outoa->Add(cBW);
 577   outoa->Add(hBW);
 578   outoa->Add(gBetaT);
 579   if (gCont1) outoa->Add(gCont1);
 580   if (gCont2) outoa->Add(gCont2);
 581
 582   return outoa;
 583
 584 }
 585
 586 void
 587 BGBlastWave_FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
 588 {
 589
 590   /* beta -> beta_max */
 591   Double_t beta = par[nBW+0];
 592   Double_t T = par[nBW+1];
 593   Double_t n = par[nBW+2];
 594   Double_t beta_max = 0.5 * (2. + n) * beta;
 595 #if 1
 596   /* check beta_max */
 597   if (beta_max >= 1. || beta_max <= 0.) {
 598     f = kMaxInt;
 599     return;
 600   }
 601   /* check T */
 602   if (T <= 0.) {
 603     f = kMaxInt;
 604     return;
 605   }
 606 #endif
 607
 608   Double_t pt, pte, val, vale, func, pull, chi = 0;
 609   /* loop over all the data */
 610   for (Int_t iBW = 0; iBW < nBW; iBW++) {
 611     /* set BGBW parameters */
 612     fBGBW[iBW]->SetParameter(4, par[iBW]);
 613     fBGBW[iBW]->SetParameter(1, beta_max);
 614     fBGBW[iBW]->SetParameter(2, T);
 615     fBGBW[iBW]->SetParameter(3, n);
 616     /* loop over all the points */
 617     for (Int_t ipt = 0; ipt < gBW[iBW]->GetN(); ipt++) {
 618       pt = gBW[iBW]->GetX()[ipt];
 619       pte = gBW[iBW]->GetEX()[ipt];
 620       val = gBW[iBW]->GetY()[ipt];
 621       vale = gBW[iBW]->GetEY()[ipt];
 622       func = fBGBW[iBW]->Eval(pt);
 623       //      func = fBGBW[iBW]->Integral(pt - pte, pt + pte);
 624       pull = (val - func) / vale;
 625       chi += pull * pull;
 626     }
 627   }
 628
 629   f = chi;
 630 }
 631
 632 /*****************************************************************/
 633
 634 TObjArray *
 635 BGBlastWave_GlobalFitRatio(TObjArray *data, Double_t *mass, Double_t profile = .7, Bool_t fixProfile = kFALSE)
 636 {
 637
 638   /* get data */
 639   Int_t ndf = 0;
 640   nBW = data->GetEntries();
 641   for (Int_t idata = 0; idata < nBW; idata++) {
 642     gBW[idata] = (TGraphErrors *)data->At(idata);
 643     gBW[idata]->SetName(Form("gBW%d", idata));
 644     ndf += gBW[idata]->GetN();
 645   }
 646
 647   /* init BG blast-wave functions */
 648   for (Int_t idata = 0; idata < nBW; idata++) {
 649     printf("init BG-BlastWaveRatio function #%d: mass = %f\n", idata, mass[idata]);
 650     fBGBWratio[idata] = BGBlastWaveRatio(Form("fBGBWratio%d", idata), mass[idata]);
 651   }
 652
 653   /* display data */
 654   TCanvas *cBW = new TCanvas("cBW");
 655   cBW->Divide(nBW, 1);
 656   for (Int_t idata = 0; idata < nBW; idata++) {
 657     cBW->cd(idata + 1);
 658     gBW[idata]->Draw("ap*");
 659   }
 660   cBW->Update();
 661
 662   /* init minuit: nBW normalizations + 3 (beta, T, n) BG-BlastWave params */
 663   const Int_t nbwpars = 3;
 664   const Int_t nfitpars = 2 * (nBW + nbwpars);
 665   TMinuit *minuit = new TMinuit(nfitpars);
 666   minuit->SetFCN(BGBlastWave_FCNRatio);
 667   Double_t arglist[10];
 668   Int_t ierflg = 0;
 669   arglist[0] = 1;
 670   minuit->mnexcm("SET ERR", arglist, 1, ierflg);
 671   for (Int_t idata = 0; idata < nBW; idata++) {
 672     minuit->mnparm(idata, Form("norm%d_num", idata), 1.e6, 1., 0., 0., ierflg);
 673     ndf--;
 674   }
 675   for (Int_t idata = nBW; idata < 2 * nBW; idata++) {
 676     minuit->mnparm(idata, Form("norm%d_den", idata), 1.e6, 1., 0., 0., ierflg);
 677     minuit->FixParameter(idata);
 678     ndf--;
 679   }
 680   minuit->mnparm(2 * nBW + 0, "<beta>_num", 0.65, 0.01, 0., 1., ierflg);
 681   minuit->mnparm(2 * nBW + 1, "T_num", 0.1, 0.01, 0., 1., ierflg);
 682   minuit->mnparm(2 * nBW + 2, "n_num", profile, 0.1, 0., 10., ierflg);
 683   minuit->mnparm(2 * nBW + 3, "<beta>_den", 0.65, 0.01, 0., 1., ierflg);
 684   minuit->mnparm(2 * nBW + 4, "T_den", 0.1, 0.01, 0., 1., ierflg);
 685   minuit->mnparm(2 * nBW + 5, "n_den", profile, 0.1, 0., 10., ierflg);
 686   ndf -= 3;
 687
 688   if (fixProfile) {
 689     minuit->FixParameter(nBW + 2);
 690     minuit->FixParameter(nBW + 5);
 691     ndf++;
 692   }
 693
 694   /* set strategy */
 695   arglist[0] = 1;
 696   minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
 697
 698   printf("-->> STARTING MIGRAD with %d parameters\n", nfitpars);
 699
 700   /* start MIGRAD minimization */
 701   arglist[0] = 500000;
 702   arglist[1] = 1.;
 703   minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
 704
 705   /* set strategy */
 706   arglist[0] = 2;
 707   minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
 708
 709   /* start MIGRAD minimization */
 710   arglist[0] = 500000;
 711   arglist[1] = 1.;
 712   minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
 713
 714   /* start IMPROVE minimization */
 715   arglist[0] = 500000;
 716   minuit->mnexcm("IMPROVE", arglist, 1, ierflg);
 717
 718   /* start MINOS */
 719   arglist[0] = 500000;
 720   arglist[1] = 2 * nBW + 1;
 721   arglist[2] = 2 * nBW + 2;
 722   arglist[3] = 2 * nBW + 3;
 723   arglist[4] = 2 * nBW + 4;
 724   arglist[5] = 2 * nBW + 5;
 725   arglist[6] = 2 * nBW + 6;
 726   minuit->mnexcm("MINOS", arglist, 7, ierflg);
 727
 728   /* print results */
 729   Double_t amin,edm,errdef;
 730   Int_t nvpar,nparx,icstat;
 731   minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
 732   minuit->mnprin(4, amin);
 733
 734   /* get parameters */
 735   Double_t beta_num, betae_num, betaeplus_num, betaeminus_num, betagcc_num, temp_num, tempe_num, tempeplus_num, tempeminus_num, tempgcc_num, prof_num, profe_num, profeplus_num, profeminus_num, profgcc_num;
 736   Double_t beta_den, betae_den, betaeplus_den, betaeminus_den, betagcc_den, temp_den, tempe_den, tempeplus_den, tempeminus_den, tempgcc_den, prof_den, profe_den, profeplus_den, profeminus_den, profgcc_den;
 737   minuit->GetParameter(2*nBW + 0, beta_num, betae_num);
 738   minuit->mnerrs(nBW + 0, betaeplus_num, betaeminus_num, betae_num, betagcc_num);
 739   minuit->GetParameter(2*nBW + 1, temp_num, tempe_num);
 740   minuit->mnerrs(nBW + 1, tempeplus_num, tempeminus_num, tempe_num, tempgcc_num);
 741   minuit->GetParameter(2*nBW + 2, prof_num, profe_num);
 742   minuit->mnerrs(nBW + 2, profeplus_num, profeminus_num, profe_num, profgcc_num);
 743   minuit->GetParameter(2*nBW + 3, beta_den, betae_den);
 744   minuit->mnerrs(nBW + 3, betaeplus_den, betaeminus_den, betae_den, betagcc_den);
 745   minuit->GetParameter(2*nBW + 4, temp_den, tempe_den);
 746   minuit->mnerrs(nBW + 4, tempeplus_den, tempeminus_den, tempe_den, tempgcc_den);
 747   minuit->GetParameter(2*nBW + 5, prof_den, profe_den);
 748   minuit->mnerrs(nBW + 5, profeplus_den, profeminus_den, profe_den, profgcc_den);
 749   Double_t beta_max_num, beta_max_den;
 750   beta_max_num = 0.5 * (2. + prof_num) * beta_num;
 751   beta_max_den = 0.5 * (2. + prof_den) * beta_den;
 752   Double_t norm_num[1000], norme_num[1000];
 753   Double_t norm_den[1000], norme_den[1000];
 754   for (Int_t idata = 0; idata < nBW; idata++)
 755     minuit->GetParameter(idata, norm_num[idata], norme_num[idata]);
 756   for (Int_t idata = 0; idata < nBW; idata++) {
 757     minuit->GetParameter(nBW + idata, norm_den[idata], norme_den[idata]);
 758   }
 759
 760   /* printout */
 761   printf("[x] *********************************\n");
 762   printf("[x] beta_max = %f\n", beta_max_num);
 763   printf("[x] <beta>   = %f +- %f (e+ = %f, e- = %f)\n", beta_num, betae_num, betaeplus_num, betaeminus_num);
 764   printf("[x] T        = %f +- %f (e+ = %f, e- = %f)\n", temp_num, tempe_num, tempeplus_num, tempeminus_num);
 765   printf("[x] n        = %f +- %f (e+ = %f, e- = %f)\n", prof_num, profe_num, profeplus_num, profeminus_num);
 766   printf("[x] *********************************\n");
 767   printf("[x] beta_max = %f\n", beta_max_den);
 768   printf("[x] <beta>   = %f +- %f (e+ = %f, e- = %f)\n", beta_den, betae_den, betaeplus_den, betaeminus_den);
 769   printf("[x] T        = %f +- %f (e+ = %f, e- = %f)\n", temp_den, tempe_den, tempeplus_den, tempeminus_den);
 770   printf("[x] n        = %f +- %f (e+ = %f, e- = %f)\n", prof_den, profe_den, profeplus_den, profeminus_den);
 771   printf("[x] *********************************\n");
 772   printf("[x] chi2     = %f\n", amin);
 773   printf("[x] ndf      = %f\n", ndf);
 774
 775   /* 1-sigma contour */
 776   minuit->SetErrorDef(1);
 777   TGraph *gCont1 = NULL;
 778   //  gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
 779   if (gCont1) gCont1->SetName("gCont1");
 780
 781   /* 2-sigma contour */
 782   minuit->SetErrorDef(4);
 783   TGraph *gCont2 = NULL;
 784   //  gCont2 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
 785   if (gCont2) gCont2->SetName("gCont2");
 786
 787   /* display fits */
 788   for (Int_t idata = 0; idata < nBW; idata++) {
 789     cBW->cd(idata + 1);
 790     fBGBWratio[idata]->SetParameter(4, norm_num[idata]);
 791     fBGBWratio[idata]->SetParameter(1, beta_max_num);
 792     fBGBWratio[idata]->SetParameter(2, temp_num);
 793     fBGBWratio[idata]->SetParameter(3, prof_num);
 794     fBGBWratio[idata]->SetParameter(8, norm_den[idata]);
 795     fBGBWratio[idata]->SetParameter(5, beta_max_den);
 796     fBGBWratio[idata]->SetParameter(6, temp_den);
 797     fBGBWratio[idata]->SetParameter(7, prof_den);
 798     fBGBWratio[idata]->Draw("same");
 799   }
 800   cBW->Update();
 801
 802   return;
 803
 804   /* histo params */
 805   TH1D *hBW = new TH1D("hBW", "", 3, 0., 3.);
 806   hBW->SetBinContent(1, beta);
 807   hBW->SetBinError(1, betae);
 808   hBW->SetBinContent(2, temp);
 809   hBW->SetBinError(2, tempe);
 810   hBW->SetBinContent(3, prof);
 811   hBW->SetBinError(3, profe);
 812
 813   /* BW graph */
 814   TGraphAsymmErrors *gBetaT = new TGraphAsymmErrors();
 815   gBetaT->SetName("gBetaT");
 816   gBetaT->SetPoint(0, beta, temp);
 817   gBetaT->SetPointEXlow(0, TMath::Abs(betaeminus));
 818   gBetaT->SetPointEXhigh(0, TMath::Abs(betaeplus));
 819   gBetaT->SetPointEYlow(0, TMath::Abs(tempeminus));
 820   gBetaT->SetPointEYhigh(0, TMath::Abs(tempeplus));
 821
 822   /* prepare output array */
 823   TObjArray *outoa = new TObjArray();
 824   for (Int_t idata = 0; idata < nBW; idata++) {
 825     outoa->Add(gBW[idata]);
 826     outoa->Add(fBGBWratio[idata]);
 827   }
 828   outoa->Add(cBW);
 829   outoa->Add(hBW);
 830   outoa->Add(gBetaT);
 831   if (gCont1) outoa->Add(gCont1);
 832   if (gCont2) outoa->Add(gCont2);
 833
 834   return outoa;
 835
 836 }
 837
 838 void
 839 BGBlastWave_FCNRatio(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
 840 {
 841
 842   /* beta -> beta_max */
 843   Double_t beta_num = par[2*nBW+0];
 844   Double_t T_num = par[2*nBW+1];
 845   Double_t n_num = par[2*nBW+2];
 846   Double_t beta_max_num = 0.5 * (2. + n_num) * beta_num;
 847   Double_t beta_den = par[2*nBW+3];
 848   Double_t T_den = par[2*nBW+4];
 849   Double_t n_den = par[2*nBW+5];
 850   Double_t beta_max_den = 0.5 * (2. + n_den) * beta_den;
 851 #if 0
 852   /* check beta_max */
 853   if (beta_max >= 1. || beta_max <= 0.) {
 854     f = kMaxInt;
 855     return;
 856   }
 857   /* check T */
 858   if (T <= 0.) {
 859     f = kMaxInt;
 860     return;
 861   }
 862 #endif
 863
 864   Double_t pt, pte, val, vale, func, pull, chi = 0;
 865   /* loop over all the data */
 866   for (Int_t iBW = 0; iBW < nBW; iBW++) {
 867     /* set BGBW parameters */
 868     fBGBWratio[iBW]->SetParameter(4, par[iBW]);
 869     fBGBWratio[iBW]->SetParameter(1, beta_max_num);
 870     fBGBWratio[iBW]->SetParameter(2, T_num);
 871     fBGBWratio[iBW]->SetParameter(3, n_num);
 872     fBGBWratio[iBW]->SetParameter(8, par[nBW+iBW]);
 873     fBGBWratio[iBW]->SetParameter(5, beta_max_den);
 874     fBGBWratio[iBW]->SetParameter(6, T_den);
 875     fBGBWratio[iBW]->SetParameter(7, n_den);
 876     /* loop over all the points */
 877     for (Int_t ipt = 0; ipt < gBW[iBW]->GetN(); ipt++) {
 878       pt = gBW[iBW]->GetX()[ipt];
 879       pte = gBW[iBW]->GetEX()[ipt];
 880       val = gBW[iBW]->GetY()[ipt];
 881       vale = gBW[iBW]->GetEY()[ipt];
 882       func = fBGBWratio[iBW]->Eval(pt);
 883       //      func = fBGBW[iBW]->Integral(pt - pte, pt + pte);
 884       pull = (val - func) / vale;
 885       chi += pull * pull;
 886     }
 887   }
 888
 889   f = chi;
 890 }
 891
 892 /*****************************************************************/
 893
 894 Int_t ratioPartNum[9] = {2, 3, 4, 3, 3, 3, 4, 4, 4};
 895 Int_t ratioChargeNum[9] = {1, 1, 1, 0, 1, 2, 0, 1, 2};
 896 Int_t ratioPartDen[9] = {2, 3, 4, 2, 2, 2, 2, 2, 2};
 897 Int_t ratioChargeDen[9] = {0, 0, 0, 0, 1, 2, 0, 1, 2};
 898
 899 Char_t *ratioPartNumName[9] = {"pion", "kaon", "proton", "kaon", "kaon", "kaon", "proton", "proton", "proton"};
 900 Char_t *ratioChargeNumName[9] = {"minus", "minus", "minus", "plus", "minus", "both", "plus", "minus", "both"};
 901 Char_t *ratioPartDenName[9] = {"pion", "kaon", "proton", "pion", "pion", "pion", "pion", "pion", "pion"};
 902 Char_t *ratioChargeDenName[9] = {"plus", "plus", "plus", "plus", "minus", "both", "plus", "minus", "both"};
 903
 904 IntegratedRatioError(const Char_t *spectrafilename, const Char_t *ratiosfilename)
 905 {
 906
 907   TFile *fspectra = TFile::Open(spectrafilename);
 908   TFile *fratios = TFile::Open(ratiosfilename);
 909
 910   Int_t icent = 0;
 911   for (Int_t irat = 3; irat < 9; irat++) {
 912
 913     Int_t ipnum = ratioPartNum[irat];
 914     Int_t icnum = ratioChargeNum[irat];
 915
 916     Int_t ipden = ratioPartDen[irat];
 917     Int_t icden = ratioChargeDen[irat];
 918
 919     if (icnum == 2 || icden == 2) continue;
 920
 921     printf("%s\n", Form("sys_cent%d_%s_%s_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat], ratioPartDenName[irat], ratioChargeDenName[irat]));
 922     TH1 *hrat_sys = (TH1 *)fratios->Get(Form("sys_cent%d_%s_%s_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat], ratioPartDenName[irat], ratioChargeDenName[irat]));
 923
 924     printf("%d %d %d %d\n", ipnum, icnum, ipden, icden);
 925     TH1 *hnum = (TH1 *)fspectra->Get(Form("cent%d_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat]));
 926     TH1 *hnum_stat = (TH1 *)fspectra->Get(Form("stat_cent%d_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat]));
 927     printf("%s = %p\n", Form("cent%d_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat]), hnum);
 928
 929     TH1 *hden = (TH1 *)fspectra->Get(Form("cent%d_%s_%s", icent, ratioPartDenName[irat], ratioChargeDenName[irat]));
 930     TH1 *hden_stat = (TH1 *)fspectra->Get(Form("stat_cent%d_%s_%s", icent, ratioPartDenName[irat], ratioChargeDenName[irat]));
 931     printf("%s = %p\n", Form("cent%d_%s_%s", icent, ratioPartDenName[irat], ratioChargeDenName[irat]), hden);
 932
 933     TH1 *hnum_plus = (TH1 *)hnum->Clone("hnum_plus");
 934     TH1 *hnum_minus = (TH1 *)hnum->Clone("hnum_minus");
 935     TH1 *hden_plus = (TH1 *)hden->Clone("hden_plus");
 936     TH1 *hden_minus = (TH1 *)hden->Clone("hden_minus");
 937     TH1 *hnum_stat_plus = (TH1 *)hnum_stat->Clone("hnum_stat_plus");
 938     TH1 *hnum_stat_minus = (TH1 *)hnum_stat->Clone("hnum_stat_minus");
 939     TH1 *hden_stat_plus = (TH1 *)hden_stat->Clone("hden_stat_plus");
 940     TH1 *hden_stat_minus = (TH1 *)hden_stat->Clone("hden_stat_minus");
 941     for (Int_t ibin = 0; ibin < hnum_plus->GetNbinsX(); ibin++) {
 942       Double_t val = hrat_sys->GetBinContent(ibin + 1);
 943       Double_t vale = hrat_sys->GetBinError(ibin + 1);
 944       if (vale <= 0.) continue;
 945       Double_t syse = vale / val;
 946
 947       val = hnum_plus->GetBinContent(ibin + 1);
 948       vale = hnum_plus->GetBinError(ibin + 1);
 949       val *= (1. + syse);
 950       vale *= (1. + syse);
 951       hnum_plus->SetBinContent(ibin + 1, val);
 952       hnum_plus->SetBinError(ibin + 1, vale);
 953
 954       val = hnum_minus->GetBinContent(ibin + 1);
 955       vale = hnum_minus->GetBinError(ibin + 1);
 956       val *= (1. - syse);
 957       vale *= (1. - syse);
 958       hnum_minus->SetBinContent(ibin + 1, val);
 959       hnum_minus->SetBinError(ibin + 1, vale);
 960
 961       val = hden_plus->GetBinContent(ibin + 1);
 962       vale = hden_plus->GetBinError(ibin + 1);
 963       val *= (1. + syse);
 964       vale *= (1. + syse);
 965       hden_plus->SetBinContent(ibin + 1, val);
 966       hden_plus->SetBinError(ibin + 1, vale);
 967
 968       val = hden_minus->GetBinContent(ibin + 1);
 969       vale = hden_minus->GetBinError(ibin + 1);
 970       val *= (1. - syse);
 971       vale *= (1. - syse);
 972       hden_minus->SetBinContent(ibin + 1, val);
 973       hden_minus->SetBinError(ibin + 1, vale);
 974
 975
 976
 977       val = hnum_stat_plus->GetBinContent(ibin + 1);
 978       vale = hnum_stat_plus->GetBinError(ibin + 1);
 979       val *= (1. + syse);
 980       vale *= (1. + syse);
 981       hnum_stat_plus->SetBinContent(ibin + 1, val);
 982       hnum_stat_plus->SetBinError(ibin + 1, vale);
 983
 984       val = hnum_stat_minus->GetBinContent(ibin + 1);
 985       vale = hnum_stat_minus->GetBinError(ibin + 1);
 986       val *= (1. - syse);
 987       vale *= (1. - syse);
 988       hnum_stat_minus->SetBinContent(ibin + 1, val);
 989       hnum_stat_minus->SetBinError(ibin + 1, vale);
 990
 991       val = hden_stat_plus->GetBinContent(ibin + 1);
 992       vale = hden_stat_plus->GetBinError(ibin + 1);
 993       val *= (1. + syse);
 994       vale *= (1. + syse);
 995       hden_stat_plus->SetBinContent(ibin + 1, val);
 996       hden_stat_plus->SetBinError(ibin + 1, vale);
 997
 998       val = hden_stat_minus->GetBinContent(ibin + 1);
 999       vale = hden_stat_minus->GetBinError(ibin + 1);
1000       val *= (1. - syse);
1001       vale *= (1. - syse);
1002       hden_stat_minus->SetBinContent(ibin + 1, val);
1003       hden_stat_minus->SetBinError(ibin + 1, vale);
1004     }
1005
1006     TCanvas *c = new TCanvas;
1007     hnum_stat->SetLineColor(8);
1008     hnum_stat->SetMarkerColor(8);
1009     hnum_stat->SetMarkerSize(1);
1010     hnum_stat->SetMarkerStyle(20);
1011     hnum_stat->Draw();
1012
1013     hnum_stat_plus->SetLineColor(2);
1014     hnum_stat_plus->SetMarkerColor(2);
1015     hnum_stat_plus->SetMarkerSize(1);
1016     hnum_stat_plus->SetMarkerStyle(20);
1017     hnum_stat_plus->Draw("same");
1018
1019     hnum_stat_minus->SetLineColor(4);
1020     hnum_stat_minus->SetMarkerColor(4);
1021     hnum_stat_minus->SetMarkerSize(1);
1022     hnum_stat_minus->SetMarkerStyle(20);
1023     hnum_stat_minus->Draw("same");
1024
1025     c->Update();
1026
1027     c = new TCanvas();
1028     TH1 *hrat = (TH1 *)hnum_stat->Clone("hrat");
1029     TH1 *hrat_plus = (TH1 *)hnum_stat_plus->Clone("hrat_plus");
1030     TH1 *hrat_minus = (TH1 *)hnum_stat_minus->Clone("hrat_minus");
1031
1032     hrat->Divide(hden_stat);
1033     hrat_plus->Divide(hden_stat);
1034     hrat_minus->Divide(hden_stat);
1035
1036     hrat->Draw();
1037     hrat_plus->Draw("same");
1038     hrat_minus->Draw("same");
1039
1040
1041     c->Update();
1042
1043     getchar();
1044
1045
1046     /* integrate as they are */
1047     printf("***** NUM:\n");
1048     IntegratedProduction(hnum, AliPID::ParticleMass(ratioPartNum[irat]), 0, "0q");
1049     printf("***** NUM PLUS:\n");
1050     IntegratedProduction(hnum_plus, AliPID::ParticleMass(ratioPartNum[irat]), 0, "0q");
1051     printf("***** NUM MINUS:\n");
1052     IntegratedProduction(hnum_minus, AliPID::ParticleMass(ratioPartNum[irat]), 0, "0q");
1053
1054     printf("***** DEN:\n");
1055     IntegratedProduction(hden, AliPID::ParticleMass(ratioPartDen[irat]), 0, "0q");
1056     printf("***** DEN PLUS:\n");
1057     IntegratedProduction(hden_plus, AliPID::ParticleMass(ratioPartDen[irat]), 0, "0q");
1058     printf("***** DEN MINUS:\n");
1059     IntegratedProduction(hden_minus, AliPID::ParticleMass(ratioPartDen[irat]), 0, "0q");
1060
1061   }
1062
1063 }
1064
1065 enum EIntegratedData_t {
1066   kInt_YieldData, kInt_YieldDataErr, kInt_YieldDataLin,
1067   kInt_MeanData, kInt_MeanDataErr, kInt_MeanDataLin,
1068   kInt_YieldLow, kInt_YieldLowErr,
1069   kInt_MeanLow, kInt_MeanLowErr,
1070   kInt_YieldHigh, kInt_YieldHighErr,
1071   kInt_MeanHigh, kInt_MeanHighErr,
1072   kInt_NData
1073 };
1074
1075 IntegratedProduction_measurement(const Char_t *filename, Option_t *opt = "q0R")
1076 {
1077
1078   for (Int_t ipart = 2; ipart < 5; ipart++)
1079     for (Int_t icharge = 0; icharge < 2; icharge++) {
1080       IntegratedProduction(filename, ipart, icharge, 0, 0., 10., opt);
1081       IntegratedProduction(filename, ipart, icharge, 1, 0., 10., opt);
1082     }
1083
1084 }
1085
1086 IntegratedProduction_systematics(const Char_t *filename, Option_t *opt = "q0R")
1087 {
1088
1089   Double_t min[5] = {0., 0., 0., 0., 0.};
1090   Double_t max[5] = {0., 0., 0.5, 1.0, 2.0};
1091   for (Int_t ipart = 2; ipart < 5; ipart++)
1092     for (Int_t icharge = 0; icharge < 2; icharge++)
1093       for (Int_t ifunc = 1; ifunc < 6; ifunc++)
1094         IntegratedProduction(filename, ipart, icharge, ifunc, min[ipart], max[ipart], opt);
1095
1096 }
1097
1098 IntegratedProduction_check(const Char_t *filename, Option_t *opt = "q0R")
1099 {
1100
1101   for (Int_t ipart = 2; ipart < 5; ipart++)
1102     for (Int_t icharge = 0; icharge < 2; icharge++)
1103       for (Int_t ifunc = 1; ifunc < 6; ifunc++)
1104         IntegratedProduction(filename, ipart, icharge, ifunc, 0., 10., opt);
1105
1106 }
1107
1108 IntegratedProduction(const Char_t *filename, Int_t ipart, Int_t icharge, Int_t ifunc = 0, Float_t min = 0., Float_t max = 10., Option_t *opt = "q0R")
1109 {
1110
1111   Char_t *centName[10] = {
1112     "0-5%", "5-10%", "10-20%", "20-40%", "40-60%", "60-80%", "80-100%", "minimum bias"
1113   };
1114
1115
1116   Char_t *chargeName[3] = {"plus", "minus", "sum"};
1117
1118   Char_t *partChargeName[5][2] = {
1119     "", "",
1120     "", "",
1121     "#pi^{+}", "#pi^{-}",
1122     "K^{+}", "K^{-}",
1123     "p", "#bar{p}"
1124   };
1125
1126   gStyle->SetOptStat(0);
1127   gStyle->SetOptFit(1);
1128
1129   /* PWGTools */
1130   gSystem->Load("libANALYSIS");
1131   gSystem->Load("libANALYSISalice");
1132   gSystem->Load("libCORRFW");
1133   gSystem->Load("libPWGTools");
1134   AliPWGFunc pwgfunc;
1135   pwgfunc.SetVarType(AliPWGFunc::kdNdpt);
1136
1137   Double_t mass = AliPID::ParticleMass(ipart);
1138
1139   TF1 *f = NULL;
1140   switch (ifunc) {
1141   case 0:
1142     f = BGBlastWave("BlastWave", mass);
1143     f->SetLineColor(2);
1144     break;
1145   case 1:
1146     f = LevyTsallis("Levy-Tsallis", mass, 5., mass);
1147     f->SetLineColor(8);
1148     break;
1149   case 2:
1150     f = Boltzmann("Boltzmann", mass);
1151     f->SetLineColor(kPink+1);
1152     break;
1153   case 3:
1154     f = pwgfunc.GetMTExp(mass, 0.1, 1., "m_{T}-exponential");
1155     f->SetLineColor(kViolet+1);
1156     break;
1157   case 4:
1158     f = pwgfunc.GetPTExp(0.1, 1., "p_{T}-exponential");
1159     f->SetLineColor(kAzure+1);
1160     break;
1161   case 5:
1162     f = pwgfunc.GetBoseEinstein(mass, 0.1, 1., "Bose-Einstein");
1163     f->SetLineColor(kYellow+1);
1164     break;
1165   default:
1166   };
1167   f->SetTitle(f->GetName());
1168   f->SetLineWidth(2);
1169   f->SetFillColor(0);
1170   f->SetMarkerColor(f->GetLineColor());
1171
1172   /* define low-pt fit range according to the number of free parameters */
1173   Int_t nfreeparams = f->GetNumberFreeParameters();
1174
1175   TCanvas *cc = new TCanvas("cCanvasFit", "", 1200, 800);
1176   cc->Divide(4, 2);
1177
1178   //  if (max < 0)
1179   //  TFile *fout = TFile::Open(Form("%s.IntegratedProduction.LowPtFit_%s_%s_%s.root", filename, f->GetName(), AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
1180   //  else if (min < 0)
1181   //  TFile *fout = TFile::Open(Form("%s.IntegratedProduction.HighPtFit_%s_%s_%s.root", filename, f->GetName(), AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
1182   //  else
1183     TFile *fout = TFile::Open(Form("%s.IntegratedProduction_%s_%s_%s.root", filename, f->GetName(), AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
1184
1185   Double_t integrated_data[kInt_NData];
1186
1187   TFile *filein = TFile::Open(filename);
1188   TH1F *hyield_data_stat = new TH1F(Form("hPartYield_data_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1189   TH1F *hyield_data_sys = new TH1F(Form("hPartYield_data_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1190   TH1F *hyield_data = new TH1F(Form("hPartYield_data_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1191   TH1F *hyield_low_stat = new TH1F(Form("hPartYield_low_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1192   TH1F *hyield_low_sys = new TH1F(Form("hPartYield_low_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1193   TH1F *hyield_low = new TH1F(Form("hPartYield_low_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1194   TH1F *hyield_high_stat = new TH1F(Form("hPartYield_high_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1195   TH1F *hyield_high_sys = new TH1F(Form("hPartYield_high_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1196   TH1F *hyield_high = new TH1F(Form("hPartYield_high_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1197   TH1F *hmean_data_stat = new TH1F(Form("hPartMean_data_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1198   TH1F *hmean_data_sys = new TH1F(Form("hPartMean_data_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1199   TH1F *hmean_data = new TH1F(Form("hPartMean_data_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1200   TH1F *hmean_low_stat = new TH1F(Form("hPartMean_low_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1201   TH1F *hmean_low_sys = new TH1F(Form("hPartMean_low_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1202   TH1F *hmean_low = new TH1F(Form("hPartMean_low_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1203   TH1F *hmean_high_stat = new TH1F(Form("hPartMean_high_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1204   TH1F *hmean_high_sys = new TH1F(Form("hPartMean_high_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1205   TH1F *hmean_high = new TH1F(Form("hPartMean_high_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1206   TH1F *hyield_stat = new TH1F(Form("hYield_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1207   TH1F *hyield_sys = new TH1F(Form("hYield_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1208   TH1F *hyield = new TH1F(Form("hYield_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1209   TH1F *hmean_stat = new TH1F(Form("hMean_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1210   TH1F *hmean_sys = new TH1F(Form("hMean_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1211   TH1F *hmean = new TH1F(Form("hMean_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
1212   for (Int_t icent = 0; icent < 7; icent++) {
1213     cc->cd(icent + 1)->SetLogy();
1214     printf("*************************\n");
1215     printf("cent=%d part=%d charge=%d\n", icent, ipart, icharge);
1216     TH1 *h = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1217     //    h = Convert_dNdy_dNdeta(h, AliPID::ParticleMass(ipart), 0.5);
1218     h->SetTitle(Form("%s (%s);p_{T} (GeV/c);dN/dp_{T}", partChargeName[ipart][icharge], centName[icent]));
1219     h->SetMarkerStyle(20);
1220     h->SetMarkerSize(1);
1221     h->SetLineWidth(1);
1222     h->SetFillColor(0);
1223     h->SetFillStyle(0);
1224     h->SetLineColor(1);
1225     h->SetMarkerColor(4);
1226     TH1 *hstat = (TH1 *)filein->Get(Form("stat_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1227     //    hstat = Convert_dNdy_dNdeta(hstat, AliPID::ParticleMass(ipart), 0.5);
1228     hstat->SetMarkerStyle(20);
1229     hstat->SetMarkerSize(1);
1230     hstat->SetLineWidth(1);
1231     hstat->SetLineColor(1);
1232     hstat->SetMarkerColor(4);
1233     TH1 *hsys = (TH1 *)filein->Get(Form("sys_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1234     //    hsys = Convert_dNdy_dNdeta(hsys, AliPID::ParticleMass(ipart), 0.5);
1235     hsys->SetMarkerStyle(20);
1236     hsys->SetMarkerSize(1);
1237     hsys->SetLineWidth(1);
1238     hsys->SetLineColor(1);
1239     hsys->SetMarkerColor(4);
1240
1241     if (max < 0.) {
1242       /* define fit range according to number of free parameter */
1243       Int_t gotpoints = 0;
1244       for (Int_t ipt = 0; ipt < h->GetNbinsX(); ipt++) {
1245         if (h->GetBinError(ipt + 1) <= 0.) continue;
1246         gotpoints++;
1247         if (gotpoints < nfreeparams*2) continue;
1248         max = h->GetXaxis()->GetBinUpEdge(ipt + 1);
1249         break;
1250       }
1251       printf("fitting up to %f\n", max);
1252     }
1253     if (min < 0.) {
1254       /* define fit range according to number of free parameter */
1255       Int_t gotpoints = 0;
1256       for (Int_t ipt = h->GetNbinsX(); ipt > 0; ipt--) {
1257         if (h->GetBinError(ipt) <= 0.) continue;
1258         gotpoints++;
1259         if (gotpoints < nfreeparams+3) continue;
1260         min = h->GetXaxis()->GetBinLowEdge(ipt);
1261         break;
1262       }
1263       printf("fitting from %f\n", min);
1264     }
1265
1266     printf("-> processing stat-only data\n");
1267
1268     IntegratedProduction(hstat, f, opt, min, max, integrated_data);
1269     fout->cd();
1270     hstat->Write(Form("IntegratedProduction_stat_cent%d_%s_%s.root", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1271
1272     hyield_data_stat->SetBinContent(icent + 1, integrated_data[kInt_YieldData]);
1273     hyield_data_stat->SetBinError(icent + 1, integrated_data[kInt_YieldDataErr]);
1274     hyield_low_stat->SetBinContent(icent + 1, integrated_data[kInt_YieldLow]);
1275     hyield_low_stat->SetBinError(icent + 1, integrated_data[kInt_YieldLowErr]);
1276
1277     hyield_high_stat->SetBinContent(icent + 1, integrated_data[kInt_YieldHigh]);
1278     hyield_high_stat->SetBinError(icent + 1, integrated_data[kInt_YieldHighErr]);
1279
1280     hmean_data_stat->SetBinContent(icent + 1, integrated_data[kInt_MeanData]);
1281     hmean_data_stat->SetBinError(icent + 1, integrated_data[kInt_MeanDataErr]);
1282
1283     hmean_low_stat->SetBinContent(icent + 1, integrated_data[kInt_MeanLow]);
1284     hmean_low_stat->SetBinError(icent + 1, integrated_data[kInt_MeanLowErr]);
1285
1286     printf("-> processing sys-only data\n");
1287
1288     IntegratedProduction(hsys, f, opt, min, max, integrated_data);
1289     fout->cd();
1290     hsys->Write(Form("IntegratedProduction_sys_cent%d_%s_%s.root", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1291
1292     hyield_data_sys->SetBinContent(icent + 1, integrated_data[kInt_YieldData]);
1293     hyield_data_sys->SetBinError(icent + 1, integrated_data[kInt_YieldDataLin]);
1294
1295     hyield_low_sys->SetBinContent(icent + 1, integrated_data[kInt_YieldLow]);
1296     hyield_low_sys->SetBinError(icent + 1, integrated_data[kInt_YieldLowErr]);
1297
1298     hyield_high_sys->SetBinContent(icent + 1, integrated_data[kInt_YieldHigh]);
1299     hyield_high_sys->SetBinError(icent + 1, integrated_data[kInt_YieldHighErr]);
1300
1301     hmean_data_sys->SetBinContent(icent + 1, integrated_data[kInt_MeanData]);
1302     hmean_data_sys->SetBinError(icent + 1, integrated_data[kInt_MeanDataLin]);
1303
1304     hmean_low_sys->SetBinContent(icent + 1, integrated_data[kInt_MeanLow]);
1305     hmean_low_sys->SetBinError(icent + 1, integrated_data[kInt_MeanLowErr]);
1306
1307     printf("-> processing stat+sys data\n");
1308
1309     IntegratedProduction(h, f, opt, min, max, integrated_data);
1310     fout->cd();
1311     h->Write(Form("IntegratedProduction_cent%d_%s_%s.root", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1312     f->SetRange(min, max);
1313     f->Write(Form("%s_cent%d_%s_%s.root", f->GetName(), icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1314     h->DrawCopy();
1315     f->DrawCopy("same");
1316     TLegend *l = gPad->BuildLegend(0.15, 0.12, 0.5, 0.3);
1317     l->SetBorderSize(0);
1318     l->SetFillColor(0);
1319     l->SetFillStyle(0);
1320
1321     hyield_data->SetBinContent(icent + 1, integrated_data[kInt_YieldData]);
1322     hyield_data->SetBinError(icent + 1, integrated_data[kInt_YieldDataLin]);
1323
1324     hyield_low->SetBinContent(icent + 1, integrated_data[kInt_YieldLow]);
1325     hyield_low->SetBinError(icent + 1, integrated_data[kInt_YieldLowErr]);
1326
1327     hyield_high->SetBinContent(icent + 1, integrated_data[kInt_YieldHigh]);
1328     hyield_high->SetBinError(icent + 1, integrated_data[kInt_YieldHighErr]);
1329
1330     hmean_data->SetBinContent(icent + 1, integrated_data[kInt_MeanData]);
1331     hmean_data->SetBinError(icent + 1, integrated_data[kInt_MeanDataLin]);
1332
1333     hmean_low->SetBinContent(icent + 1, integrated_data[kInt_MeanLow]);
1334     hmean_low->SetBinError(icent + 1, integrated_data[kInt_MeanLowErr]);
1335
1336     hmean_high->SetBinContent(icent + 1, integrated_data[kInt_MeanHigh]);
1337     hmean_high->SetBinError(icent + 1, integrated_data[kInt_MeanHighErr]);
1338
1339
1340   }
1341
1342   hyield_stat->Add(hyield_data_stat);
1343   hyield_stat->Add(hyield_low_stat);
1344   hyield_stat->Add(hyield_high_stat);
1345
1346   hyield_sys->Add(hyield_data_sys);
1347   hyield_sys->Add(hyield_low_sys);
1348   hyield_sys->Add(hyield_high_sys);
1349
1350   hyield->Add(hyield_data);
1351   hyield->Add(hyield_low);
1352   hyield->Add(hyield_high);
1353
1354
1355   hmean_stat->Add(hmean_data_stat);
1356   hmean_stat->Add(hmean_low_stat);
1357   hmean_stat->Add(hmean_high_stat);
1358   //      hmean_stat->Divide(hyield_stat);
1359   hmean_sys->Add(hmean_data_sys);
1360   hmean_sys->Add(hmean_low_sys);
1361   hmean_sys->Add(hmean_high_sys);
1362   //      hmean_sys->Divide(hyield_sys);
1363   hmean->Add(hmean_data);
1364   hmean->Add(hmean_low);
1365   hmean->Add(hmean_high);
1366
1367   fout->cd();
1368
1369   hyield->Write();
1370   hmean->Write();
1371
1372   hyield_stat->Write();
1373   hmean_stat->Write();
1374
1375   hyield_sys->Write();
1376   hmean_sys->Write();
1377
1378   hyield_data->Write();
1379   hmean_data->Write();
1380
1381   hyield_data_stat->Write();
1382   hmean_data_stat->Write();
1383
1384   hyield_data_sys->Write();
1385   hmean_data_sys->Write();
1386
1387   hyield_low->Write();
1388   hmean_low->Write();
1389
1390   hyield_high->Write();
1391   hmean_high->Write();
1392
1393   hyield_low_stat->Write();
1394   hmean_low_stat->Write();
1395
1396   hyield_high_stat->Write();
1397   hmean_high_stat->Write();
1398
1399   hmean_low_sys->Write();
1400   hyield_low_sys->Write();
1401
1402   hmean_high_sys->Write();
1403   hyield_high_sys->Write();
1404
1405   cc->Write();
1406
1407   fout->Close();
1408 }
1409
1410 IntegratedProduction_pp(const Char_t *filename, Char_t *what = "", Int_t ifunc = 1, Option_t *opt = "0qI")
1411 {
1412   TFile *filein = TFile::Open(filename);
1413   for (Int_t ipart = 0; ipart < 6; ipart++) {
1414     printf("*************************\n");
1415     printf("part=%d\n", ipart);
1416     TH1 *h = (TH1 *)filein->Get(Form("hComb_ITSsa%d_TPC%d_TOF%d_HMPID%d", ipart, ipart, ipart, ipart));
1417     IntegratedProduction(h, AliPID::ParticleMass(ipart % 3 + 2), ifunc, opt);
1418     if (gPad) gPad->SaveAs(Form("IntegratedProduction_pp_%d.root", ipart));
1419   }
1420 }
1421
1422 IntegratedProduction(TH1 *h, TF1 *f = NULL, Option_t *opt = "0q", Float_t min = 0., Float_t max = 10., Double_t *integrated_data = NULL, Bool_t verbose = kFALSE)
1423 {
1424
1425   Double_t yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, partyield[3], partyielderr[3], partyielderrcorr[3], partmean[3], partmeanerr[3], partmeanerrcorr[3];
1426
1427   TVirtualFitter::SetMaxIterations(1000000);
1428
1429   if (f) {
1430     Int_t fres = 1;
1431     Int_t fatt = 0;
1432     while (fres != 0 && fatt < 10) {
1433       fres = h->Fit(f, opt, "", min, max);
1434       fres = h->Fit(f, opt, "", min, max);
1435       printf("fit res = %d\n", fres);
1436       fatt++;
1437     }
1438   }
1439
1440   // return;
1441
1442   GetYieldAndMean(h, f, yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, 0., 10., partyield, partyielderr, partyielderrcorr, partmean, partmeanerr, partmeanerrcorr);
1443
1444   //  Double_t fint = f ? f->Integral(0.,10.) : 0.;
1445   //  Double_t finte = f ? f->IntegralError(0.,10.) : 0.;
1446   //  Double_t fmean = f ? f->Mean(0., 10.) : 0.;
1447
1448   if (verbose) {
1449   printf("----\n");
1450   printf("dN/dy        = %f +- %f (%f)\n", yield, yielderr, yielderrcorr);
1451   printf("<pt>         = %f +- %f (%f)\n", mean, meanerr, meanerrcorr);
1452   printf("----\n");
1453   printf("dN/dy (data) = %f +- %f (%f)\n", partyield[0], partyielderr[0], partyielderrcorr[0]);
1454   printf("dN/dy (low)  = %f +- %f (%f)\n", partyield[1], partyielderr[1], partyielderrcorr[1]);
1455   printf("dN/dy (high) = %f +- %f (%f)\n", partyield[2], partyielderr[2], partyielderrcorr[2]);
1456   printf("<pt> (data)  = %f +- %f (%f)\n", partmean[0], partmeanerr[0], partmeanerrcorr[0]);
1457   printf("<pt> (low)   = %f +- %f (%f)\n", partmean[1], partmeanerr[1], partmeanerrcorr[1]);
1458   printf("<pt> (high)  = %f +- %f (%f)\n", partmean[2], partmeanerr[2], partmeanerrcorr[2]);
1459   printf("----\n");
1460   //  printf("dN/dy (func) = %f +- %f\n", fint, finte);
1461   //  printf("<pT> (func)  = %f +- %f\n", fmean, 0.);
1462   }
1463
1464   if (!integrated_data) return;
1465
1466   integrated_data[kInt_YieldData] = partyield[0];
1467   integrated_data[kInt_YieldDataErr] = partyielderr[0];
1468   integrated_data[kInt_YieldDataLin] = partyielderrcorr[0];
1469   integrated_data[kInt_MeanData] = partmean[0];
1470   integrated_data[kInt_MeanDataErr] = partmeanerr[0];
1471   integrated_data[kInt_MeanDataLin] = partmeanerrcorr[0];
1472   integrated_data[kInt_YieldLow] = partyield[1];
1473   integrated_data[kInt_YieldLowErr] = partyielderr[1];
1474   integrated_data[kInt_MeanLow] = partmean[1];
1475   integrated_data[kInt_MeanLowErr] = partmeanerr[1];
1476   integrated_data[kInt_YieldHigh] = partyield[2];
1477   integrated_data[kInt_YieldHighErr] = partyielderr[2];
1478   integrated_data[kInt_MeanHigh] = partmean[2];
1479   integrated_data[kInt_MeanHighErr] = partmeanerr[2];
1480
1481   //  TH1 *hr = (TH1 *)h->Clone("hr");
1482   //  hr->Divide(f);
1483   //  new TCanvas;
1484   //  hr->Draw();
1485
1486   //  TProfile *p = new TProfile("p", "", 100, 0., 10.);
1487   //  gROOT->LoadMacro("HistoUtils.C");
1488   //  HistoUtils_Function2Profile(f, p);
1489   //  p->Draw();
1490 }
1491
1492 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, Double_t *partmean, Double_t *partmeanerr, Double_t *partmeanerrcorr)
1493 {
1494
1495   /* find lowest edge in histo */
1496   Int_t binlo;
1497   Double_t lo;
1498   for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
1499     if (h->GetBinContent(ibin) != 0.) {
1500       binlo = ibin;
1501       lo = h->GetBinLowEdge(ibin);
1502       break;
1503     }
1504   }
1505
1506   /* find highest edge in histo */
1507   Int_t binhi;
1508   Double_t hi;
1509   for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
1510     if (h->GetBinContent(ibin) != 0.) {
1511       binhi = ibin + 1;
1512       hi = h->GetBinLowEdge(ibin + 1);
1513       break;
1514     }
1515   }
1516
1517   /* integrate the data */
1518   Double_t cont, err, width, cent, integral_data = 0., integralerr_data = 0., integralerrcorr_data = 0., meanintegral_data = 0., meanintegralerr_data = 0., meanintegralerrcorr_data = 0.;
1519   for (Int_t ibin = binlo; ibin < binhi; ibin++) {
1520     cent = h->GetBinCenter(ibin);
1521     width = h->GetBinWidth(ibin);
1522     cont = h->GetBinContent(ibin);
1523     err = h->GetBinError(ibin);
1524     /* check we didn't get an empty bin in between */
1525     if (cont != 0. && err != 0.) {
1526       /* all right, use data */
1527       integral_data += cont * width;
1528       integralerr_data += err * err * width * width;
1529       integralerrcorr_data += err * width;
1530       meanintegral_data += cont * width * cent;
1531       meanintegralerr_data += err * err * width * width * cent * cent;
1532       meanintegralerrcorr_data += err * width * cent;
1533     }
1534     else {
1535       /* missing data-point, complain and use function */
1536       printf("WARNING: missing data-point at %f\n", cent);
1537       printf("         using function as a patch\n");
1538       integral_data += f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
1539       integralerr_data += f->IntegralError(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1), 0, 0, 1.e-6);
1540       meanintegral_data += f->Mean(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1)) * f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
1541       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);
1542     }
1543   }
1544   integralerr_data = TMath::Sqrt(integralerr_data);
1545   meanintegralerr_data = TMath::Sqrt(meanintegralerr_data);
1546
1547   /* integrate below the data */
1548   if (!f) {
1549     min = lo;
1550     max = hi;
1551     printf("WARNING: no function provided! Use only data with no extrapolation\n");
1552   }
1553   Double_t integral_lo = min < lo ? f->Integral(min, lo, (Double_t *)0, 1.e-6) : 0.;
1554   Double_t integralerr_lo = min < lo ? f->IntegralError(min, lo, 0, 0, 1.e-6) : 0.;
1555   Double_t meanintegral_lo = min < lo ? f->Mean(min, lo, (Double_t *)0, 1.e-6) * integral_lo : 0.;
1556   Double_t meanintegralerr_lo = min < lo ? f->Mean(min, lo, (Double_t *)0, 1.e-6) * integralerr_lo : 0.;
1557
1558   /* integrate above the data */
1559   Double_t integral_hi = max > hi ? f->Integral(hi, max, (Double_t *)0, 1.e-6) : 0.;
1560   Double_t integralerr_hi = max > hi ? f->IntegralError(hi, max, 0, 0, 1.e-6) : 0.;
1561   Double_t meanintegral_hi = max > hi ? f->Mean(hi, max, (Double_t *)0, 1.e-6) * integral_hi : 0.;
1562   Double_t meanintegralerr_hi = max > hi ? f->Mean(hi, max, (Double_t *)0, 1.e-6) * integralerr_hi : 0.;
1563
1564   /* compute integrated yield */
1565   yield = integral_data + integral_lo + integral_hi;
1566   yielderr = TMath::Sqrt(integralerr_data * integralerr_data +
1567                          integralerr_lo * integralerr_lo +
1568                          integralerr_hi * integralerr_hi);
1569   yielderrcorr = TMath::Sqrt(integralerrcorr_data * integralerrcorr_data +
1570                              integralerr_lo * integralerr_lo +
1571                              integralerr_hi * integralerr_hi);
1572
1573   /* compute integrated mean */
1574   mean = (meanintegral_data + meanintegral_lo + meanintegral_hi) / yield;
1575   meanerr = TMath::Sqrt(meanintegralerr_data * meanintegralerr_data +
1576                         meanintegralerr_lo * meanintegralerr_lo +
1577                         meanintegralerr_hi * meanintegralerr_hi) / yield;
1578   meanerrcorr = TMath::Sqrt(meanintegralerrcorr_data * meanintegralerrcorr_data +
1579                             meanintegralerr_lo * meanintegralerr_lo +
1580                             meanintegralerr_hi * meanintegralerr_hi) / yield;
1581
1582   /* set partial yields */
1583   partyield[0] = integral_data;
1584   partyielderr[0] = integralerr_data;
1585   partyielderrcorr[0] = integralerrcorr_data;
1586   partyield[1] = integral_lo;
1587   partyielderr[1] = integralerr_lo;
1588   partyielderrcorr[1] = integralerr_lo;
1589   partyield[2] = integral_hi;
1590   partyielderr[2] = integralerr_hi;
1591   partyielderrcorr[2] = integralerr_hi;
1592
1593   /* set partial means */
1594   partmean[0] = meanintegral_data;
1595   partmeanerr[0] = meanintegralerr_data;
1596   partmeanerrcorr[0] = meanintegralerrcorr_data;
1597   partmean[1] = meanintegral_lo;
1598   partmeanerr[1] = meanintegralerr_lo;
1599   partmeanerrcorr[1] = meanintegralerr_lo;
1600   partmean[2] = meanintegral_hi;
1601   partmeanerr[2] = meanintegralerr_hi;
1602   partmeanerrcorr[2] = meanintegralerr_hi;
1603
1604 }
1605
1606 /*****************************************************************/
1607
1608 Double_t
1609 y2eta(Double_t pt, Double_t mass, Double_t y){
1610   Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
1611   return TMath::ASinH(mt / pt * TMath::SinH(y));
1612 }
1613 Double_t
1614 eta2y(Double_t pt, Double_t mass, Double_t eta){
1615   Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
1616   return TMath::ASinH(pt / mt * TMath::SinH(eta));
1617 }
1618
1619 TH1 *
1620 Convert_dNdy_1over2pipt_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
1621 {
1622
1623   TH1 *hout = hin->Clone("hout");
1624   hout->Reset();
1625   Double_t pt, mt, conv, val, vale;
1626   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
1627     pt = hin->GetBinCenter(ibin + 1);
1628     conv = eta2y(pt, mass, eta) / eta;
1629     val = hin->GetBinContent(ibin + 1);
1630     vale = hin->GetBinError(ibin + 1);
1631     val /= (2. * TMath::Pi() * pt);
1632     vale /= (2. * TMath::Pi() * pt);
1633     val *= conv;
1634     vale *= conv;
1635     hout->SetBinContent(ibin + 1, val);
1636     hout->SetBinError(ibin + 1, vale);
1637   }
1638   return hout;
1639 }
1640
1641 TH1 *
1642 Convert_dNdy_1over2pipt_dNdy(TH1 *hin)
1643 {
1644
1645   TH1 *hout = hin->Clone("hout");
1646   hout->Reset();
1647   Double_t pt, mt, conv, val, vale;
1648   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
1649     pt = hin->GetBinCenter(ibin + 1);
1650     val = hin->GetBinContent(ibin + 1);
1651     vale = hin->GetBinError(ibin + 1);
1652     val /= (2. * TMath::Pi() * pt);
1653     vale /= (2. * TMath::Pi() * pt);
1654     hout->SetBinContent(ibin + 1, val);
1655     hout->SetBinError(ibin + 1, vale);
1656   }
1657   return hout;
1658 }
1659
1660 TH1 *
1661 Convert_dNdy_1overpt_dNdy(TH1 *hin)
1662 {
1663
1664   TH1 *hout = hin->Clone("hout");
1665   hout->Reset();
1666   Double_t pt, mt, conv, val, vale;
1667   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
1668     pt = hin->GetBinCenter(ibin + 1);
1669     val = hin->GetBinContent(ibin + 1);
1670     vale = hin->GetBinError(ibin + 1);
1671     val /= pt;
1672     vale /= pt;
1673     hout->SetBinContent(ibin + 1, val);
1674     hout->SetBinError(ibin + 1, vale);
1675   }
1676   return hout;
1677 }
1678
1679 TH1 *
1680 Convert_dNdy_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
1681 {
1682
1683   TH1 *hout = hin->Clone("hout");
1684   hout->Reset();
1685   Double_t pt, mt, conv, val, vale;
1686   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
1687     pt = hin->GetBinCenter(ibin + 1);
1688     conv = eta2y(pt, mass, eta) / eta;
1689     val = hin->GetBinContent(ibin + 1);
1690     vale = hin->GetBinError(ibin + 1);
1691     val *= conv;
1692     vale *= conv;
1693     hout->SetBinContent(ibin + 1, val);
1694     hout->SetBinError(ibin + 1, vale);
1695   }
1696   return hout;
1697 }
1698
1699 TGraph *
1700 Convert_dNdy_dNdeta(TGraph *hin, Double_t mass, Double_t eta = 0.8)
1701 {
1702
1703   TGraph *hout = hin->Clone("hout");
1704   //  hout->Reset();
1705   Double_t pt, mt, conv, val, valelo, valehi;
1706   for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
1707     pt = hin->GetX()[ibin];
1708     conv = eta2y(pt, mass, eta) / eta;
1709     val = hin->GetY()[ibin];
1710     valelo = hin->GetEYlow()[ibin];
1711     valehi = hin->GetEYhigh()[ibin];
1712     val *= conv;
1713     valelo *= conv;
1714     valehi *= conv;
1715     hout->GetX()[ibin] = pt;
1716     hout->GetY()[ibin] = val;
1717     hout->GetEYlow()[ibin] = valelo;
1718     hout->GetEYhigh()[ibin] = valehi;
1719   }
1720   return hout;
1721 }
1722
1723 TH1 *
1724 SummedId_1over2pipt_dNdeta(const Char_t *filename, Int_t icent, Float_t etarange, Float_t scalePi = 1.)
1725 {
1726
1727   const Char_t *chargeName[2] = {
1728     "plus", "minus"
1729   };
1730
1731   TFile *filein = TFile::Open(filename);
1732   TH1 *hy[AliPID::kSPECIES][2];
1733   TH1 *heta[AliPID::kSPECIES][2];
1734   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
1735     for (Int_t icharge = 0; icharge < 2; icharge++) {
1736       hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1737       if (!hy[ipart][icharge]) {
1738         printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
1739         return NULL;
1740       }
1741       heta[ipart][icharge] = Convert_dNdy_1over2pipt_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart), etarange);
1742       if (ipart == 2) heta[ipart][icharge]->Scale(scalePi);
1743     }
1744
1745   /* sum */
1746   TH1D *hsum = heta[2][0]->Clone("hsum");
1747   hsum->Reset();
1748   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
1749     for (Int_t icharge = 0; icharge < 2; icharge++) {
1750       hsum->Add(heta[ipart][icharge]);
1751     }
1752   for (Int_t ipt = 0; ipt < hsum->GetNbinsX(); ipt++) {
1753     Double_t err = 0.;
1754     for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
1755       for (Int_t icharge = 0; icharge < 2; icharge++) {
1756         err += heta[ipart][icharge]->GetBinError(ipt + 1);
1757       }
1758     hsum->SetBinError(ipt + 1, err);
1759   }
1760   return hsum;
1761 }
1762
1763 TH1 *
1764 SummedId_dNdeta(const Char_t *filename, Int_t icent)
1765 {
1766
1767   const Char_t *chargeName[2] = {
1768     "plus", "minus"
1769   };
1770
1771   TFile *filein = TFile::Open(filename);
1772   TH1 *hy[AliPID::kSPECIES][2];
1773   TH1 *heta[AliPID::kSPECIES][2];
1774   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
1775     for (Int_t icharge = 0; icharge < 2; icharge++) {
1776       hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1777       if (!hy[ipart][icharge]) {
1778         printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
1779         return NULL;
1780       }
1781       heta[ipart][icharge] = Convert_dNdy_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
1782     }
1783
1784   /* sum */
1785   TH1D *hsum = heta[2][0]->Clone("hsum");
1786   hsum->Reset();
1787   for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
1788     for (Int_t icharge = 0; icharge < 2; icharge++)
1789       hsum->Add(heta[ipart][icharge]);
1790
1791   return hsum;
1792 }
1793
1794 /*****************************************************************/
1795
1796 TH1 *
1797 ReturnExtremeHighHisto(TH1 *hin, TH1 *herr = NULL)
1798 {
1799   TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehigh", hin->GetName()));
1800   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
1801     if (hin->GetBinError(ibin + 1) <= 0.) continue;
1802     Double_t val = hin->GetBinContent(ibin + 1);
1803     Double_t err = hin->GetBinError(ibin + 1);
1804     if (herr) err = herr->GetBinError(ibin + 1);
1805     hout->SetBinContent(ibin + 1, val + err);
1806   }
1807   return hout;
1808 }
1809
1810 TH1 *
1811 ReturnExtremeLowHisto(TH1 *hin, TH1 *herr = NULL)
1812 {
1813   TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremelow", hin->GetName()));
1814   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
1815     if (hin->GetBinError(ibin + 1) <= 0.) continue;
1816     Double_t val = hin->GetBinContent(ibin + 1);
1817     Double_t err = hin->GetBinError(ibin + 1);
1818     if (herr) err = herr->GetBinError(ibin + 1);
1819     hout->SetBinContent(ibin + 1, val - err);
1820   }
1821   return hout;
1822 }
1823
1824 TH1 *
1825 ReturnExtremeSoftHisto(TH1 *hin, TH1 *herr = NULL)
1826 {
1827   return ReturnExtremeHisto(hin, herr, -1.);
1828 }
1829
1830 TH1 *
1831 ReturnExtremeHardHisto(TH1 *hin, TH1 *herr = NULL)
1832 {
1833   return ReturnExtremeHisto(hin, herr, 1.);
1834 }
1835
1836 TH1 *
1837 ReturnExtremeHisto(TH1 *hin, TH1 *herr = NULL, Float_t sign = 1.)
1838 {
1839   Double_t ptlow, pthigh;
1840   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
1841     if (hin->GetBinError(ibin + 1) <= 0.) continue;
1842     ptlow = hin->GetBinLowEdge(ibin + 1);
1843     break;
1844   }
1845   for (Int_t ibin = hin->GetNbinsX(); ibin >= 0; ibin--) {
1846     if (hin->GetBinError(ibin + 1) <= 0.) continue;
1847     pthigh = hin->GetBinLowEdge(ibin + 2);
1848     break;
1849   }
1850
1851   Double_t mean = hin->GetMean();
1852   Double_t maxdiff = 0.;
1853   TH1 *hmax = NULL;
1854   for (Int_t inode = 0; inode < hin->GetNbinsX(); inode++) {
1855
1856     Double_t ptnode = hin->GetBinCenter(inode + 1);
1857     TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehard", hin->GetName()));
1858
1859     for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
1860       if (hin->GetBinError(ibin + 1) <= 0.) continue;
1861       Double_t val = hin->GetBinContent(ibin + 1);
1862       Double_t err = hin->GetBinError(ibin + 1);
1863       if (herr) err = herr->GetBinError(ibin + 1);
1864       Double_t cen = hin->GetBinCenter(ibin + 1);
1865       //    err *= -1. - (cen - ptlow) * (-1. - 1.) / (pthigh - ptlow);
1866       if (cen < ptnode)
1867         err *= -1. + (cen - ptlow) / (ptnode - ptlow);
1868       else
1869         err *= (cen - ptnode) / (pthigh - ptnode);
1870
1871       hout->SetBinContent(ibin + 1, val + sign * err);
1872     }
1873
1874     Double_t diff = TMath::Abs(mean - hout->GetMean());
1875     if (diff > maxdiff) {
1876       //      printf("found max at %f\n", ptnode);
1877       if (hmax) delete hmax;
1878       hmax = (TH1 *)hout->Clone("hmax");
1879       maxdiff = diff;
1880     }
1881     delete hout;
1882   }
1883   return hmax;
1884 }
1885
1886 TGraphErrors *
1887 ReturnExtremeSoftGraph(TGraphErrors *hin, TGraphErrors *herr = NULL)
1888 {
1889   TGraphErrors *hout = (TGraphErrors *)hin->Clone(Form("%s_extremesoft", hin->GetName()));
1890   Double_t ptlow, pthigh;
1891   Double_t ptlow = hin->GetX()[0];
1892   Double_t pthigh = hin->GetX()[hin->GetN() - 1];
1893
1894   for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
1895     Double_t val = hin->GetY()[ibin];
1896     Double_t err = hin->GetEY()[ibin];
1897     if (herr) err = herr->GetEY()[ibin];
1898     Double_t cen = hin->GetX()[ibin];
1899     err *= 1. + (cen - ptlow) * (-1. - 1.) / (pthigh - ptlow);
1900     hout->SetPoint(ibin, cen, val + err);
1901   }
1902   return hout;
1903 }
1904
1905
1906 TGraphErrors *
1907 ReturnExtremeHardGraph(TGraphErrors *hin, TGraphErrors *herr = NULL)
1908 {
1909   TGraphErrors *hout = (TGraphErrors *)hin->Clone(Form("%s_extremehard", hin->GetName()));
1910   Double_t ptlow, pthigh;
1911   Double_t ptlow = hin->GetX()[0];
1912   Double_t pthigh = hin->GetX()[hin->GetN() - 1];
1913
1914   for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
1915     Double_t val = hin->GetY()[ibin];
1916     Double_t err = hin->GetEY()[ibin];
1917     if (herr) err = herr->GetEY()[ibin];
1918     Double_t cen = hin->GetX()[ibin];
1919     err *= -1. - (cen - ptlow) * (-1. - 1.) / (pthigh - ptlow);
1920     hout->SetPoint(ibin, cen, val + err);
1921   }
1922   return hout;
1923 }
1924