]>
Commit | Line | Data |
---|---|---|
1f1d6be7 | 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 */ | |
4070f709 | 1130 | gSystem->Load("libANALYSIS"); |
1131 | gSystem->Load("libANALYSISalice"); | |
1132 | gSystem->Load("libCORRFW"); | |
1133 | gSystem->Load("libPWGTools"); | |
1f1d6be7 | 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 |