added UTILS dir with first macros
authorrpreghen <rpreghen@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 13 Nov 2013 15:31:27 +0000 (15:31 +0000)
committerrpreghen <rpreghen@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 13 Nov 2013 15:31:27 +0000 (15:31 +0000)
PWGLF/SPECTRA/UTILS/SpectraUtils.C [new file with mode: 0644]
PWGLF/SPECTRA/UTILS/YieldMean.C [new file with mode: 0644]

diff --git a/PWGLF/SPECTRA/UTILS/SpectraUtils.C b/PWGLF/SPECTRA/UTILS/SpectraUtils.C
new file mode 100644 (file)
index 0000000..a6f3eb0
--- /dev/null
@@ -0,0 +1,1924 @@
+/*
+  several function used for PbPb combined spectra
+  Blast Wave is also implemented here
+  further documentation will come
+  
+  author: Roberto Preghenella
+  email : preghenella@bo.infn.it
+*/
+
+
+/*****************************************************************/
+/* BOLTZMANN
+/*****************************************************************/
+
+Double_t
+Boltzmann_Func(const Double_t *x, const Double_t *p)
+{
+  /* dN/dpt */
+  
+  Double_t pt = x[0];
+  Double_t mass = p[0];
+  Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
+  Double_t T = p[1];
+  Double_t norm = p[2];
+
+  return pt * norm * mt * TMath::Exp(-mt / T);
+}
+
+TF1 *
+Boltzmann(const Char_t *name, Double_t mass, Double_t T = 0.1, Double_t norm = 1.)
+{
+  
+  TF1 *fBoltzmann = new TF1(name, Boltzmann_Func, 0., 10., 3);
+  fBoltzmann->SetParameters(mass, T, norm);
+  fBoltzmann->SetParNames("mass", "T", "norm");
+  fBoltzmann->FixParameter(0, mass);
+  return fBoltzmann;
+}
+
+/*****************************************************************/
+/* LEVY-TSALLIS */
+/*****************************************************************/
+
+Double_t
+LevyTsallis_Func(const Double_t *x, const Double_t *p)
+{
+  /* dN/dpt */
+  
+  Double_t pt = x[0];
+  Double_t mass = p[0];
+  Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
+  Double_t n = p[1];
+  Double_t C = p[2];
+  Double_t norm = p[3];
+
+  Double_t part1 = (n - 1.) * (n - 2.);
+  Double_t part2 = n * C * (n * C + mass * (n - 2.));
+  Double_t part3 = part1 / part2;
+  Double_t part4 = 1. + (mt - mass) / n / C;
+  Double_t part5 = TMath::Power(part4, -n);
+  return pt * norm * part3 * part5;
+}
+
+TF1 *
+LevyTsallis(const Char_t *name, Double_t mass, Double_t n = 5., Double_t C = 0.1, Double_t norm = 1.)
+{
+  
+  TF1 *fLevyTsallis = new TF1(name, LevyTsallis_Func, 0., 10., 4);
+  fLevyTsallis->SetParameters(mass, n, C, norm);
+  fLevyTsallis->SetParNames("mass", "n", "C", "norm");
+  fLevyTsallis->FixParameter(0, mass);
+  fLevyTsallis->SetParLimits(1, 1.e-3, 1.e3);
+  fLevyTsallis->SetParLimits(2, 1.e-3, 1.e3);
+  fLevyTsallis->SetParLimits(3, 1.e-6, 1.e6);
+  return fLevyTsallis;
+}
+
+/*****************************************************************/
+/* BOLTZMANN-GIBBS BLAST-WAVE */
+/*****************************************************************/
+
+static TF1 *fBGBlastWave_Integrand = NULL;
+static TF1 *fBGBlastWave_Integrand_num = NULL;
+static TF1 *fBGBlastWave_Integrand_den = NULL;
+Double_t
+BGBlastWave_Integrand(const Double_t *x, const Double_t *p)
+{
+  
+  /* 
+     x[0] -> r (radius)
+     p[0] -> mT (transverse mass)
+     p[1] -> pT (transverse momentum)
+     p[2] -> beta_max (surface velocity)
+     p[3] -> T (freezout temperature)
+     p[4] -> n (velocity profile)
+  */
+  
+  Double_t r = x[0];
+  Double_t mt = p[0];
+  Double_t pt = p[1];
+  Double_t beta_max = p[2];
+  Double_t temp_1 = 1. / p[3];
+  Double_t n = p[4];
+
+  Double_t beta = beta_max * TMath::Power(r, n);
+  if (beta > 0.9999999999999999) beta = 0.9999999999999999;
+  Double_t rho = TMath::ATanH(beta);
+  Double_t argI0 = pt * TMath::SinH(rho) * temp_1;
+  if (argI0 > 700.) argI0 = 700.;
+  Double_t argK1 = mt * TMath::CosH(rho) * temp_1;
+  //  if (argI0 > 100 || argI0 < -100)
+  //    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);
+  return r * mt * TMath::BesselI0(argI0) * TMath::BesselK1(argK1);
+  
+}
+
+Double_t
+BGBlastWave_Func(const Double_t *x, const Double_t *p)
+{
+  /* dN/dpt */
+  
+  Double_t pt = x[0];
+  Double_t mass = p[0];
+  Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
+  Double_t beta_max = p[1];
+  Double_t temp = p[2];
+  Double_t n = p[3];
+  Double_t norm = p[4];
+  
+  if (!fBGBlastWave_Integrand)
+    fBGBlastWave_Integrand = new TF1("fBGBlastWave_Integrand", BGBlastWave_Integrand, 0., 1., 5);
+  fBGBlastWave_Integrand->SetParameters(mt, pt, beta_max, temp, n);
+  Double_t integral = fBGBlastWave_Integrand->Integral(0., 1., (Double_t *)0, 1.e-6);
+  return norm * pt * integral;
+}
+
+Double_t
+BGBlastWaveRatio_Func(const Double_t *x, const Double_t *p)
+{
+  /* dN/dpt */
+  
+  Double_t pt = x[0];
+  Double_t mass = p[0];
+  Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
+  Double_t beta_max_num = p[1];
+  Double_t temp_num = p[2];
+  Double_t n_num = p[3];
+  Double_t norm_num = p[4];
+  Double_t beta_max_den = p[5];
+  Double_t temp_den = p[6];
+  Double_t n_den = p[7];
+  Double_t norm_den = p[8];
+  
+  if (!fBGBlastWave_Integrand_num)
+    fBGBlastWave_Integrand_num = new TF1("fBGBlastWave_Integrand_num", BGBlastWave_Integrand, 0., 1., 5);
+  fBGBlastWave_Integrand_num->SetParameters(mt, pt, beta_max_num, temp_num, n_num);
+  Double_t integral_num = fBGBlastWave_Integrand_num->Integral(0., 1.);
+
+  if (!fBGBlastWave_Integrand_den)
+    fBGBlastWave_Integrand_den = new TF1("fBGBlastWave_Integrand_den", BGBlastWave_Integrand, 0., 1., 5);
+  fBGBlastWave_Integrand_den->SetParameters(mt, pt, beta_max_den, temp_den, n_den);
+  Double_t integral_den = fBGBlastWave_Integrand_den->Integral(0., 1.);
+
+  return (norm_num / norm_den) * (integral_num / integral_den);
+}
+
+Double_t
+BGBlastWaveParticleRatio_Func(const Double_t *x, const Double_t *p)
+{
+  /* dN/dpt */
+  
+  Double_t pt = x[0];
+  Double_t mass_num = p[0];
+  Double_t mass_den = p[1];
+  Double_t mt_num = TMath::Sqrt(pt * pt + mass_num * mass_num);
+  Double_t mt_den = TMath::Sqrt(pt * pt + mass_den * mass_den);
+  Double_t beta_max = p[2];
+  Double_t temp = p[3];
+  Double_t n = p[4];
+  Double_t norm_num = p[5];
+  Double_t norm_den = p[6];
+  
+  if (!fBGBlastWave_Integrand_num)
+    fBGBlastWave_Integrand_num = new TF1("fBGBlastWave_Integrand_num", BGBlastWave_Integrand, 0., 1., 5);
+  fBGBlastWave_Integrand_num->SetParameters(mt_num, pt, beta_max, temp, n);
+  Double_t integral_num = fBGBlastWave_Integrand_num->Integral(0., 1.);
+  
+  if (!fBGBlastWave_Integrand_den)
+    fBGBlastWave_Integrand_den = new TF1("fBGBlastWave_Integrand_den", BGBlastWave_Integrand, 0., 1., 5);
+  fBGBlastWave_Integrand_den->SetParameters(mt_den, pt, beta_max, temp, n);
+  Double_t integral_den = fBGBlastWave_Integrand_den->Integral(0., 1.);
+
+  return (norm_num / norm_den) * (integral_num / integral_den);
+}
+
+Double_t
+BGBlastWave_Func_OneOverPt(const Double_t *x, const Double_t *p)
+{
+  /* 1/pt dN/dpt */
+  
+  Double_t pt = x[0];
+  Double_t mass = p[0];
+  Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
+  Double_t beta_max = p[1];
+  Double_t temp = p[2];
+  Double_t n = p[3];
+  Double_t norm = p[4];
+  
+  if (!fBGBlastWave_Integrand)
+    fBGBlastWave_Integrand = new TF1("fBGBlastWave_Integrand", BGBlastWave_Integrand, 0., 1., 5);
+  fBGBlastWave_Integrand->SetParameters(mt, pt, beta_max, temp, n);
+  Double_t integral = fBGBlastWave_Integrand->Integral(0., 1., (Double_t *)0, 1.e-3);
+
+  return norm * integral;
+}
+
+TF1 *
+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)
+{
+  
+  TF1 *fBGBlastWave = new TF1(name, BGBlastWave_Func, 0., 10., 5);
+  fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm);
+  fBGBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
+  fBGBlastWave->FixParameter(0, mass);
+  fBGBlastWave->SetParLimits(1, 0.01, 0.99);
+  fBGBlastWave->SetParLimits(2, 0.01, 1.);
+  fBGBlastWave->SetParLimits(3, 0.01, 50.);
+  return fBGBlastWave;
+}
+
+TF1 *
+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)
+{
+  
+  TF1 *fBGBlastWave = new TF1(name, BGBlastWaveRatio_Func, 0., 10., 9);
+  fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm, beta_max, temp, n, norm);
+  fBGBlastWave->SetParNames("mass", "beta_max_num", "T_num", "n_num", "norm_num", "beta_max_den", "T_den", "n_den", "norm_den");
+  fBGBlastWave->FixParameter(0, mass);
+  fBGBlastWave->SetParLimits(1, 0.01, 0.99);
+  fBGBlastWave->SetParLimits(2, 0.01, 1.);
+  fBGBlastWave->SetParLimits(3, 0.01, 10.);
+  fBGBlastWave->SetParLimits(5, 0.01, 0.99);
+  fBGBlastWave->SetParLimits(6, 0.01, 1.);
+  fBGBlastWave->SetParLimits(7, 0.01, 10.);
+  return fBGBlastWave;
+}
+
+TF1 *
+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)
+{
+  
+  TF1 *fBGBlastWave = new TF1(name, BGBlastWaveParticleRatio_Func, 0., 10., 7);
+  fBGBlastWave->SetParameters(mass_num, mass_den, beta_max, temp, n, norm_num, norm_den);
+  fBGBlastWave->SetParNames("mass_num", "mass_den", "beta_max", "T", "n", "norm_num", "norm_den");
+  fBGBlastWave->FixParameter(0, mass_num);
+  fBGBlastWave->FixParameter(1, mass_den);
+  fBGBlastWave->SetParLimits(2, 0.01, 0.99);
+  fBGBlastWave->SetParLimits(3, 0.01, 1.);
+  fBGBlastWave->SetParLimits(4, 0.01, 10.);
+  return fBGBlastWave;
+}
+
+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)
+{
+  
+  TF1 *fBGBlastWave = new TF1(name, BGBlastWave_Func_OneOverPt, 0., 10., 5);
+  fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm);
+  fBGBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
+  fBGBlastWave->FixParameter(0, mass);
+  fBGBlastWave->SetParLimits(1, 0.01, 0.99);
+  fBGBlastWave->SetParLimits(2, 0.01, 1.);
+  fBGBlastWave->SetParLimits(3, 0.01, 50.);
+  return fBGBlastWave;
+}
+
+/*****************************************************************/
+/* TSALLIS BLAST-WAVE */
+/*****************************************************************/
+
+static TF1 *fTsallisBlastWave_Integrand_r = NULL;
+Double_t
+TsallisBlastWave_Integrand_r(const Double_t *x, const Double_t *p)
+{
+  /* 
+     x[0] -> r (radius)
+     p[0] -> mT (transverse mass)
+     p[1] -> pT (transverse momentum)
+     p[2] -> beta_max (surface velocity)
+     p[3] -> T (freezout temperature)
+     p[4] -> n (velocity profile)
+     p[5] -> q
+     p[6] -> y (rapidity)
+     p[7] -> phi (azimuthal angle)
+  */
+  
+  Double_t r = x[0];
+  Double_t mt = p[0];
+  Double_t pt = p[1];
+  Double_t beta_max = p[2];
+  Double_t temp_1 = 1. / p[3];
+  Double_t n = p[4];
+  Double_t q = p[5];
+  Double_t y = p[6];
+  Double_t phi = p[7];
+
+  if (q <= 1.) return r;
+
+  Double_t beta = beta_max * TMath::Power(r, n);
+  Double_t rho = TMath::ATanH(beta);
+  
+  Double_t part1 = mt * TMath::CosH(y) * TMath::CosH(rho);
+  Double_t part2 = pt * TMath::SinH(rho) * TMath::Cos(phi);
+  Double_t part3 = part1 - part2;
+  Double_t part4 = 1 + (q - 1.) * temp_1 * part3;
+  Double_t expo = -1. / (q - 1.);
+  //  printf("part1=%f, part2=%f, part3=%f, part4=%f, expo=%f\n", part1, part2, part3, part4, expo);
+  Double_t part5 = TMath::Power(part4, expo);
+
+  return r * part5;
+}
+
+static TF1 *fTsallisBlastWave_Integrand_phi = NULL;
+Double_t
+TsallisBlastWave_Integrand_phi(const Double_t *x, const Double_t *p)
+{
+  /* 
+     x[0] -> phi (azimuthal angle)
+  */
+  
+  Double_t phi = x[0];
+  fTsallisBlastWave_Integrand_r->SetParameter(7, phi);
+  Double_t integral = fTsallisBlastWave_Integrand_r->Integral(0., 1.);
+  return integral;
+}
+
+static TF1 *fTsallisBlastWave_Integrand_y = NULL;
+Double_t
+TsallisBlastWave_Integrand_y(const Double_t *x, const Double_t *p)
+{
+  /* 
+     x[0] -> y (rapidity)
+  */
+
+  Double_t y = x[0];
+  fTsallisBlastWave_Integrand_r->SetParameter(6, y);
+  Double_t integral = fTsallisBlastWave_Integrand_phi->Integral(-TMath::Pi(), TMath::Pi());
+  return TMath::CosH(y) * integral;
+}
+
+Double_t
+TsallisBlastWave_Func(const Double_t *x, const Double_t *p)
+{
+  /* dN/dpt */
+  
+  Double_t pt = x[0];
+  Double_t mass = p[0];
+  Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
+  Double_t beta_max = p[1];
+  Double_t temp = p[2];
+  Double_t n = p[3];
+  Double_t q = p[4];
+  Double_t norm = p[5];
+
+  if (!fTsallisBlastWave_Integrand_r)
+    fTsallisBlastWave_Integrand_r = new TF1("fTsallisBlastWave_Integrand_r", TsallisBlastWave_Integrand_r, 0., 1., 8);
+  if (!fTsallisBlastWave_Integrand_phi)
+    fTsallisBlastWave_Integrand_phi = new TF1("fTsallisBlastWave_Integrand_phi", TsallisBlastWave_Integrand_phi, -TMath::Pi(), TMath::Pi(), 0);
+  if (!fTsallisBlastWave_Integrand_y)
+    fTsallisBlastWave_Integrand_y = new TF1("fTsallisBlastWave_Integrand_y", TsallisBlastWave_Integrand_y, -0.5, 0.5, 0);
+
+  fTsallisBlastWave_Integrand_r->SetParameters(mt, pt, beta_max, temp, n, q, 0., 0.);
+  Double_t integral = fTsallisBlastWave_Integrand_y->Integral(-0.5, 0.5);
+  return norm * pt * integral;
+}
+
+TF1 *
+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)
+{
+  
+  TF1 *fTsallisBlastWave = new TF1(name, TsallisBlastWave_Func, 0., 10., 6);
+  fTsallisBlastWave->SetParameters(mass, beta_max, temp, n, q, norm);
+  fTsallisBlastWave->SetParNames("mass", "beta_max", "T", "n", "q", "norm");
+  fTsallisBlastWave->FixParameter(0, mass);
+  fTsallisBlastWave->SetParLimits(1, 0.01, 0.99);
+  fTsallisBlastWave->SetParLimits(2, 0.01, 1.);
+  fTsallisBlastWave->SetParLimits(3, 0.1, 10.);
+  fTsallisBlastWave->SetParLimits(4, 1., 10.);
+  return fTsallisBlastWave;
+}
+
+/*****************************************************************/
+/*****************************************************************/
+/*****************************************************************/
+
+
+TF1 *
+BGBlastWave_SingleFit(TH1 *h, Double_t mass, Option_t *opt = "")
+{
+
+  TF1 *f = BGBlastWave(Form("fBGBW_%s", h->GetName()), mass);
+  h->Fit(f);
+  h->Fit(f);
+  h->Fit(f, opt);
+  return f;
+  
+}
+
+Int_t nBW;
+TF1 *fBGBW[1000];
+TF1 *fBGBWratio[1000];
+TGraphErrors *gBW[1000];
+
+TObjArray *
+BGBlastWave_GlobalFit(TObjArray *data, Double_t *mass, Double_t profile = 0.5, Bool_t computeCont = kFALSE, Bool_t fixProfile = kFALSE)
+{
+
+  /* get data */
+  Int_t ndf = 0;
+  nBW = data->GetEntries();
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    gBW[idata] = (TGraphErrors *)data->At(idata);
+    gBW[idata]->SetName(Form("gBW%d", idata));
+    ndf += gBW[idata]->GetN();
+  }
+
+  /* init BG blast-wave functions */
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    printf("init BG-BlastWave function #%d: mass = %f\n", idata, mass[idata]);
+    fBGBW[idata] = BGBlastWave(Form("fBGBW%d", idata), mass[idata]);
+  }
+
+  if (computeCont)
+    printf("-> compute contours requested\n");
+
+  /* display data */
+  TCanvas *cBW = new TCanvas("cBW");
+  cBW->Divide(nBW, 1);
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    cBW->cd(idata + 1);
+    gBW[idata]->Draw("ap*");
+  }
+  cBW->Update();
+
+  /* init minuit: nBW normalizations + 3 (beta, T, n) BG-BlastWave params */
+  const Int_t nbwpars = 3;
+  const Int_t nfitpars = nBW + nbwpars;
+  TMinuit *minuit = new TMinuit(nfitpars);
+  minuit->SetFCN(BGBlastWave_FCN);
+  Double_t arglist[10];
+  Int_t ierflg = 0;
+  arglist[0] = 1;
+  minuit->mnexcm("SET ERR", arglist, 1, ierflg);
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    minuit->mnparm(idata, Form("norm%d", idata), 1.e6, 1., 0., 0., ierflg);
+    ndf--;
+  }
+  //  minuit->mnparm(nBW + 0, "<beta>", 0.55, 0.01, 0., 1., ierflg);
+  //  minuit->mnparm(nBW + 1, "T", 0.14, 0.01, 0., 1., ierflg);
+  //  minuit->mnparm(nBW + 2, "n", profile, 0.1, 0., 10., ierflg);
+  
+  minuit->mnparm(nBW + 0, "<beta>", 0.7, 0.01, 0.2, 0.7, ierflg);
+  minuit->mnparm(nBW + 1, "T", 0.07, 0.001, 0.07, 0.2, ierflg);
+  minuit->mnparm(nBW + 2, "n", profile, 0.1, 0.6, 5., ierflg);
+
+  ndf -= 3;
+  if (fixProfile) {
+    minuit->FixParameter(nBW + 2);
+    ndf++;
+  }
+
+  /* set strategy */
+  arglist[0] = 1;
+  minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
+
+  /* start MIGRAD minimization */
+  arglist[0] = 500000;
+  arglist[1] = 1.;
+  minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
+
+  /* set strategy */
+  arglist[0] = 2;
+  minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
+
+  /* start MIGRAD minimization */
+  arglist[0] = 500000;
+  arglist[1] = 1.;
+  minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
+
+  /* start IMPROVE minimization */
+  arglist[0] = 500000;
+  minuit->mnexcm("IMPROVE", arglist, 1, ierflg);
+
+  /* start MINOS */
+  arglist[0] = 500000;
+  arglist[1] = nBW + 1;
+  arglist[2] = nBW + 2;
+  arglist[3] = nBW + 3;
+  minuit->mnexcm("MINOS", arglist, 4, ierflg);
+
+  /* print results */
+  Double_t amin,edm,errdef;
+  Int_t nvpar,nparx,icstat;
+  minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
+  minuit->mnprin(4, amin);
+
+  /* get parameters */
+  Double_t beta, betae, betaeplus, betaeminus, betagcc, temp, tempe, tempeplus, tempeminus, tempgcc, prof, profe, profeplus, profeminus, profgcc;
+  minuit->GetParameter(nBW + 0, beta, betae);
+  minuit->mnerrs(nBW + 0, betaeplus, betaeminus, betae, betagcc);
+  minuit->GetParameter(nBW + 1, temp, tempe);
+  minuit->mnerrs(nBW + 1, tempeplus, tempeminus, tempe, tempgcc);
+  minuit->GetParameter(nBW + 2, prof, profe);
+  minuit->mnerrs(nBW + 2, profeplus, profeminus, profe, profgcc);
+  Double_t beta_max = 0.5 * (2. + prof) * beta;
+  Double_t norm[1000], norme[1000];
+  for (Int_t idata = 0; idata < nBW; idata++)
+    minuit->GetParameter(idata, norm[idata], norme[idata]);
+
+  /* printout */
+  printf("[x] *********************************\n");
+  printf("[x] beta_max = %f\n", beta_max);
+  printf("[x] <beta>   = %f +- %f (e+ = %f, e- = %f)\n", beta, betae, betaeplus, betaeminus);
+  printf("[x] T        = %f +- %f (e+ = %f, e- = %f)\n", temp, tempe, tempeplus, tempeminus);
+  printf("[x] n        = %f +- %f (e+ = %f, e- = %f)\n", prof, profe, profeplus, profeminus);
+  printf("[x] chi2     = %f\n", amin);
+  printf("[x] ndf      = %f\n", ndf);
+
+  /* 1-sigma contour */
+  minuit->SetErrorDef(1);
+  TGraph *gCont1 = NULL;
+  if (computeCont) gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
+  if (gCont1) gCont1->SetName("gCont1");
+
+  /* 2-sigma contour */
+  minuit->SetErrorDef(4);
+  TGraph *gCont2 = NULL;
+  //  if (computeCont) gCont2 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
+  if (gCont2) gCont2->SetName("gCont2");
+
+  /* display fits */
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    cBW->cd(idata + 1);
+    fBGBW[idata]->SetParameter(4, norm[idata]);
+    fBGBW[idata]->SetParameter(1, beta_max);
+    fBGBW[idata]->SetParameter(2, temp);
+    fBGBW[idata]->SetParameter(3, prof);
+    fBGBW[idata]->Draw("same");
+  }
+  cBW->Update();
+
+  /* histo params */
+  TH1D *hBW = new TH1D("hBW", "", 4, 0., 4.);
+  hBW->SetBinContent(1, beta);
+  hBW->SetBinError(1, betae);
+  hBW->SetBinContent(2, temp);
+  hBW->SetBinError(2, tempe);
+  hBW->SetBinContent(3, prof);
+  hBW->SetBinError(3, profe);
+  hBW->SetBinContent(4, amin/ndf);
+
+  /* BW graph */
+  TGraphAsymmErrors *gBetaT = new TGraphAsymmErrors();
+  gBetaT->SetName("gBetaT");
+  gBetaT->SetPoint(0, beta, temp);
+  gBetaT->SetPointEXlow(0, TMath::Abs(betaeminus));
+  gBetaT->SetPointEXhigh(0, TMath::Abs(betaeplus));
+  gBetaT->SetPointEYlow(0, TMath::Abs(tempeminus));
+  gBetaT->SetPointEYhigh(0, TMath::Abs(tempeplus));
+
+  /* prepare output array */
+  TObjArray *outoa = new TObjArray();
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    outoa->Add(gBW[idata]);
+    outoa->Add(fBGBW[idata]);
+  }
+  outoa->Add(cBW);
+  outoa->Add(hBW);
+  outoa->Add(gBetaT);
+  if (gCont1) outoa->Add(gCont1);
+  if (gCont2) outoa->Add(gCont2);
+
+  return outoa;
+
+}
+
+void 
+BGBlastWave_FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+{
+
+  /* beta -> beta_max */
+  Double_t beta = par[nBW+0];
+  Double_t T = par[nBW+1];
+  Double_t n = par[nBW+2];
+  Double_t beta_max = 0.5 * (2. + n) * beta;
+#if 1
+  /* check beta_max */
+  if (beta_max >= 1. || beta_max <= 0.) {
+    f = kMaxInt;
+    return;
+  }
+  /* check T */
+  if (T <= 0.) {
+    f = kMaxInt;
+    return;
+  }
+#endif
+
+  Double_t pt, pte, val, vale, func, pull, chi = 0;
+  /* loop over all the data */
+  for (Int_t iBW = 0; iBW < nBW; iBW++) {
+    /* set BGBW parameters */
+    fBGBW[iBW]->SetParameter(4, par[iBW]);
+    fBGBW[iBW]->SetParameter(1, beta_max);
+    fBGBW[iBW]->SetParameter(2, T);
+    fBGBW[iBW]->SetParameter(3, n);
+    /* loop over all the points */
+    for (Int_t ipt = 0; ipt < gBW[iBW]->GetN(); ipt++) {
+      pt = gBW[iBW]->GetX()[ipt];
+      pte = gBW[iBW]->GetEX()[ipt];
+      val = gBW[iBW]->GetY()[ipt];
+      vale = gBW[iBW]->GetEY()[ipt];
+      func = fBGBW[iBW]->Eval(pt);
+      //      func = fBGBW[iBW]->Integral(pt - pte, pt + pte);
+      pull = (val - func) / vale;
+      chi += pull * pull;
+    }
+  }
+
+  f = chi;
+}
+
+/*****************************************************************/
+
+TObjArray *
+BGBlastWave_GlobalFitRatio(TObjArray *data, Double_t *mass, Double_t profile = .7, Bool_t fixProfile = kFALSE)
+{
+
+  /* get data */
+  Int_t ndf = 0;
+  nBW = data->GetEntries();
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    gBW[idata] = (TGraphErrors *)data->At(idata);
+    gBW[idata]->SetName(Form("gBW%d", idata));
+    ndf += gBW[idata]->GetN();
+  }
+
+  /* init BG blast-wave functions */
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    printf("init BG-BlastWaveRatio function #%d: mass = %f\n", idata, mass[idata]);
+    fBGBWratio[idata] = BGBlastWaveRatio(Form("fBGBWratio%d", idata), mass[idata]);
+  }
+
+  /* display data */
+  TCanvas *cBW = new TCanvas("cBW");
+  cBW->Divide(nBW, 1);
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    cBW->cd(idata + 1);
+    gBW[idata]->Draw("ap*");
+  }
+  cBW->Update();
+
+  /* init minuit: nBW normalizations + 3 (beta, T, n) BG-BlastWave params */
+  const Int_t nbwpars = 3;
+  const Int_t nfitpars = 2 * (nBW + nbwpars);
+  TMinuit *minuit = new TMinuit(nfitpars);
+  minuit->SetFCN(BGBlastWave_FCNRatio);
+  Double_t arglist[10];
+  Int_t ierflg = 0;
+  arglist[0] = 1;
+  minuit->mnexcm("SET ERR", arglist, 1, ierflg);
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    minuit->mnparm(idata, Form("norm%d_num", idata), 1.e6, 1., 0., 0., ierflg);
+    ndf--;
+  }
+  for (Int_t idata = nBW; idata < 2 * nBW; idata++) {
+    minuit->mnparm(idata, Form("norm%d_den", idata), 1.e6, 1., 0., 0., ierflg);
+    minuit->FixParameter(idata);
+    ndf--;
+  }
+  minuit->mnparm(2 * nBW + 0, "<beta>_num", 0.65, 0.01, 0., 1., ierflg);
+  minuit->mnparm(2 * nBW + 1, "T_num", 0.1, 0.01, 0., 1., ierflg);
+  minuit->mnparm(2 * nBW + 2, "n_num", profile, 0.1, 0., 10., ierflg);
+  minuit->mnparm(2 * nBW + 3, "<beta>_den", 0.65, 0.01, 0., 1., ierflg);
+  minuit->mnparm(2 * nBW + 4, "T_den", 0.1, 0.01, 0., 1., ierflg);
+  minuit->mnparm(2 * nBW + 5, "n_den", profile, 0.1, 0., 10., ierflg);
+  ndf -= 3;
+
+  if (fixProfile) {
+    minuit->FixParameter(nBW + 2);
+    minuit->FixParameter(nBW + 5);
+    ndf++;
+  }
+
+  /* set strategy */
+  arglist[0] = 1;
+  minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
+
+  printf("-->> STARTING MIGRAD with %d parameters\n", nfitpars);
+
+  /* start MIGRAD minimization */
+  arglist[0] = 500000;
+  arglist[1] = 1.;
+  minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
+
+  /* set strategy */
+  arglist[0] = 2;
+  minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
+
+  /* start MIGRAD minimization */
+  arglist[0] = 500000;
+  arglist[1] = 1.;
+  minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
+
+  /* start IMPROVE minimization */
+  arglist[0] = 500000;
+  minuit->mnexcm("IMPROVE", arglist, 1, ierflg);
+
+  /* start MINOS */
+  arglist[0] = 500000;
+  arglist[1] = 2 * nBW + 1;
+  arglist[2] = 2 * nBW + 2;
+  arglist[3] = 2 * nBW + 3;
+  arglist[4] = 2 * nBW + 4;
+  arglist[5] = 2 * nBW + 5;
+  arglist[6] = 2 * nBW + 6;
+  minuit->mnexcm("MINOS", arglist, 7, ierflg);
+
+  /* print results */
+  Double_t amin,edm,errdef;
+  Int_t nvpar,nparx,icstat;
+  minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
+  minuit->mnprin(4, amin);
+
+  /* get parameters */
+  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;
+  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;
+  minuit->GetParameter(2*nBW + 0, beta_num, betae_num);
+  minuit->mnerrs(nBW + 0, betaeplus_num, betaeminus_num, betae_num, betagcc_num);
+  minuit->GetParameter(2*nBW + 1, temp_num, tempe_num);
+  minuit->mnerrs(nBW + 1, tempeplus_num, tempeminus_num, tempe_num, tempgcc_num);
+  minuit->GetParameter(2*nBW + 2, prof_num, profe_num);
+  minuit->mnerrs(nBW + 2, profeplus_num, profeminus_num, profe_num, profgcc_num);
+  minuit->GetParameter(2*nBW + 3, beta_den, betae_den);
+  minuit->mnerrs(nBW + 3, betaeplus_den, betaeminus_den, betae_den, betagcc_den);
+  minuit->GetParameter(2*nBW + 4, temp_den, tempe_den);
+  minuit->mnerrs(nBW + 4, tempeplus_den, tempeminus_den, tempe_den, tempgcc_den);
+  minuit->GetParameter(2*nBW + 5, prof_den, profe_den);
+  minuit->mnerrs(nBW + 5, profeplus_den, profeminus_den, profe_den, profgcc_den);
+  Double_t beta_max_num, beta_max_den;
+  beta_max_num = 0.5 * (2. + prof_num) * beta_num;
+  beta_max_den = 0.5 * (2. + prof_den) * beta_den;
+  Double_t norm_num[1000], norme_num[1000];
+  Double_t norm_den[1000], norme_den[1000];
+  for (Int_t idata = 0; idata < nBW; idata++)
+    minuit->GetParameter(idata, norm_num[idata], norme_num[idata]);
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    minuit->GetParameter(nBW + idata, norm_den[idata], norme_den[idata]);
+  }
+
+  /* printout */
+  printf("[x] *********************************\n");
+  printf("[x] beta_max = %f\n", beta_max_num);
+  printf("[x] <beta>   = %f +- %f (e+ = %f, e- = %f)\n", beta_num, betae_num, betaeplus_num, betaeminus_num);
+  printf("[x] T        = %f +- %f (e+ = %f, e- = %f)\n", temp_num, tempe_num, tempeplus_num, tempeminus_num);
+  printf("[x] n        = %f +- %f (e+ = %f, e- = %f)\n", prof_num, profe_num, profeplus_num, profeminus_num);
+  printf("[x] *********************************\n");
+  printf("[x] beta_max = %f\n", beta_max_den);
+  printf("[x] <beta>   = %f +- %f (e+ = %f, e- = %f)\n", beta_den, betae_den, betaeplus_den, betaeminus_den);
+  printf("[x] T        = %f +- %f (e+ = %f, e- = %f)\n", temp_den, tempe_den, tempeplus_den, tempeminus_den);
+  printf("[x] n        = %f +- %f (e+ = %f, e- = %f)\n", prof_den, profe_den, profeplus_den, profeminus_den);
+  printf("[x] *********************************\n");
+  printf("[x] chi2     = %f\n", amin);
+  printf("[x] ndf      = %f\n", ndf);
+
+  /* 1-sigma contour */
+  minuit->SetErrorDef(1);
+  TGraph *gCont1 = NULL;
+  //  gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
+  if (gCont1) gCont1->SetName("gCont1");
+
+  /* 2-sigma contour */
+  minuit->SetErrorDef(4);
+  TGraph *gCont2 = NULL;
+  //  gCont2 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
+  if (gCont2) gCont2->SetName("gCont2");
+
+  /* display fits */
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    cBW->cd(idata + 1);
+    fBGBWratio[idata]->SetParameter(4, norm_num[idata]);
+    fBGBWratio[idata]->SetParameter(1, beta_max_num);
+    fBGBWratio[idata]->SetParameter(2, temp_num);
+    fBGBWratio[idata]->SetParameter(3, prof_num);
+    fBGBWratio[idata]->SetParameter(8, norm_den[idata]);
+    fBGBWratio[idata]->SetParameter(5, beta_max_den);
+    fBGBWratio[idata]->SetParameter(6, temp_den);
+    fBGBWratio[idata]->SetParameter(7, prof_den);
+    fBGBWratio[idata]->Draw("same");
+  }
+  cBW->Update();
+
+  return;
+
+  /* histo params */
+  TH1D *hBW = new TH1D("hBW", "", 3, 0., 3.);
+  hBW->SetBinContent(1, beta);
+  hBW->SetBinError(1, betae);
+  hBW->SetBinContent(2, temp);
+  hBW->SetBinError(2, tempe);
+  hBW->SetBinContent(3, prof);
+  hBW->SetBinError(3, profe);
+
+  /* BW graph */
+  TGraphAsymmErrors *gBetaT = new TGraphAsymmErrors();
+  gBetaT->SetName("gBetaT");
+  gBetaT->SetPoint(0, beta, temp);
+  gBetaT->SetPointEXlow(0, TMath::Abs(betaeminus));
+  gBetaT->SetPointEXhigh(0, TMath::Abs(betaeplus));
+  gBetaT->SetPointEYlow(0, TMath::Abs(tempeminus));
+  gBetaT->SetPointEYhigh(0, TMath::Abs(tempeplus));
+
+  /* prepare output array */
+  TObjArray *outoa = new TObjArray();
+  for (Int_t idata = 0; idata < nBW; idata++) {
+    outoa->Add(gBW[idata]);
+    outoa->Add(fBGBWratio[idata]);
+  }
+  outoa->Add(cBW);
+  outoa->Add(hBW);
+  outoa->Add(gBetaT);
+  if (gCont1) outoa->Add(gCont1);
+  if (gCont2) outoa->Add(gCont2);
+
+  return outoa;
+
+}
+
+void 
+BGBlastWave_FCNRatio(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+{
+
+  /* beta -> beta_max */
+  Double_t beta_num = par[2*nBW+0];
+  Double_t T_num = par[2*nBW+1];
+  Double_t n_num = par[2*nBW+2];
+  Double_t beta_max_num = 0.5 * (2. + n_num) * beta_num;
+  Double_t beta_den = par[2*nBW+3];
+  Double_t T_den = par[2*nBW+4];
+  Double_t n_den = par[2*nBW+5];
+  Double_t beta_max_den = 0.5 * (2. + n_den) * beta_den;
+#if 0
+  /* check beta_max */
+  if (beta_max >= 1. || beta_max <= 0.) {
+    f = kMaxInt;
+    return;
+  }
+  /* check T */
+  if (T <= 0.) {
+    f = kMaxInt;
+    return;
+  }
+#endif
+
+  Double_t pt, pte, val, vale, func, pull, chi = 0;
+  /* loop over all the data */
+  for (Int_t iBW = 0; iBW < nBW; iBW++) {
+    /* set BGBW parameters */
+    fBGBWratio[iBW]->SetParameter(4, par[iBW]);
+    fBGBWratio[iBW]->SetParameter(1, beta_max_num);
+    fBGBWratio[iBW]->SetParameter(2, T_num);
+    fBGBWratio[iBW]->SetParameter(3, n_num);
+    fBGBWratio[iBW]->SetParameter(8, par[nBW+iBW]);
+    fBGBWratio[iBW]->SetParameter(5, beta_max_den);
+    fBGBWratio[iBW]->SetParameter(6, T_den);
+    fBGBWratio[iBW]->SetParameter(7, n_den);
+    /* loop over all the points */
+    for (Int_t ipt = 0; ipt < gBW[iBW]->GetN(); ipt++) {
+      pt = gBW[iBW]->GetX()[ipt];
+      pte = gBW[iBW]->GetEX()[ipt];
+      val = gBW[iBW]->GetY()[ipt];
+      vale = gBW[iBW]->GetEY()[ipt];
+      func = fBGBWratio[iBW]->Eval(pt);
+      //      func = fBGBW[iBW]->Integral(pt - pte, pt + pte);
+      pull = (val - func) / vale;
+      chi += pull * pull;
+    }
+  }
+
+  f = chi;
+}
+
+/*****************************************************************/
+
+Int_t ratioPartNum[9] = {2, 3, 4, 3, 3, 3, 4, 4, 4};
+Int_t ratioChargeNum[9] = {1, 1, 1, 0, 1, 2, 0, 1, 2};
+Int_t ratioPartDen[9] = {2, 3, 4, 2, 2, 2, 2, 2, 2};
+Int_t ratioChargeDen[9] = {0, 0, 0, 0, 1, 2, 0, 1, 2};
+
+Char_t *ratioPartNumName[9] = {"pion", "kaon", "proton", "kaon", "kaon", "kaon", "proton", "proton", "proton"};
+Char_t *ratioChargeNumName[9] = {"minus", "minus", "minus", "plus", "minus", "both", "plus", "minus", "both"};
+Char_t *ratioPartDenName[9] = {"pion", "kaon", "proton", "pion", "pion", "pion", "pion", "pion", "pion"};
+Char_t *ratioChargeDenName[9] = {"plus", "plus", "plus", "plus", "minus", "both", "plus", "minus", "both"};
+
+IntegratedRatioError(const Char_t *spectrafilename, const Char_t *ratiosfilename)
+{
+
+  TFile *fspectra = TFile::Open(spectrafilename);
+  TFile *fratios = TFile::Open(ratiosfilename);
+  
+  Int_t icent = 0;
+  for (Int_t irat = 3; irat < 9; irat++) {
+
+    Int_t ipnum = ratioPartNum[irat];
+    Int_t icnum = ratioChargeNum[irat];
+
+    Int_t ipden = ratioPartDen[irat];
+    Int_t icden = ratioChargeDen[irat];
+
+    if (icnum == 2 || icden == 2) continue;
+
+    printf("%s\n", Form("sys_cent%d_%s_%s_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat], ratioPartDenName[irat], ratioChargeDenName[irat]));
+    TH1 *hrat_sys = (TH1 *)fratios->Get(Form("sys_cent%d_%s_%s_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat], ratioPartDenName[irat], ratioChargeDenName[irat]));
+
+    printf("%d %d %d %d\n", ipnum, icnum, ipden, icden);
+    TH1 *hnum = (TH1 *)fspectra->Get(Form("cent%d_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat]));
+    TH1 *hnum_stat = (TH1 *)fspectra->Get(Form("stat_cent%d_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat]));
+    printf("%s = %p\n", Form("cent%d_%s_%s", icent, ratioPartNumName[irat], ratioChargeNumName[irat]), hnum);
+  
+    TH1 *hden = (TH1 *)fspectra->Get(Form("cent%d_%s_%s", icent, ratioPartDenName[irat], ratioChargeDenName[irat]));
+    TH1 *hden_stat = (TH1 *)fspectra->Get(Form("stat_cent%d_%s_%s", icent, ratioPartDenName[irat], ratioChargeDenName[irat]));
+    printf("%s = %p\n", Form("cent%d_%s_%s", icent, ratioPartDenName[irat], ratioChargeDenName[irat]), hden);
+    
+    TH1 *hnum_plus = (TH1 *)hnum->Clone("hnum_plus");
+    TH1 *hnum_minus = (TH1 *)hnum->Clone("hnum_minus");
+    TH1 *hden_plus = (TH1 *)hden->Clone("hden_plus");
+    TH1 *hden_minus = (TH1 *)hden->Clone("hden_minus");
+    TH1 *hnum_stat_plus = (TH1 *)hnum_stat->Clone("hnum_stat_plus");
+    TH1 *hnum_stat_minus = (TH1 *)hnum_stat->Clone("hnum_stat_minus");
+    TH1 *hden_stat_plus = (TH1 *)hden_stat->Clone("hden_stat_plus");
+    TH1 *hden_stat_minus = (TH1 *)hden_stat->Clone("hden_stat_minus");
+    for (Int_t ibin = 0; ibin < hnum_plus->GetNbinsX(); ibin++) {
+      Double_t val = hrat_sys->GetBinContent(ibin + 1);
+      Double_t vale = hrat_sys->GetBinError(ibin + 1);
+      if (vale <= 0.) continue;
+      Double_t syse = vale / val;
+
+      val = hnum_plus->GetBinContent(ibin + 1);
+      vale = hnum_plus->GetBinError(ibin + 1);
+      val *= (1. + syse);
+      vale *= (1. + syse);
+      hnum_plus->SetBinContent(ibin + 1, val);
+      hnum_plus->SetBinError(ibin + 1, vale);
+
+      val = hnum_minus->GetBinContent(ibin + 1);
+      vale = hnum_minus->GetBinError(ibin + 1);
+      val *= (1. - syse);
+      vale *= (1. - syse);
+      hnum_minus->SetBinContent(ibin + 1, val);
+      hnum_minus->SetBinError(ibin + 1, vale);
+
+      val = hden_plus->GetBinContent(ibin + 1);
+      vale = hden_plus->GetBinError(ibin + 1);
+      val *= (1. + syse);
+      vale *= (1. + syse);
+      hden_plus->SetBinContent(ibin + 1, val);
+      hden_plus->SetBinError(ibin + 1, vale);
+
+      val = hden_minus->GetBinContent(ibin + 1);
+      vale = hden_minus->GetBinError(ibin + 1);
+      val *= (1. - syse);
+      vale *= (1. - syse);
+      hden_minus->SetBinContent(ibin + 1, val);
+      hden_minus->SetBinError(ibin + 1, vale);
+
+
+
+      val = hnum_stat_plus->GetBinContent(ibin + 1);
+      vale = hnum_stat_plus->GetBinError(ibin + 1);
+      val *= (1. + syse);
+      vale *= (1. + syse);
+      hnum_stat_plus->SetBinContent(ibin + 1, val);
+      hnum_stat_plus->SetBinError(ibin + 1, vale);
+
+      val = hnum_stat_minus->GetBinContent(ibin + 1);
+      vale = hnum_stat_minus->GetBinError(ibin + 1);
+      val *= (1. - syse);
+      vale *= (1. - syse);
+      hnum_stat_minus->SetBinContent(ibin + 1, val);
+      hnum_stat_minus->SetBinError(ibin + 1, vale);
+
+      val = hden_stat_plus->GetBinContent(ibin + 1);
+      vale = hden_stat_plus->GetBinError(ibin + 1);
+      val *= (1. + syse);
+      vale *= (1. + syse);
+      hden_stat_plus->SetBinContent(ibin + 1, val);
+      hden_stat_plus->SetBinError(ibin + 1, vale);
+
+      val = hden_stat_minus->GetBinContent(ibin + 1);
+      vale = hden_stat_minus->GetBinError(ibin + 1);
+      val *= (1. - syse);
+      vale *= (1. - syse);
+      hden_stat_minus->SetBinContent(ibin + 1, val);
+      hden_stat_minus->SetBinError(ibin + 1, vale);
+    }
+
+    TCanvas *c = new TCanvas;
+    hnum_stat->SetLineColor(8);
+    hnum_stat->SetMarkerColor(8);
+    hnum_stat->SetMarkerSize(1);
+    hnum_stat->SetMarkerStyle(20);
+    hnum_stat->Draw();
+    
+    hnum_stat_plus->SetLineColor(2);
+    hnum_stat_plus->SetMarkerColor(2);
+    hnum_stat_plus->SetMarkerSize(1);
+    hnum_stat_plus->SetMarkerStyle(20);
+    hnum_stat_plus->Draw("same");
+    
+    hnum_stat_minus->SetLineColor(4);
+    hnum_stat_minus->SetMarkerColor(4);
+    hnum_stat_minus->SetMarkerSize(1);
+    hnum_stat_minus->SetMarkerStyle(20);
+    hnum_stat_minus->Draw("same");
+    
+    c->Update();
+
+    c = new TCanvas();
+    TH1 *hrat = (TH1 *)hnum_stat->Clone("hrat");
+    TH1 *hrat_plus = (TH1 *)hnum_stat_plus->Clone("hrat_plus");
+    TH1 *hrat_minus = (TH1 *)hnum_stat_minus->Clone("hrat_minus");
+
+    hrat->Divide(hden_stat);
+    hrat_plus->Divide(hden_stat);
+    hrat_minus->Divide(hden_stat);
+
+    hrat->Draw();
+    hrat_plus->Draw("same");
+    hrat_minus->Draw("same");
+
+
+    c->Update();
+    
+    getchar();
+    
+    
+    /* integrate as they are */
+    printf("***** NUM:\n");
+    IntegratedProduction(hnum, AliPID::ParticleMass(ratioPartNum[irat]), 0, "0q");
+    printf("***** NUM PLUS:\n");
+    IntegratedProduction(hnum_plus, AliPID::ParticleMass(ratioPartNum[irat]), 0, "0q");
+    printf("***** NUM MINUS:\n");
+    IntegratedProduction(hnum_minus, AliPID::ParticleMass(ratioPartNum[irat]), 0, "0q");
+    
+    printf("***** DEN:\n");
+    IntegratedProduction(hden, AliPID::ParticleMass(ratioPartDen[irat]), 0, "0q");
+    printf("***** DEN PLUS:\n");
+    IntegratedProduction(hden_plus, AliPID::ParticleMass(ratioPartDen[irat]), 0, "0q");
+    printf("***** DEN MINUS:\n");
+    IntegratedProduction(hden_minus, AliPID::ParticleMass(ratioPartDen[irat]), 0, "0q");
+
+  }
+  
+}
+
+enum EIntegratedData_t {
+  kInt_YieldData, kInt_YieldDataErr, kInt_YieldDataLin,
+  kInt_MeanData, kInt_MeanDataErr, kInt_MeanDataLin,
+  kInt_YieldLow, kInt_YieldLowErr,
+  kInt_MeanLow, kInt_MeanLowErr,
+  kInt_YieldHigh, kInt_YieldHighErr,
+  kInt_MeanHigh, kInt_MeanHighErr,
+  kInt_NData
+};
+
+IntegratedProduction_measurement(const Char_t *filename, Option_t *opt = "q0R")
+{
+
+  for (Int_t ipart = 2; ipart < 5; ipart++)
+    for (Int_t icharge = 0; icharge < 2; icharge++) {
+      IntegratedProduction(filename, ipart, icharge, 0, 0., 10., opt);
+      IntegratedProduction(filename, ipart, icharge, 1, 0., 10., opt);
+    }
+
+}
+
+IntegratedProduction_systematics(const Char_t *filename, Option_t *opt = "q0R")
+{
+
+  Double_t min[5] = {0., 0., 0., 0., 0.};
+  Double_t max[5] = {0., 0., 0.5, 1.0, 2.0};
+  for (Int_t ipart = 2; ipart < 5; ipart++)
+    for (Int_t icharge = 0; icharge < 2; icharge++)
+      for (Int_t ifunc = 1; ifunc < 6; ifunc++)
+       IntegratedProduction(filename, ipart, icharge, ifunc, min[ipart], max[ipart], opt);
+
+}
+
+IntegratedProduction_check(const Char_t *filename, Option_t *opt = "q0R")
+{
+
+  for (Int_t ipart = 2; ipart < 5; ipart++)
+    for (Int_t icharge = 0; icharge < 2; icharge++)
+      for (Int_t ifunc = 1; ifunc < 6; ifunc++)
+       IntegratedProduction(filename, ipart, icharge, ifunc, 0., 10., opt);
+
+}
+
+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")
+{
+
+  Char_t *centName[10] = {
+    "0-5%", "5-10%", "10-20%", "20-40%", "40-60%", "60-80%", "80-100%", "minimum bias"
+  };
+  
+  
+  Char_t *chargeName[3] = {"plus", "minus", "sum"};
+  
+  Char_t *partChargeName[5][2] = {
+    "", "",
+    "", "",
+    "#pi^{+}", "#pi^{-}",
+    "K^{+}", "K^{-}",
+    "p", "#bar{p}"
+  };
+  
+  gStyle->SetOptStat(0);
+  gStyle->SetOptFit(1);
+  
+  /* PWGTools */
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libCORRFW.so");
+  gSystem->Load("libPWGTools.so");
+  AliPWGFunc pwgfunc;
+  pwgfunc.SetVarType(AliPWGFunc::kdNdpt);
+
+  Double_t mass = AliPID::ParticleMass(ipart);
+
+  TF1 *f = NULL;
+  switch (ifunc) {
+  case 0:
+    f = BGBlastWave("BlastWave", mass);
+    f->SetLineColor(2);
+    break;
+  case 1:
+    f = LevyTsallis("Levy-Tsallis", mass, 5., mass);
+    f->SetLineColor(8);
+    break;
+  case 2:
+    f = Boltzmann("Boltzmann", mass);
+    f->SetLineColor(kPink+1);
+    break;
+  case 3:
+    f = pwgfunc.GetMTExp(mass, 0.1, 1., "m_{T}-exponential");
+    f->SetLineColor(kViolet+1);
+    break;
+  case 4:
+    f = pwgfunc.GetPTExp(0.1, 1., "p_{T}-exponential");
+    f->SetLineColor(kAzure+1);
+    break;
+  case 5:
+    f = pwgfunc.GetBoseEinstein(mass, 0.1, 1., "Bose-Einstein");
+    f->SetLineColor(kYellow+1);
+    break;
+  default:
+  };
+  f->SetTitle(f->GetName());
+  f->SetLineWidth(2);
+  f->SetFillColor(0);
+  f->SetMarkerColor(f->GetLineColor());
+
+  /* define low-pt fit range according to the number of free parameters */
+  Int_t nfreeparams = f->GetNumberFreeParameters();
+
+  TCanvas *cc = new TCanvas("cCanvasFit", "", 1200, 800);
+  cc->Divide(4, 2);
+
+  //  if (max < 0)
+  //  TFile *fout = TFile::Open(Form("%s.IntegratedProduction.LowPtFit_%s_%s_%s.root", filename, f->GetName(), AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
+  //  else if (min < 0)
+  //  TFile *fout = TFile::Open(Form("%s.IntegratedProduction.HighPtFit_%s_%s_%s.root", filename, f->GetName(), AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
+  //  else
+    TFile *fout = TFile::Open(Form("%s.IntegratedProduction_%s_%s_%s.root", filename, f->GetName(), AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
+  
+  Double_t integrated_data[kInt_NData];
+
+  TFile *filein = TFile::Open(filename);
+  TH1F *hyield_data_stat = new TH1F(Form("hPartYield_data_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_data_sys = new TH1F(Form("hPartYield_data_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_data = new TH1F(Form("hPartYield_data_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_low_stat = new TH1F(Form("hPartYield_low_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_low_sys = new TH1F(Form("hPartYield_low_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_low = new TH1F(Form("hPartYield_low_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_high_stat = new TH1F(Form("hPartYield_high_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_high_sys = new TH1F(Form("hPartYield_high_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_high = new TH1F(Form("hPartYield_high_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_data_stat = new TH1F(Form("hPartMean_data_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_data_sys = new TH1F(Form("hPartMean_data_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_data = new TH1F(Form("hPartMean_data_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_low_stat = new TH1F(Form("hPartMean_low_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_low_sys = new TH1F(Form("hPartMean_low_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_low = new TH1F(Form("hPartMean_low_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_high_stat = new TH1F(Form("hPartMean_high_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_high_sys = new TH1F(Form("hPartMean_high_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_high = new TH1F(Form("hPartMean_high_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_stat = new TH1F(Form("hYield_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield_sys = new TH1F(Form("hYield_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hyield = new TH1F(Form("hYield_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_stat = new TH1F(Form("hMean_stat_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean_sys = new TH1F(Form("hMean_sys_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  TH1F *hmean = new TH1F(Form("hMean_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", 7, 0, 7);
+  for (Int_t icent = 0; icent < 7; icent++) {
+    cc->cd(icent + 1)->SetLogy();
+    printf("*************************\n");
+    printf("cent=%d part=%d charge=%d\n", icent, ipart, icharge);
+    TH1 *h = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+    //    h = Convert_dNdy_dNdeta(h, AliPID::ParticleMass(ipart), 0.5);
+    h->SetTitle(Form("%s (%s);p_{T} (GeV/c);dN/dp_{T}", partChargeName[ipart][icharge], centName[icent]));
+    h->SetMarkerStyle(20);
+    h->SetMarkerSize(1);
+    h->SetLineWidth(1);
+    h->SetFillColor(0);
+    h->SetFillStyle(0);
+    h->SetLineColor(1);
+    h->SetMarkerColor(4);
+    TH1 *hstat = (TH1 *)filein->Get(Form("stat_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+    //    hstat = Convert_dNdy_dNdeta(hstat, AliPID::ParticleMass(ipart), 0.5);
+    hstat->SetMarkerStyle(20);
+    hstat->SetMarkerSize(1);
+    hstat->SetLineWidth(1);
+    hstat->SetLineColor(1);
+    hstat->SetMarkerColor(4);
+    TH1 *hsys = (TH1 *)filein->Get(Form("sys_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+    //    hsys = Convert_dNdy_dNdeta(hsys, AliPID::ParticleMass(ipart), 0.5);
+    hsys->SetMarkerStyle(20);
+    hsys->SetMarkerSize(1);
+    hsys->SetLineWidth(1);
+    hsys->SetLineColor(1);
+    hsys->SetMarkerColor(4);
+
+    if (max < 0.) {
+      /* define fit range according to number of free parameter */
+      Int_t gotpoints = 0;
+      for (Int_t ipt = 0; ipt < h->GetNbinsX(); ipt++) {
+       if (h->GetBinError(ipt + 1) <= 0.) continue;
+       gotpoints++;
+       if (gotpoints < nfreeparams*2) continue;
+       max = h->GetXaxis()->GetBinUpEdge(ipt + 1);
+       break;
+      }
+      printf("fitting up to %f\n", max);
+    }    
+    if (min < 0.) {
+      /* define fit range according to number of free parameter */
+      Int_t gotpoints = 0;
+      for (Int_t ipt = h->GetNbinsX(); ipt > 0; ipt--) {
+       if (h->GetBinError(ipt) <= 0.) continue;
+       gotpoints++;
+       if (gotpoints < nfreeparams+3) continue;
+       min = h->GetXaxis()->GetBinLowEdge(ipt);
+       break;
+      }
+      printf("fitting from %f\n", min);
+    }    
+
+    printf("-> processing stat-only data\n");
+
+    IntegratedProduction(hstat, f, opt, min, max, integrated_data);
+    fout->cd();
+    hstat->Write(Form("IntegratedProduction_stat_cent%d_%s_%s.root", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+
+    hyield_data_stat->SetBinContent(icent + 1, integrated_data[kInt_YieldData]);
+    hyield_data_stat->SetBinError(icent + 1, integrated_data[kInt_YieldDataErr]);
+    hyield_low_stat->SetBinContent(icent + 1, integrated_data[kInt_YieldLow]);
+    hyield_low_stat->SetBinError(icent + 1, integrated_data[kInt_YieldLowErr]);
+
+    hyield_high_stat->SetBinContent(icent + 1, integrated_data[kInt_YieldHigh]);
+    hyield_high_stat->SetBinError(icent + 1, integrated_data[kInt_YieldHighErr]);
+
+    hmean_data_stat->SetBinContent(icent + 1, integrated_data[kInt_MeanData]);
+    hmean_data_stat->SetBinError(icent + 1, integrated_data[kInt_MeanDataErr]);
+
+    hmean_low_stat->SetBinContent(icent + 1, integrated_data[kInt_MeanLow]);
+    hmean_low_stat->SetBinError(icent + 1, integrated_data[kInt_MeanLowErr]);
+    
+    printf("-> processing sys-only data\n");
+
+    IntegratedProduction(hsys, f, opt, min, max, integrated_data);
+    fout->cd();
+    hsys->Write(Form("IntegratedProduction_sys_cent%d_%s_%s.root", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+
+    hyield_data_sys->SetBinContent(icent + 1, integrated_data[kInt_YieldData]);
+    hyield_data_sys->SetBinError(icent + 1, integrated_data[kInt_YieldDataLin]);
+
+    hyield_low_sys->SetBinContent(icent + 1, integrated_data[kInt_YieldLow]);
+    hyield_low_sys->SetBinError(icent + 1, integrated_data[kInt_YieldLowErr]);
+
+    hyield_high_sys->SetBinContent(icent + 1, integrated_data[kInt_YieldHigh]);
+    hyield_high_sys->SetBinError(icent + 1, integrated_data[kInt_YieldHighErr]);
+
+    hmean_data_sys->SetBinContent(icent + 1, integrated_data[kInt_MeanData]);
+    hmean_data_sys->SetBinError(icent + 1, integrated_data[kInt_MeanDataLin]);
+
+    hmean_low_sys->SetBinContent(icent + 1, integrated_data[kInt_MeanLow]);
+    hmean_low_sys->SetBinError(icent + 1, integrated_data[kInt_MeanLowErr]);
+    
+    printf("-> processing stat+sys data\n");
+
+    IntegratedProduction(h, f, opt, min, max, integrated_data);
+    fout->cd();
+    h->Write(Form("IntegratedProduction_cent%d_%s_%s.root", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+    f->SetRange(min, max);
+    f->Write(Form("%s_cent%d_%s_%s.root", f->GetName(), icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+    h->DrawCopy();
+    f->DrawCopy("same");
+    TLegend *l = gPad->BuildLegend(0.15, 0.12, 0.5, 0.3);
+    l->SetBorderSize(0);
+    l->SetFillColor(0);
+    l->SetFillStyle(0);
+
+    hyield_data->SetBinContent(icent + 1, integrated_data[kInt_YieldData]);
+    hyield_data->SetBinError(icent + 1, integrated_data[kInt_YieldDataLin]);
+
+    hyield_low->SetBinContent(icent + 1, integrated_data[kInt_YieldLow]);
+    hyield_low->SetBinError(icent + 1, integrated_data[kInt_YieldLowErr]);
+
+    hyield_high->SetBinContent(icent + 1, integrated_data[kInt_YieldHigh]);
+    hyield_high->SetBinError(icent + 1, integrated_data[kInt_YieldHighErr]);
+
+    hmean_data->SetBinContent(icent + 1, integrated_data[kInt_MeanData]);
+    hmean_data->SetBinError(icent + 1, integrated_data[kInt_MeanDataLin]);
+
+    hmean_low->SetBinContent(icent + 1, integrated_data[kInt_MeanLow]);
+    hmean_low->SetBinError(icent + 1, integrated_data[kInt_MeanLowErr]);
+    
+    hmean_high->SetBinContent(icent + 1, integrated_data[kInt_MeanHigh]);
+    hmean_high->SetBinError(icent + 1, integrated_data[kInt_MeanHighErr]);
+    
+    
+  }
+
+  hyield_stat->Add(hyield_data_stat);
+  hyield_stat->Add(hyield_low_stat);
+  hyield_stat->Add(hyield_high_stat);
+
+  hyield_sys->Add(hyield_data_sys);
+  hyield_sys->Add(hyield_low_sys);
+  hyield_sys->Add(hyield_high_sys);
+
+  hyield->Add(hyield_data);
+  hyield->Add(hyield_low);
+  hyield->Add(hyield_high);
+  
+  
+  hmean_stat->Add(hmean_data_stat);
+  hmean_stat->Add(hmean_low_stat);
+  hmean_stat->Add(hmean_high_stat);
+  //      hmean_stat->Divide(hyield_stat);
+  hmean_sys->Add(hmean_data_sys);
+  hmean_sys->Add(hmean_low_sys);
+  hmean_sys->Add(hmean_high_sys);
+  //      hmean_sys->Divide(hyield_sys);
+  hmean->Add(hmean_data);
+  hmean->Add(hmean_low);
+  hmean->Add(hmean_high);
+  
+  fout->cd();
+
+  hyield->Write();
+  hmean->Write();
+
+  hyield_stat->Write();
+  hmean_stat->Write();
+
+  hyield_sys->Write();
+  hmean_sys->Write();
+
+  hyield_data->Write();
+  hmean_data->Write();
+
+  hyield_data_stat->Write();
+  hmean_data_stat->Write();
+
+  hyield_data_sys->Write();
+  hmean_data_sys->Write();
+
+  hyield_low->Write();
+  hmean_low->Write();
+
+  hyield_high->Write();
+  hmean_high->Write();
+
+  hyield_low_stat->Write();
+  hmean_low_stat->Write();
+
+  hyield_high_stat->Write();
+  hmean_high_stat->Write();
+
+  hmean_low_sys->Write();
+  hyield_low_sys->Write();
+  
+  hmean_high_sys->Write();
+  hyield_high_sys->Write();
+
+  cc->Write();
+  
+  fout->Close();
+}
+
+IntegratedProduction_pp(const Char_t *filename, Char_t *what = "", Int_t ifunc = 1, Option_t *opt = "0qI")
+{
+  TFile *filein = TFile::Open(filename);
+  for (Int_t ipart = 0; ipart < 6; ipart++) {
+    printf("*************************\n");
+    printf("part=%d\n", ipart);
+    TH1 *h = (TH1 *)filein->Get(Form("hComb_ITSsa%d_TPC%d_TOF%d_HMPID%d", ipart, ipart, ipart, ipart));
+    IntegratedProduction(h, AliPID::ParticleMass(ipart % 3 + 2), ifunc, opt);
+    if (gPad) gPad->SaveAs(Form("IntegratedProduction_pp_%d.root", ipart));
+  }
+}
+
+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)
+{
+
+  Double_t yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, partyield[3], partyielderr[3], partyielderrcorr[3], partmean[3], partmeanerr[3], partmeanerrcorr[3];
+
+  TVirtualFitter::SetMaxIterations(1000000);
+
+  if (f) {
+    Int_t fres = 1;
+    Int_t fatt = 0;
+    while (fres != 0 && fatt < 10) {
+      fres = h->Fit(f, opt, "", min, max);
+      fres = h->Fit(f, opt, "", min, max);
+      printf("fit res = %d\n", fres);
+      fatt++;
+    }
+  }
+
+  // return;
+
+  GetYieldAndMean(h, f, yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, 0., 10., partyield, partyielderr, partyielderrcorr, partmean, partmeanerr, partmeanerrcorr);
+
+  //  Double_t fint = f ? f->Integral(0.,10.) : 0.;
+  //  Double_t finte = f ? f->IntegralError(0.,10.) : 0.;
+  //  Double_t fmean = f ? f->Mean(0., 10.) : 0.;
+
+  if (verbose) {
+  printf("----\n");
+  printf("dN/dy        = %f +- %f (%f)\n", yield, yielderr, yielderrcorr);
+  printf("<pt>         = %f +- %f (%f)\n", mean, meanerr, meanerrcorr);
+  printf("----\n");
+  printf("dN/dy (data) = %f +- %f (%f)\n", partyield[0], partyielderr[0], partyielderrcorr[0]);
+  printf("dN/dy (low)  = %f +- %f (%f)\n", partyield[1], partyielderr[1], partyielderrcorr[1]);
+  printf("dN/dy (high) = %f +- %f (%f)\n", partyield[2], partyielderr[2], partyielderrcorr[2]);
+  printf("<pt> (data)  = %f +- %f (%f)\n", partmean[0], partmeanerr[0], partmeanerrcorr[0]);
+  printf("<pt> (low)   = %f +- %f (%f)\n", partmean[1], partmeanerr[1], partmeanerrcorr[1]);
+  printf("<pt> (high)  = %f +- %f (%f)\n", partmean[2], partmeanerr[2], partmeanerrcorr[2]);
+  printf("----\n");
+  //  printf("dN/dy (func) = %f +- %f\n", fint, finte);
+  //  printf("<pT> (func)  = %f +- %f\n", fmean, 0.);
+  }
+
+  if (!integrated_data) return;
+
+  integrated_data[kInt_YieldData] = partyield[0];
+  integrated_data[kInt_YieldDataErr] = partyielderr[0];
+  integrated_data[kInt_YieldDataLin] = partyielderrcorr[0];
+  integrated_data[kInt_MeanData] = partmean[0];
+  integrated_data[kInt_MeanDataErr] = partmeanerr[0];
+  integrated_data[kInt_MeanDataLin] = partmeanerrcorr[0];
+  integrated_data[kInt_YieldLow] = partyield[1];
+  integrated_data[kInt_YieldLowErr] = partyielderr[1];
+  integrated_data[kInt_MeanLow] = partmean[1];
+  integrated_data[kInt_MeanLowErr] = partmeanerr[1];
+  integrated_data[kInt_YieldHigh] = partyield[2];
+  integrated_data[kInt_YieldHighErr] = partyielderr[2];
+  integrated_data[kInt_MeanHigh] = partmean[2];
+  integrated_data[kInt_MeanHighErr] = partmeanerr[2];
+
+  //  TH1 *hr = (TH1 *)h->Clone("hr");
+  //  hr->Divide(f);
+  //  new TCanvas;
+  //  hr->Draw();
+
+  //  TProfile *p = new TProfile("p", "", 100, 0., 10.);
+  //  gROOT->LoadMacro("HistoUtils.C");
+  //  HistoUtils_Function2Profile(f, p);
+  //  p->Draw();
+}
+
+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)
+{
+
+  /* find lowest edge in histo */
+  Int_t binlo;
+  Double_t lo;
+  for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
+    if (h->GetBinContent(ibin) != 0.) {
+      binlo = ibin;
+      lo = h->GetBinLowEdge(ibin);
+      break;
+    }
+  }
+  
+  /* find highest edge in histo */
+  Int_t binhi;
+  Double_t hi;
+  for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
+    if (h->GetBinContent(ibin) != 0.) {
+      binhi = ibin + 1;
+      hi = h->GetBinLowEdge(ibin + 1);
+      break;
+    }
+  }
+  
+  /* integrate the data */
+  Double_t cont, err, width, cent, integral_data = 0., integralerr_data = 0., integralerrcorr_data = 0., meanintegral_data = 0., meanintegralerr_data = 0., meanintegralerrcorr_data = 0.;
+  for (Int_t ibin = binlo; ibin < binhi; ibin++) {
+    cent = h->GetBinCenter(ibin);
+    width = h->GetBinWidth(ibin);
+    cont = h->GetBinContent(ibin);
+    err = h->GetBinError(ibin);
+    /* check we didn't get an empty bin in between */
+    if (cont != 0. && err != 0.) {
+      /* all right, use data */
+      integral_data += cont * width;
+      integralerr_data += err * err * width * width;
+      integralerrcorr_data += err * width;
+      meanintegral_data += cont * width * cent;
+      meanintegralerr_data += err * err * width * width * cent * cent;
+      meanintegralerrcorr_data += err * width * cent;
+    }
+    else {
+      /* missing data-point, complain and use function */
+      printf("WARNING: missing data-point at %f\n", cent);
+      printf("         using function as a patch\n");
+      integral_data += f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
+      integralerr_data += f->IntegralError(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1), 0, 0, 1.e-6);
+      meanintegral_data += f->Mean(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1)) * f->Integral(h->GetBinLowEdge(ibin), h->GetBinLowEdge(ibin+1));
+      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);
+    }
+  }
+  integralerr_data = TMath::Sqrt(integralerr_data);
+  meanintegralerr_data = TMath::Sqrt(meanintegralerr_data);
+  
+  /* integrate below the data */
+  if (!f) {
+    min = lo;
+    max = hi;
+    printf("WARNING: no function provided! Use only data with no extrapolation\n");
+  }
+  Double_t integral_lo = min < lo ? f->Integral(min, lo, (Double_t *)0, 1.e-6) : 0.;
+  Double_t integralerr_lo = min < lo ? f->IntegralError(min, lo, 0, 0, 1.e-6) : 0.;
+  Double_t meanintegral_lo = min < lo ? f->Mean(min, lo, (Double_t *)0, 1.e-6) * integral_lo : 0.;
+  Double_t meanintegralerr_lo = min < lo ? f->Mean(min, lo, (Double_t *)0, 1.e-6) * integralerr_lo : 0.;
+  
+  /* integrate above the data */
+  Double_t integral_hi = max > hi ? f->Integral(hi, max, (Double_t *)0, 1.e-6) : 0.;
+  Double_t integralerr_hi = max > hi ? f->IntegralError(hi, max, 0, 0, 1.e-6) : 0.;
+  Double_t meanintegral_hi = max > hi ? f->Mean(hi, max, (Double_t *)0, 1.e-6) * integral_hi : 0.;
+  Double_t meanintegralerr_hi = max > hi ? f->Mean(hi, max, (Double_t *)0, 1.e-6) * integralerr_hi : 0.;
+
+  /* compute integrated yield */
+  yield = integral_data + integral_lo + integral_hi;
+  yielderr = TMath::Sqrt(integralerr_data * integralerr_data + 
+                        integralerr_lo * integralerr_lo + 
+                        integralerr_hi * integralerr_hi);
+  yielderrcorr = TMath::Sqrt(integralerrcorr_data * integralerrcorr_data + 
+                            integralerr_lo * integralerr_lo + 
+                            integralerr_hi * integralerr_hi);
+  
+  /* compute integrated mean */
+  mean = (meanintegral_data + meanintegral_lo + meanintegral_hi) / yield;
+  meanerr = TMath::Sqrt(meanintegralerr_data * meanintegralerr_data + 
+                       meanintegralerr_lo * meanintegralerr_lo + 
+                       meanintegralerr_hi * meanintegralerr_hi) / yield;
+  meanerrcorr = TMath::Sqrt(meanintegralerrcorr_data * meanintegralerrcorr_data + 
+                           meanintegralerr_lo * meanintegralerr_lo + 
+                           meanintegralerr_hi * meanintegralerr_hi) / yield;
+
+  /* set partial yields */
+  partyield[0] = integral_data;
+  partyielderr[0] = integralerr_data;
+  partyielderrcorr[0] = integralerrcorr_data;
+  partyield[1] = integral_lo;
+  partyielderr[1] = integralerr_lo;
+  partyielderrcorr[1] = integralerr_lo;
+  partyield[2] = integral_hi;
+  partyielderr[2] = integralerr_hi;
+  partyielderrcorr[2] = integralerr_hi;
+
+  /* set partial means */
+  partmean[0] = meanintegral_data;
+  partmeanerr[0] = meanintegralerr_data;
+  partmeanerrcorr[0] = meanintegralerrcorr_data;
+  partmean[1] = meanintegral_lo;
+  partmeanerr[1] = meanintegralerr_lo;
+  partmeanerrcorr[1] = meanintegralerr_lo;
+  partmean[2] = meanintegral_hi;
+  partmeanerr[2] = meanintegralerr_hi;
+  partmeanerrcorr[2] = meanintegralerr_hi;
+  
+}
+
+/*****************************************************************/
+
+Double_t 
+y2eta(Double_t pt, Double_t mass, Double_t y){
+  Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
+  return TMath::ASinH(mt / pt * TMath::SinH(y));
+}
+Double_t 
+eta2y(Double_t pt, Double_t mass, Double_t eta){
+  Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
+  return TMath::ASinH(pt / mt * TMath::SinH(eta));
+}
+
+TH1 *
+Convert_dNdy_1over2pipt_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
+{
+
+  TH1 *hout = hin->Clone("hout");
+  hout->Reset();
+  Double_t pt, mt, conv, val, vale;
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    pt = hin->GetBinCenter(ibin + 1);
+    conv = eta2y(pt, mass, eta) / eta;
+    val = hin->GetBinContent(ibin + 1);
+    vale = hin->GetBinError(ibin + 1);
+    val /= (2. * TMath::Pi() * pt);
+    vale /= (2. * TMath::Pi() * pt);
+    val *= conv;
+    vale *= conv;
+    hout->SetBinContent(ibin + 1, val);
+    hout->SetBinError(ibin + 1, vale);
+  }
+  return hout;
+}
+
+TH1 *
+Convert_dNdy_1over2pipt_dNdy(TH1 *hin)
+{
+
+  TH1 *hout = hin->Clone("hout");
+  hout->Reset();
+  Double_t pt, mt, conv, val, vale;
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    pt = hin->GetBinCenter(ibin + 1);
+    val = hin->GetBinContent(ibin + 1);
+    vale = hin->GetBinError(ibin + 1);
+    val /= (2. * TMath::Pi() * pt);
+    vale /= (2. * TMath::Pi() * pt);
+    hout->SetBinContent(ibin + 1, val);
+    hout->SetBinError(ibin + 1, vale);
+  }
+  return hout;
+}
+
+TH1 *
+Convert_dNdy_1overpt_dNdy(TH1 *hin)
+{
+
+  TH1 *hout = hin->Clone("hout");
+  hout->Reset();
+  Double_t pt, mt, conv, val, vale;
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    pt = hin->GetBinCenter(ibin + 1);
+    val = hin->GetBinContent(ibin + 1);
+    vale = hin->GetBinError(ibin + 1);
+    val /= pt;
+    vale /= pt;
+    hout->SetBinContent(ibin + 1, val);
+    hout->SetBinError(ibin + 1, vale);
+  }
+  return hout;
+}
+
+TH1 *
+Convert_dNdy_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
+{
+
+  TH1 *hout = hin->Clone("hout");
+  hout->Reset();
+  Double_t pt, mt, conv, val, vale;
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    pt = hin->GetBinCenter(ibin + 1);
+    conv = eta2y(pt, mass, eta) / eta;
+    val = hin->GetBinContent(ibin + 1);
+    vale = hin->GetBinError(ibin + 1);
+    val *= conv;
+    vale *= conv;
+    hout->SetBinContent(ibin + 1, val);
+    hout->SetBinError(ibin + 1, vale);
+  }
+  return hout;
+}
+
+TGraph *
+Convert_dNdy_dNdeta(TGraph *hin, Double_t mass, Double_t eta = 0.8)
+{
+
+  TGraph *hout = hin->Clone("hout");
+  //  hout->Reset();
+  Double_t pt, mt, conv, val, valelo, valehi;
+  for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
+    pt = hin->GetX()[ibin];
+    conv = eta2y(pt, mass, eta) / eta;
+    val = hin->GetY()[ibin];
+    valelo = hin->GetEYlow()[ibin];
+    valehi = hin->GetEYhigh()[ibin];
+    val *= conv;
+    valelo *= conv;
+    valehi *= conv;
+    hout->GetX()[ibin] = pt;
+    hout->GetY()[ibin] = val;
+    hout->GetEYlow()[ibin] = valelo;
+    hout->GetEYhigh()[ibin] = valehi;
+  }
+  return hout;
+}
+
+TH1 *
+SummedId_1over2pipt_dNdeta(const Char_t *filename, Int_t icent, Float_t etarange, Float_t scalePi = 1.)
+{
+
+  const Char_t *chargeName[2] = {
+    "plus", "minus"
+  };
+
+  TFile *filein = TFile::Open(filename);
+  TH1 *hy[AliPID::kSPECIES][2];
+  TH1 *heta[AliPID::kSPECIES][2];
+  for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
+    for (Int_t icharge = 0; icharge < 2; icharge++) {
+      hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+      if (!hy[ipart][icharge]) {
+       printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
+       return NULL;
+      }
+      heta[ipart][icharge] = Convert_dNdy_1over2pipt_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart), etarange);
+      if (ipart == 2) heta[ipart][icharge]->Scale(scalePi);
+    }
+
+  /* sum */
+  TH1D *hsum = heta[2][0]->Clone("hsum");
+  hsum->Reset();
+  for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
+    for (Int_t icharge = 0; icharge < 2; icharge++) {
+      hsum->Add(heta[ipart][icharge]);
+    }
+  for (Int_t ipt = 0; ipt < hsum->GetNbinsX(); ipt++) {
+    Double_t err = 0.;
+    for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
+      for (Int_t icharge = 0; icharge < 2; icharge++) {
+       err += heta[ipart][icharge]->GetBinError(ipt + 1);
+      }
+    hsum->SetBinError(ipt + 1, err);
+  }
+  return hsum;
+}
+
+TH1 *
+SummedId_dNdeta(const Char_t *filename, Int_t icent)
+{
+
+  const Char_t *chargeName[2] = {
+    "plus", "minus"
+  };
+
+  TFile *filein = TFile::Open(filename);
+  TH1 *hy[AliPID::kSPECIES][2];
+  TH1 *heta[AliPID::kSPECIES][2];
+  for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
+    for (Int_t icharge = 0; icharge < 2; icharge++) {
+      hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+      if (!hy[ipart][icharge]) {
+       printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
+       return NULL;
+      }
+      heta[ipart][icharge] = Convert_dNdy_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
+    }
+
+  /* sum */
+  TH1D *hsum = heta[2][0]->Clone("hsum");
+  hsum->Reset();
+  for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
+    for (Int_t icharge = 0; icharge < 2; icharge++)
+      hsum->Add(heta[ipart][icharge]);
+
+  return hsum;
+}
+
+/*****************************************************************/
+
+TH1 *
+ReturnExtremeHighHisto(TH1 *hin, TH1 *herr = NULL)
+{
+  TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehigh", hin->GetName()));
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    Double_t val = hin->GetBinContent(ibin + 1);
+    Double_t err = hin->GetBinError(ibin + 1);
+    if (herr) err = herr->GetBinError(ibin + 1);
+    hout->SetBinContent(ibin + 1, val + err);
+  }
+  return hout;
+}
+
+TH1 *
+ReturnExtremeLowHisto(TH1 *hin, TH1 *herr = NULL)
+{
+  TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremelow", hin->GetName()));
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    Double_t val = hin->GetBinContent(ibin + 1);
+    Double_t err = hin->GetBinError(ibin + 1);
+    if (herr) err = herr->GetBinError(ibin + 1);
+    hout->SetBinContent(ibin + 1, val - err);
+  }
+  return hout;
+}
+
+TH1 *
+ReturnExtremeSoftHisto(TH1 *hin, TH1 *herr = NULL)
+{
+  return ReturnExtremeHisto(hin, herr, -1.);
+}
+
+TH1 *
+ReturnExtremeHardHisto(TH1 *hin, TH1 *herr = NULL)
+{
+  return ReturnExtremeHisto(hin, herr, 1.);
+}
+
+TH1 *
+ReturnExtremeHisto(TH1 *hin, TH1 *herr = NULL, Float_t sign = 1.)
+{
+  Double_t ptlow, pthigh;
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    ptlow = hin->GetBinLowEdge(ibin + 1);
+    break;
+  }
+  for (Int_t ibin = hin->GetNbinsX(); ibin >= 0; ibin--) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    pthigh = hin->GetBinLowEdge(ibin + 2);
+    break;
+  }
+
+  Double_t mean = hin->GetMean();
+  Double_t maxdiff = 0.;
+  TH1 *hmax = NULL;
+  for (Int_t inode = 0; inode < hin->GetNbinsX(); inode++) {
+
+    Double_t ptnode = hin->GetBinCenter(inode + 1);
+    TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehard", hin->GetName()));
+    
+    for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+      if (hin->GetBinError(ibin + 1) <= 0.) continue;
+      Double_t val = hin->GetBinContent(ibin + 1);
+      Double_t err = hin->GetBinError(ibin + 1);
+      if (herr) err = herr->GetBinError(ibin + 1);
+      Double_t cen = hin->GetBinCenter(ibin + 1);
+      //    err *= -1. - (cen - ptlow) * (-1. - 1.) / (pthigh - ptlow);
+      if (cen < ptnode)
+        err *= -1. + (cen - ptlow) / (ptnode - ptlow);
+      else
+        err *= (cen - ptnode) / (pthigh - ptnode);
+
+      hout->SetBinContent(ibin + 1, val + sign * err);
+    }
+
+    Double_t diff = TMath::Abs(mean - hout->GetMean());
+    if (diff > maxdiff) {
+      //      printf("found max at %f\n", ptnode);
+      if (hmax) delete hmax;
+      hmax = (TH1 *)hout->Clone("hmax");
+      maxdiff = diff;
+    }
+    delete hout;
+  }
+  return hmax;
+}
+
+TGraphErrors *
+ReturnExtremeSoftGraph(TGraphErrors *hin, TGraphErrors *herr = NULL)
+{
+  TGraphErrors *hout = (TGraphErrors *)hin->Clone(Form("%s_extremesoft", hin->GetName()));
+  Double_t ptlow, pthigh;
+  Double_t ptlow = hin->GetX()[0];
+  Double_t pthigh = hin->GetX()[hin->GetN() - 1];
+
+  for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
+    Double_t val = hin->GetY()[ibin];
+    Double_t err = hin->GetEY()[ibin];
+    if (herr) err = herr->GetEY()[ibin];
+    Double_t cen = hin->GetX()[ibin];
+    err *= 1. + (cen - ptlow) * (-1. - 1.) / (pthigh - ptlow);
+    hout->SetPoint(ibin, cen, val + err);
+  }
+  return hout;
+}
+
+
+TGraphErrors *
+ReturnExtremeHardGraph(TGraphErrors *hin, TGraphErrors *herr = NULL)
+{
+  TGraphErrors *hout = (TGraphErrors *)hin->Clone(Form("%s_extremehard", hin->GetName()));
+  Double_t ptlow, pthigh;
+  Double_t ptlow = hin->GetX()[0];
+  Double_t pthigh = hin->GetX()[hin->GetN() - 1];
+
+  for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
+    Double_t val = hin->GetY()[ibin];
+    Double_t err = hin->GetEY()[ibin];
+    if (herr) err = herr->GetEY()[ibin];
+    Double_t cen = hin->GetX()[ibin];
+    err *= -1. - (cen - ptlow) * (-1. - 1.) / (pthigh - ptlow);
+    hout->SetPoint(ibin, cen, val + err);
+  }
+  return hout;
+}
+
diff --git a/PWGLF/SPECTRA/UTILS/YieldMean.C b/PWGLF/SPECTRA/UTILS/YieldMean.C
new file mode 100644 (file)
index 0000000..a2312b4
--- /dev/null
@@ -0,0 +1,419 @@
+/* definition of the fields in the histogram returned */
+enum EValue_t {
+  kYield = 1,
+  kYieldStat,
+  kYieldSysHi,
+  kYieldSysLo,
+  kMean,
+  kMeanStat,
+  kMeanSysHi,
+  kMeanSysLo
+};
+
+TH1 *
+YieldMean(TH1 *hstat, TH1 *hsys, TF1 *f = NULL, Double_t min = 0., Double_t max = 10., Double_t loprecision = 0.01, Double_t hiprecision = 0.1, Option_t *opt = "0q")
+{
+  /* set many iterations when fitting the data so we don't
+     stop minimization with MAX_CALLS */
+  TVirtualFitter::SetMaxIterations(1000000);
+
+  /* create output histo */
+  Double_t integral, mean;
+  TH1 *hout = new TH1D("hout", "", 8, 0, 8);
+  TH1 *hlo, *hhi;
+  
+  /* create histo with stat+sys errors */
+  TH1 *htot = (TH1 *)hstat->Clone("htot");
+  for (Int_t ibin = 0; ibin < htot->GetNbinsX(); ibin++) {
+    htot->SetBinError(ibin + 1, TMath::Sqrt(hsys->GetBinError(ibin + 1) * hsys->GetBinError(ibin + 1) + hstat->GetBinError(ibin + 1) * hstat->GetBinError(ibin + 1)));
+  }
+
+  /*
+   *   measure the central value 
+   */
+  
+  do Int_t fitres = htot->Fit(f, opt);
+  while (fitres != 0);
+  hlo = YieldMean_LowExtrapolationHisto(htot, f, min, loprecision);
+  hhi = YieldMean_HighExtrapolationHisto(htot, f, max, hiprecision);
+  YieldMean_IntegralMean(htot, hlo, hhi, integral, mean);
+  hout->SetBinContent(kYield, integral);
+  hout->SetBinContent(kMean, mean);
+
+  /*
+   * STATISTICS
+   */
+  
+  TCanvas *cCanvasStat = new TCanvas("cCanvasStat");
+  cCanvasStat->Divide(2, 1);
+  
+  /*
+   * measure statistical error
+   */
+
+  /* fit with stat error */
+  do Int_t fitres = hstat->Fit(f, opt);
+  while (fitres != 0);
+  hlo = YieldMean_LowExtrapolationHisto(hstat, f, min, loprecision);
+  hhi = YieldMean_HighExtrapolationHisto(hstat, f, max, hiprecision);
+  
+  /* random generation with integration (coarse) */
+  TH1 *hIntegral_tmp = new TH1F("hIntegral_tmp", "", 1000, 0.75 * integral, 1.25 * integral);
+  TH1 *hMean_tmp = new TH1F("hMean_tmp", "", 1000, 0.75 * mean, 1.25 * mean);
+  for (Int_t irnd = 0; irnd < 100; irnd++) {
+    /* get random histogram */
+    TH1 *hrnd = YieldMean_ReturnRandom(hstat);
+    /* fit */
+    TH1 *hrndlo = YieldMean_ReturnCoherentRandom(hlo);
+    TH1 *hrndhi = YieldMean_ReturnCoherentRandom(hhi);
+    /* integrate */
+    YieldMean_IntegralMean(hrnd, hrndlo, hrndhi, integral, mean);
+    hIntegral_tmp->Fill(integral);
+    hMean_tmp->Fill(mean);
+    delete hrnd;
+    delete hrndlo;
+    delete hrndhi;
+  }
+  /* random generation with integration (fine) */
+  TH1 *hIntegral = new TH1F("hIntegral", "", 100, 
+                            hIntegral_tmp->GetMean() - 10. * hIntegral_tmp->GetRMS(),
+                            hIntegral_tmp->GetMean() + 10. * hIntegral_tmp->GetRMS());
+  TH1 *hMean = new TH1F("hMean", "", 100,
+                        hMean_tmp->GetMean() - 10. * hMean_tmp->GetRMS(),
+                        hMean_tmp->GetMean() + 10. * hMean_tmp->GetRMS());
+  for (Int_t irnd = 0; irnd < 1000; irnd++) {
+    /* get random histogram */
+    TH1 *hrnd = YieldMean_ReturnRandom(hstat);
+    /* fit */
+    TH1 *hrndlo = YieldMean_ReturnCoherentRandom(hlo);
+    TH1 *hrndhi = YieldMean_ReturnCoherentRandom(hhi);
+    /* integrate */
+    YieldMean_IntegralMean(hrnd, hrndlo, hrndhi, integral, mean);
+    hIntegral->Fill(integral);
+    hMean->Fill(mean);
+    delete hrnd;
+    delete hrndlo;
+    delete hrndhi;
+  }
+  TF1 *gaus = (TF1 *)gROOT->GetFunction("gaus");
+  
+  cCanvasStat->cd(1);
+  hIntegral->Fit(gaus, "q");
+  integral = hout->GetBinContent(kYield) * gaus->GetParameter(2) / gaus->GetParameter(1);
+  hout->SetBinContent(kYieldStat, integral);
+  
+  cCanvasStat->cd(2);
+  hMean->Fit(gaus, "q");
+  mean = hout->GetBinContent(kMean) * gaus->GetParameter(2) / gaus->GetParameter(1);
+  hout->SetBinContent(kMeanStat, mean);
+  
+  /*
+   * SYSTEMATICS
+   */
+
+  TCanvas *cCanvasSys = new TCanvas("cCanvasYieldSys");
+  cCanvasSys->Divide(2, 1);
+  cCanvasSys->cd(1)->DrawFrame(min, 1.e-3, max, 1.e3);
+  hsys->SetMarkerStyle(20);
+  hsys->SetMarkerColor(1);
+  hsys->SetMarkerSize(1);
+  hsys->Draw("same");
+  cCanvasSys->cd(2)->DrawFrame(min, 1.e-3, max, 1.e3);
+  hsys->Draw("same");
+  
+  /*
+   * systematic error high
+   */
+
+  TH1 *hhigh = YieldMean_ReturnExtremeHighHisto(hsys);
+  do Int_t fitres = hhigh->Fit(f, opt);
+  while (fitres != 0);
+  hlo = YieldMean_LowExtrapolationHisto(hhigh, f, min, loprecision);
+  hhi = YieldMean_HighExtrapolationHisto(hhigh, f, max, hiprecision);
+  YieldMean_IntegralMean(hhigh, hlo, hhi, integral, mean);
+  integral = TMath::Abs(integral - hout->GetBinContent(kYield));
+  hout->SetBinContent(kYieldSysHi, integral);
+
+  cCanvasSys->cd(1);
+  f->SetLineColor(2);
+  f->DrawCopy("same");
+  
+  /*
+   * systematic error hard
+   */
+
+  TH1 *hhard = YieldMean_ReturnExtremeHardHisto(hsys);
+  do Int_t fitres = hhard->Fit(f, opt);
+  while (fitres != 0);
+  hlo = YieldMean_LowExtrapolationHisto(hhard, f, min, loprecision);
+  hhi = YieldMean_HighExtrapolationHisto(hhard, f, max, hiprecision);
+  YieldMean_IntegralMean(hhard, hlo, hhi, integral, mean);
+  mean = TMath::Abs(mean - hout->GetBinContent(kMean));
+  hout->SetBinContent(kMeanSysHi, mean);
+
+  cCanvasSys->cd(2);
+  f->SetLineColor(2);
+  f->DrawCopy("same");
+  
+  /*
+   * systematic error low
+   */
+
+  TH1 *hlow = YieldMean_ReturnExtremeLowHisto(hsys);
+  do Int_t fitres = hlow->Fit(f, opt);
+  while (fitres != 0);
+  hlo = YieldMean_LowExtrapolationHisto(hlow, f, min, loprecision);
+  hhi = YieldMean_HighExtrapolationHisto(hlow, f, max, hiprecision);
+  YieldMean_IntegralMean(hlow, hlo, hhi, integral, mean);
+  integral = TMath::Abs(integral - hout->GetBinContent(kYield));
+  hout->SetBinContent(kYieldSysLo, integral);
+
+  cCanvasSys->cd(1);
+  f->SetLineColor(4);
+  f->DrawCopy("same");
+
+  /*
+   * systematic error soft
+   */
+
+  TH1 *hsoft = YieldMean_ReturnExtremeSoftHisto(hsys);
+  do Int_t fitres = hsoft->Fit(f, opt);
+  while (fitres != 0);
+  hlo = YieldMean_LowExtrapolationHisto(hsoft, f, min, loprecision);
+  hhi = YieldMean_HighExtrapolationHisto(hsoft, f, max, hiprecision);
+  YieldMean_IntegralMean(hsoft, hlo, hhi, integral, mean);
+  mean = TMath::Abs(mean - hout->GetBinContent(kMean));
+  hout->SetBinContent(kMeanSysLo, mean);
+
+  cCanvasSys->cd(2);
+  f->SetLineColor(4);
+  f->DrawCopy("same");
+
+  return hout;
+}
+
+TH1 *
+YieldMean_LowExtrapolationHisto(TH1 *h, TF1 *f, Double_t min, Double_t binwidth = 0.01)
+{
+  /* find lowest edge in histo */
+  Int_t binlo;
+  Double_t lo;
+  for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
+    if (h->GetBinContent(ibin) != 0.) {
+      binlo = ibin;
+      lo = h->GetBinLowEdge(ibin);
+      break;
+    }
+  }
+  
+  Int_t nbins = (lo - min) / binwidth;
+  TH1 *hlo = new TH1F("hlo", "", nbins, min, lo);
+  
+  /* integrate function in histogram bins */
+  Double_t cont, err, width;
+  for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) {
+    width = hlo->GetBinWidth(ibin + 1);
+    cont = f->Integral(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
+    err = f->IntegralError(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
+    hlo->SetBinContent(ibin + 1, cont / width);
+    hlo->SetBinError(ibin + 1, err / width);
+  }
+
+  return hlo;
+}
+
+TH1 *
+YieldMean_HighExtrapolationHisto(TH1 *h, TF1 *f, Double_t max, Double_t binwidth = 0.1)
+{
+  /* find highest edge in histo */
+  Int_t binhi;
+  Double_t hi;
+  for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
+    if (h->GetBinContent(ibin) != 0.) {
+      binhi = ibin + 1;
+      hi = h->GetBinLowEdge(ibin + 1);
+      break;
+    }
+  }
+  
+  Int_t nbins = (max - hi) / binwidth;
+  TH1 *hhi = new TH1F("hhi", "", nbins, hi, max);
+  
+  /* integrate function in histogram bins */
+  Double_t cont, err, width;
+  for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) {
+    width = hhi->GetBinWidth(ibin + 1);
+    cont = f->Integral(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
+    err = f->IntegralError(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
+    hhi->SetBinContent(ibin + 1, cont / width);
+    hhi->SetBinError(ibin + 1, err / width);
+  }
+
+  return hhi;
+}
+
+TH1 *
+YieldMean_ReturnRandom(TH1 *hin)
+{
+  TH1 *hout = (TH1 *)hin->Clone("hout");
+  hout->Reset();
+  Double_t cont, err;
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    cont = hin->GetBinContent(ibin + 1);
+    err = hin->GetBinError(ibin + 1);
+    hout->SetBinContent(ibin + 1, gRandom->Gaus(cont, err));
+    hout->SetBinError(ibin + 1, err);
+  }
+  return hout;
+}
+
+TH1 *
+YieldMean_ReturnCoherentRandom(TH1 *hin)
+{
+  TH1 *hout = (TH1 *)hin->Clone("hout");
+  hout->Reset();
+  Double_t cont, err, cohe;
+  cohe = gRandom->Gaus(0., 1.);
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    cont = hin->GetBinContent(ibin + 1);
+    err = hin->GetBinError(ibin + 1);
+    hout->SetBinContent(ibin + 1, cont + cohe * err);
+    hout->SetBinError(ibin + 1, err);
+  }
+  return hout;
+}
+
+TH1 *
+YieldMean_ReturnExtremeHighHisto(TH1 *hin)
+{
+  TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehigh", hin->GetName()));
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    Double_t val = hin->GetBinContent(ibin + 1);
+    Double_t err = hin->GetBinError(ibin + 1);
+    hout->SetBinContent(ibin + 1, val + err);
+  }
+  return hout;
+}
+
+TH1 *
+YieldMean_ReturnExtremeLowHisto(TH1 *hin)
+{
+  TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremelow", hin->GetName()));
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    Double_t val = hin->GetBinContent(ibin + 1);
+    Double_t err = hin->GetBinError(ibin + 1);
+    hout->SetBinContent(ibin + 1, val - err);
+  }
+  return hout;
+}
+
+TH1 *
+YieldMean_ReturnExtremeSoftHisto(TH1 *hin)
+{
+  return YieldMean_ReturnExtremeHisto(hin, -1.);
+}
+
+TH1 *
+YieldMean_ReturnExtremeHardHisto(TH1 *hin)
+{
+  return YieldMean_ReturnExtremeHisto(hin, 1.);
+}
+
+TH1 *
+YieldMean_ReturnExtremeHisto(TH1 *hin, Float_t sign = 1.)
+{
+  Double_t ptlow, pthigh;
+  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    ptlow = hin->GetBinLowEdge(ibin + 1);
+    break;
+  }
+  for (Int_t ibin = hin->GetNbinsX(); ibin >= 0; ibin--) {
+    if (hin->GetBinError(ibin + 1) <= 0.) continue;
+    pthigh = hin->GetBinLowEdge(ibin + 2);
+    break;
+  }
+
+  Double_t mean = hin->GetMean();
+  Double_t maxdiff = 0.;
+  TH1 *hmax = NULL;
+  for (Int_t inode = 0; inode < hin->GetNbinsX(); inode++) {
+
+    Double_t ptnode = hin->GetBinCenter(inode + 1);
+    TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehard", hin->GetName()));
+    
+    for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+      if (hin->GetBinError(ibin + 1) <= 0.) continue;
+      Double_t val = hin->GetBinContent(ibin + 1);
+      Double_t err = hin->GetBinError(ibin + 1);
+      Double_t cen = hin->GetBinCenter(ibin + 1);
+      if (cen < ptnode)
+        err *= -1. + (cen - ptlow) / (ptnode - ptlow);
+      else
+        err *= (cen - ptnode) / (pthigh - ptnode);
+
+      hout->SetBinContent(ibin + 1, val + sign * err);
+    }
+
+    Double_t diff = TMath::Abs(mean - hout->GetMean());
+    if (diff > maxdiff) {
+      //      printf("found max at %f\n", ptnode);
+      if (hmax) delete hmax;
+      hmax = (TH1 *)hout->Clone("hmax");
+      maxdiff = diff;
+    }
+    delete hout;
+  }
+  return hmax;
+}
+
+YieldMean_IntegralMean(TH1 *hdata, TH1 *hlo, TH1 *hhi, Double_t &integral, Double_t &mean)
+{
+  
+  /*
+   * compute integrals
+   */
+
+  Double_t cont, err, width, cent;
+  Double_t I = 0., IX = 0., Ierr = 0., IXerr = 0., Ilerr = 0., IXlerr = 0.;
+  Double_t M = 0., Merr = 0., Mlerr = 0., C;
+
+  /* integrate the data */
+  for (Int_t ibin = 0; ibin < hdata->GetNbinsX(); ibin++) {
+    cent = hdata->GetBinCenter(ibin + 1);
+    width = hdata->GetBinWidth(ibin + 1);
+    cont = width * hdata->GetBinContent(ibin + 1);
+    err = width * hdata->GetBinError(ibin + 1);
+    if (err <= 0.) continue;
+    I += cont;
+    IX += cont * cent;
+  }
+  /* integrate low */
+  for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) {
+    cent = hlo->GetBinCenter(ibin + 1);
+    width = hlo->GetBinWidth(ibin + 1);
+    cont = width * hlo->GetBinContent(ibin + 1);
+    err = width * hlo->GetBinError(ibin + 1);
+    if (err <= 0.) continue;
+    I += cont;
+    IX += cont * cent;
+  }
+  /* integrate high */
+  for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) {
+    cent = hhi->GetBinCenter(ibin + 1);
+    width = hhi->GetBinWidth(ibin + 1);
+    cont = width * hhi->GetBinContent(ibin + 1);
+    err = width * hhi->GetBinError(ibin + 1);
+    if (err <= 0.) continue;
+    I += cont;
+    IX += cont * cent;
+  }
+
+  /* set values */
+  integral = I;
+  mean = IX / I;
+}