.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / 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