]>
Commit | Line | Data |
---|---|---|
c3fc8cd6 | 1 | #ifndef ALILANDAUGAUS_H |
2 | #define ALILANDAUGAUS_H | |
3 | /** | |
4 | * @file AliLandauGaus.h | |
5 | * @author Christian Holm Christensen <cholm@nbi.dk> | |
6 | * @date Tue Mar 11 08:53:47 2014 | |
7 | * | |
8 | * @brief Declaration and implementation of Landau-Gauss distributions. | |
9 | * | |
10 | * @ingroup pwglf_forward | |
11 | */ | |
12 | ||
13 | #include <TObject.h> | |
14 | #include <TF1.h> | |
15 | #include <TMath.h> | |
16 | ||
17 | /** | |
18 | * This class contains static member functions to calculate the energy | |
19 | * loss stragling - most notably the N-particle energy loss as a sum | |
20 | * of convolutions of a Landau and Gauss distribution. | |
21 | * | |
22 | * That is, for a single particle we have the function @f$ f(x)@f$: | |
23 | * | |
24 | * @f[ | |
25 | * f(x;\Delta_p,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}} | |
26 | * \int_{-\infty}^{+\infty} dx' f'_{L}(x',\Delta_p,\xi) | |
27 | * \exp{-\frac{(x-x')^2}{2\sigma'^2}}} | |
28 | * @f] | |
29 | * | |
30 | * where @f$ f'_{L}@f$ is the Landau distribution, @f$\Delta_p@f$ the | |
31 | * most probable energy loss, @f$ \xi@f$ the width of the Landau, and | |
32 | * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$. Here, @f$\sigma@f$ is the | |
33 | * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter | |
34 | * modelling noise in the detector. | |
35 | * | |
36 | * for @f$ i@f$ particles this is modified to | |
37 | * | |
38 | * @f[ | |
39 | * f_i(x;\Delta_{p},\xi,\sigma')=f(x;\Delta_{p,i},\xi_i,\sigma_i') | |
40 | * @f] | |
41 | * | |
42 | * corresponding to @f$ i@f$ particles i.e., with the substitutions | |
43 | * @f{eqnarray*}{ | |
44 | * \Delta_p \rightarrow \Delta_{p,i}&=& i(\Delta_p + \xi\log(i))\\ | |
45 | * \xi \rightarrow \xi_i &=& i \xi\\ | |
46 | * \sigma \rightarrow \sigma_i &=& \sqrt{i}\sigma\\ | |
47 | * \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2 | |
48 | * @f} | |
49 | * | |
50 | * Because of the convolution with a Gaussian, the most-probable-value | |
51 | * @f$\Delta_p'@f$ of the resulting distribution is not really at the | |
52 | * Landau most-probable-value @f$\Delta_p@f$. In fact we find that | |
53 | * @f$\Delta_p' > \Delta_p@f$. | |
54 | * | |
55 | * Ideally, one would find an analytic expression for this shift by | |
56 | * solving | |
57 | * | |
58 | * @f{eqnarray} | |
59 | * 0 &=& \frac{d f_i(x;\Delta_p,\xi,\sigma)}{d\Delta}\\ | |
60 | * &=& \frac{d}{d\Delta}\frac{1}{\sigma' \sqrt{2 \pi}} | |
61 | * \int_{-\infty}^{+\infty} dx' f'_{L}(x',\Delta_p,\xi) | |
62 | * \exp{-\frac{(x-x')^2}{2\sigma'^2}}} | |
63 | * @f{eqnarray} | |
64 | * | |
65 | * for @f$x@f$ as a function of @f$\Delta_p,\xi,\sigma,i@f$. However, | |
66 | * do to the complex nature of the Landau distribution this is not | |
67 | * really feasible. | |
68 | * | |
69 | * Instead, the shift was studied numerically. Landau-Gauss | |
70 | * distributions for @f$i=1,\ldots@f$ where generated with varying | |
71 | * @f$\xi@f$ and @f$\sigma@f$. The distributions was then numerically | |
72 | * differentiated and the root @f$\Delta_p'@f$ of that derivative | |
73 | * found numerically. The difference | |
74 | * @f$\delta\Delta_p=\Delta_p'-\Delta_p@f$ was then studied as a | |
75 | * function of the @f$\sigma,\xi@f$ parameters and an approximate | |
76 | * expression was found | |
77 | * | |
78 | * @f[ | |
79 | * \delta\Delta_p \approx \frac{c \sigma u}{(1+1/i)^{p u^{3/2}}} | |
80 | * @f] | |
81 | * | |
82 | * where @f$ u=\sigma/\xi@f$. The parameters @f$c@f$ and @f$p@f$ is | |
83 | * found to depend on @f$ u@f$ only weakly, and for practical | |
84 | * applications where @f$u\approx1@f$, we set @f$ c=p=1/2@f$. | |
85 | * | |
86 | * For the evaluating the full energy loss distribution from f@$ | |
87 | * 1+2+\ldots,n@f$ particles, we evaluate | |
88 | * | |
89 | * @f[ | |
90 | * f_N(x;\Delta_p,\xi,\sigma')=\sum_{i=1}^N a_i f_i(x;\Delta_p,\xi,\sigma',a) | |
91 | * @f] | |
92 | * | |
93 | * where @f$ f(x;\Delta_p,\xi,\sigma')@f$ is the convolution of a | |
94 | * Landau with a Gaussian (see LandauGaus), and @f$ a@f$ is a vector of | |
95 | * weights for each @f$ f_i@f$. Note that @f$ a_1 = 1@f$. | |
96 | * | |
97 | * Everything is defined in this header file to make it easy to move | |
98 | * this code around. Nothing here's meant to be persistent, so we | |
99 | * can easily do that. | |
100 | * | |
101 | * References: | |
102 | * - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5"> | |
103 | * Nucl.Instrum.Meth.B1:16</a> | |
104 | * - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a> | |
105 | * - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html"> | |
106 | * ROOT implementation</a> | |
107 | * | |
108 | * @ingroup pwglf_forward | |
109 | */ | |
110 | class AliLandauGaus | |
111 | { | |
112 | public: | |
113 | /** Enumeration of parameters */ | |
114 | enum { | |
115 | kC = 0, | |
116 | kDelta, | |
117 | kXi, | |
118 | kSigma, | |
119 | kSigmaN, | |
120 | kN, | |
121 | kA | |
122 | }; | |
123 | /** Enumeration of colors */ | |
124 | ||
125 | //__________________________________________________________________ | |
126 | /** | |
127 | * @{ | |
128 | * @name Constants | |
129 | */ | |
130 | //------------------------------------------------------------------ | |
131 | /** | |
132 | * The shift of the most probable value for the ROOT function | |
133 | * TMath::Landau | |
134 | * | |
135 | * @return Shift of TMath::Landau | |
136 | */ | |
137 | static Double_t MPShift() { return -0.22278298; } | |
138 | /** | |
139 | * Constant of @f$\sigma@f$ shift | |
140 | * | |
141 | * @return Constant of @f$\sigma@f$ shift | |
142 | */ | |
143 | static Double_t SigmaShiftC() { return 0.5446; } | |
144 | /** | |
145 | * Power factor of @f$\sigma@f$ shift | |
146 | * | |
147 | * @return Power factor of @f$\sigma@f$ shift | |
148 | */ | |
149 | static Double_t SigmaShiftP() { return 0.5746; } | |
150 | /** | |
151 | * Normalization constant | |
152 | * | |
153 | * @return The landau-gauss normalization constant | |
154 | */ | |
155 | static Double_t InvSq2Pi() { return 1. / TMath::Sqrt(TMath::TwoPi()); } | |
156 | /** | |
157 | * How many sigma's of the Gaussian in the Landau, Gaussian | |
158 | * convolution to integrate over | |
159 | */ | |
160 | static Double_t NSigma() { return 5; } | |
161 | /** | |
162 | * Number of steps to do in the Landau, Gaussiam convolution | |
163 | */ | |
164 | static Int_t NSteps() { return 100; } | |
165 | /* @} */ | |
166 | ||
167 | //__________________________________________________________________ | |
168 | /** | |
169 | * @{ | |
170 | * @name Function calculations | |
171 | */ | |
172 | //------------------------------------------------------------------ | |
173 | /** | |
174 | * Calculate the shifted Landau | |
175 | * @f[ | |
176 | * f'_{L}(x;\Delta_p,\xi) = f_L(x;\Delta_p+0.22278298\xi) | |
177 | * @f] | |
178 | * | |
179 | * where @f$ f_{L}@f$ is the ROOT implementation of the Landau | |
180 | * distribution (known to have @f$\Delta_{p}=-0.22278298@f$ for | |
181 | * @f$\Delta_p=0,\xi=1@f$. | |
182 | * | |
183 | * @param x Where to evaluate @f$ f'_{L}@f$ | |
184 | * @param delta Most probable value | |
185 | * @param xi The 'width' of the distribution | |
186 | * | |
187 | * @return @f$ f'_{L}(x;\Delta,\xi) @f$ | |
188 | */ | |
189 | static Double_t Fl(Double_t x, Double_t delta, Double_t xi); | |
190 | //------------------------------------------------------------------ | |
191 | /** | |
192 | * Calculate the value of a Landau convolved with a Gaussian | |
193 | * | |
194 | * @f[ | |
195 | * f(x;\Delta_p,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}} | |
196 | * \int_{-\infty}^{+\infty} dx' f'_{L}(x';\Delta_p,\xi) | |
197 | * \exp{-\frac{(x-x')^2}{2\sigma'^2}} | |
198 | * @f] | |
199 | * | |
200 | * Note that this function uses the constants NSteps() and | |
201 | * NSigma() | |
202 | * | |
203 | * @param x where to evaluate @f$ f@f$ | |
204 | * @param delta @f$ \Delta_p@f$ of @f$ f(x;\Delta_p,\xi,\sigma')@f$ | |
205 | * @param xi @f$ \xi@f$ of @f$ f(x;\Delta_p,\xi,\sigma')@f$ | |
206 | * @param sigma @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$ | |
207 | * @param sigma_n @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$ | |
208 | * | |
209 | * @return @f$ f@f$ evaluated at @f$ x@f$. | |
210 | */ | |
211 | static Double_t F(Double_t x, Double_t delta, Double_t xi, | |
212 | Double_t sigma, Double_t sigma_n); | |
213 | //------------------------------------------------------------------ | |
214 | /** | |
215 | * Evaluate | |
216 | * @f[ | |
217 | * f_i(x;\Delta_p,\xi,\sigma') = f(x;\Delta_{p,i},\xi_i,\sigma_i') | |
218 | * @f] | |
219 | * corresponding to @f$ i@f$ particles. | |
220 | * | |
221 | * @param x Where to evaluate | |
222 | * @param delta @f$ \Delta@f$ | |
223 | * @param xi @f$ \xi@f$ | |
224 | * @param sigma @f$ \sigma@f$ | |
225 | * @param sigma_n @f$ \sigma_n@f$ | |
226 | * @param i @f$ i @f$ | |
227 | * | |
228 | * @return @f$ f_i @f$ evaluated | |
229 | */ | |
230 | static Double_t Fi(Double_t x, Double_t delta, Double_t xi, | |
231 | Double_t sigma, Double_t sigma_n, Int_t i); | |
232 | ||
233 | //------------------------------------------------------------------ | |
234 | /** | |
235 | * Numerically evaluate | |
236 | * @f[ | |
237 | * \left.\frac{\partial f_i}{\partial p_i}\right|_{x} | |
238 | * @f] | |
239 | * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter. The mapping | |
240 | * of the parameters is given by | |
241 | * | |
242 | * - 0: @f$\Delta@f$ | |
243 | * - 1: @f$\xi@f$ | |
244 | * - 2: @f$\sigma@f$ | |
245 | * - 3: @f$\sigma_n@f$ | |
246 | * | |
247 | * @param x Where to evaluate | |
248 | * @param ipar Parameter number | |
249 | * @param dp @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$ | |
250 | * @param delta @f$ \Delta_p@f$ | |
251 | * @param xi @f$ \xi@f$ | |
252 | * @param sigma @f$ \sigma@f$ | |
253 | * @param sigma_n @f$ \sigma_n@f$ | |
254 | * @param i @f$ i@f$ | |
255 | * | |
256 | * @return @f$ f_i@f$ evaluated | |
257 | */ | |
258 | static Double_t DFidPar(Double_t x, UShort_t ipar, Double_t dp, | |
259 | Double_t delta, Double_t xi, | |
260 | Double_t sigma, Double_t sigma_n, Int_t i); | |
261 | ||
262 | //------------------------------------------------------------------ | |
263 | /** | |
264 | * Evaluate | |
265 | * @f[ | |
266 | * f_N(x;\Delta_p,\xi,\sigma')=\sum_{i=1}^N a_i f_i(x;\Delta_p,\xi,\sigma',a) | |
267 | * @f] | |
268 | * | |
269 | * where @f$ f(x;\Delta_p,\xi,\sigma')@f$ is the convolution of a | |
270 | * Landau with a Gaussian (see LandauGaus). | |
271 | * | |
272 | * @param x Where to evaluate @f$ f_N@f$ | |
273 | * @param delta @f$ \Delta_1@f$ | |
274 | * @param xi @f$ \xi_1@f$ | |
275 | * @param sigma @f$ \sigma_1@f$ | |
276 | * @param sigma_n @f$ \sigma_n@f$ | |
277 | * @param n @f$ N@f$ in the sum above. | |
278 | * @param a Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for | |
279 | * @f$ i > 1@f$ | |
280 | * | |
281 | * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$ | |
282 | */ | |
283 | static Double_t Fn(Double_t x, Double_t delta, Double_t xi, | |
284 | Double_t sigma, Double_t sigma_n, Int_t n, | |
285 | const Double_t* a); | |
286 | /** | |
287 | * Get parameters for the @f$ i@f$ particle response. | |
288 | * | |
289 | * @f{eqnarray*}{ | |
290 | * \Delta \rightarrow \Delta_i &=& i(\Delta + \xi\log(i))\\ | |
291 | * \xi \rightarrow \xi_i &=& i \xi\\ | |
292 | * \sigma \rightarrow \sigma_i &=& \sqrt{i}\sigma\\ | |
293 | * \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2 | |
294 | * @f} | |
295 | * | |
296 | * | |
297 | * @param i Number of particles | |
298 | * @param delta Input: single particle, output @f$ i@f$ particle | |
299 | * @param xi Input: single particle, output @f$ i@f$ particle | |
300 | * @param sigma Input: single particle, output @f$ i@f$ particle | |
301 | */ | |
302 | static void IPars(Int_t i, Double_t& delta, Double_t& xi, Double_t& sigma); | |
303 | /** | |
304 | * Set and check if sigma shift is enabled | |
305 | * | |
306 | * @param val if <0, then only check. Otherwise set enabled (>0) or not (=0) | |
307 | * | |
308 | * @return whether the sigma shift is enabled or not | |
309 | */ | |
310 | static Bool_t EnableSigmaShift(Short_t val=-1); | |
311 | /** | |
312 | * Get the shift of the MPV due to convolution with a Gaussian. | |
313 | * | |
314 | * @f[ | |
315 | * \delta\Delta_p \approx \frac{c \sigma u}{(1+1/i)^{p u^{3/2}}} | |
316 | * @f] | |
317 | * | |
318 | * where @f$ u=\sigma/\xi@f$. | |
319 | * | |
320 | * @param i Number of particles | |
321 | * @param xi Landau width @f$\xi@f$ | |
322 | * @param sigma Gaussian variance @f$\sigma@f$ | |
323 | * | |
324 | * @return The shift | |
325 | */ | |
326 | static Double_t SigmaShift(Int_t i, Double_t xi, Double_t sigma); | |
327 | /* @} */ | |
328 | ||
329 | ||
330 | //__________________________________________________________________ | |
331 | /** | |
332 | * @{ | |
333 | * @name Utilities for defining TF1 objects | |
334 | */ | |
335 | //------------------------------------------------------------------ | |
336 | /** | |
337 | * Generate a TF1 object of @f$ f_1@f$ | |
338 | * | |
339 | * @param c Constant | |
340 | * @param delta @f$ \Delta_1@f$ | |
341 | * @param xi @f$ \xi_1@f$ | |
342 | * @param sigma @f$ \sigma_1@f$ | |
343 | * @param sigma_n @f$ \sigma_n@f$ | |
344 | * @param xmin Least value of range | |
345 | * @param xmax Largest value of range | |
346 | * | |
347 | * @return Newly allocated TF1 object | |
348 | */ | |
349 | static TF1* MakeF1(Double_t c, | |
350 | Double_t delta, Double_t xi, | |
351 | Double_t sigma, Double_t sigma_n, | |
352 | Double_t xmin, Double_t xmax); | |
353 | //------------------------------------------------------------------ | |
354 | /** | |
355 | * Generate a TF1 object of @f$ f_I@f$ | |
356 | * | |
357 | * @param c Constant | |
358 | * @param delta @f$ \Delta_1@f$ | |
359 | * @param xi @f$ \xi_1@f$ | |
360 | * @param sigma @f$ \sigma_1@f$ | |
361 | * @param sigma_n @f$ \sigma_n@f$ | |
362 | * @param i @f$ i@f$ - the number of particles | |
363 | * @param xmin Least value of range | |
364 | * @param xmax Largest value of range | |
365 | * | |
366 | * @return Newly allocated TF1 object | |
367 | */ | |
368 | static TF1* MakeFi(Double_t c, | |
369 | Double_t delta, Double_t xi, | |
370 | Double_t sigma, Double_t sigma_n, | |
371 | Int_t i, | |
372 | Double_t xmin, Double_t xmax); | |
373 | //------------------------------------------------------------------ | |
374 | /** | |
375 | * Generate a TF1 object of @f$ f_N@f$ | |
376 | * | |
377 | * @param c Constant | |
378 | * @param delta @f$ \Delta_1@f$ | |
379 | * @param xi @f$ \xi_1@f$ | |
380 | * @param sigma @f$ \sigma_1@f$ | |
381 | * @param sigma_n @f$ \sigma_n@f$ | |
382 | * @param n @f$ N@f$ - how many particles to sum to | |
383 | * @param a Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for | |
384 | * @f$ i > 1@f$ | |
385 | * @param xmin Least value of range | |
386 | * @param xmax Largest value of range | |
387 | * | |
388 | * @return Newly allocated TF1 object | |
389 | */ | |
390 | static TF1* MakeFn(Double_t c, | |
391 | Double_t delta, Double_t xi, | |
392 | Double_t sigma, Double_t sigma_n, | |
393 | Int_t n, const Double_t* a, | |
394 | Double_t xmin, Double_t xmax); | |
395 | /** | |
396 | * Make a TF1 object that describes the convoluted energy loss of | |
397 | * primary and secondary particles. | |
398 | * | |
399 | * @param c1 Primary weight | |
400 | * @param delta Most-probable value | |
401 | * @param xi1 Primary width | |
402 | * @param sigma Gaussian variance | |
403 | * @param c2 Secondary weight | |
404 | * @param xi2 Secondary with | |
405 | * | |
406 | * @return Pointer to newly allocated TF1 object | |
407 | */ | |
408 | static TF1* MakeComposite(Double_t c1, | |
409 | Double_t delta, | |
410 | Double_t xi1, | |
411 | Double_t sigma, | |
412 | Double_t c2, | |
413 | Double_t xi2, | |
414 | Double_t xmin, | |
415 | Double_t xmax); | |
416 | ||
417 | //------------------------------------------------------------------ | |
418 | /** | |
419 | * Get the color of the @f$ i@f$ particle response | |
420 | * | |
421 | * @param i Particle number | |
422 | * | |
423 | * @return Color | |
424 | */ | |
425 | static Color_t GetIColor(Int_t i); | |
426 | //------------------------------------------------------------------ | |
427 | /** | |
428 | * Utility function for TF1 definition | |
429 | * | |
430 | * @param pp Array of parameters | |
431 | * @param xp Pointer to independent variable | |
432 | * | |
433 | * @see AliLandauGaus::F | |
434 | * | |
435 | * @return Landau convolved with a Gauss | |
436 | */ | |
437 | static Double_t F1Func(Double_t* xp, Double_t* pp); | |
438 | //------------------------------------------------------------------ | |
439 | /** | |
440 | * Utility function for TF1 definition | |
441 | * | |
442 | * @param pp Array of parameters | |
443 | * @param xp Pointer to independent variable | |
444 | * | |
445 | * @see AliLandauGaus::Fn | |
446 | * | |
447 | * @return Landau convolved with a Gauss | |
448 | */ | |
449 | static Double_t FnFunc(Double_t* xp, Double_t* pp); | |
450 | //------------------------------------------------------------------ | |
451 | /** | |
452 | * Utility function for TF1 definition | |
453 | * | |
454 | * @param pp Array of parameters | |
455 | * @param xp Pointer to independent variable | |
456 | * | |
457 | * @see AliLandauGaus::Fn | |
458 | * | |
459 | * @return Landau convolved with a Gauss | |
460 | */ | |
461 | static Double_t FiFunc(Double_t* xp, Double_t* pp); | |
462 | //------------------------------------------------------------------ | |
463 | /** | |
464 | * Utility function for TF1 definition | |
465 | * | |
466 | * @param pp Array of parameters | |
467 | * @param xp Pointer to independent variable | |
468 | * | |
469 | * @see AliLandauGaus::F | |
470 | * | |
471 | * @return Landau convolved with a Gauss | |
472 | */ | |
473 | static Double_t CompFunc(Double_t* xp, Double_t* pp); | |
474 | /* @} */ | |
475 | }; | |
476 | //____________________________________________________________________ | |
477 | inline Bool_t | |
478 | AliLandauGaus::EnableSigmaShift(Short_t val) | |
479 | { | |
480 | static Bool_t enabled = true; | |
481 | if (val >= 0) enabled = val == 1; | |
482 | return enabled; | |
483 | } | |
484 | //____________________________________________________________________ | |
485 | inline void | |
486 | AliLandauGaus::IPars(Int_t i, Double_t& delta, Double_t& xi, Double_t& sigma) | |
487 | { | |
488 | #ifndef NO_SIGMA_SHIFT | |
489 | Double_t dDelta = EnableSigmaShift() ? SigmaShift(i, xi, sigma) : 0; | |
490 | #else | |
491 | // This is for testing | |
492 | Double_t dDelta = 0; | |
493 | #endif | |
494 | if (i == 1) { | |
495 | delta -= dDelta; | |
496 | return; | |
497 | } | |
498 | ||
499 | delta = i * (delta + xi * TMath::Log(i)) - dDelta; | |
500 | xi = i * xi; | |
501 | sigma = TMath::Sqrt(Double_t(i)) * sigma; | |
502 | } | |
503 | //____________________________________________________________________ | |
504 | inline Double_t | |
505 | AliLandauGaus::SigmaShift(Int_t i, Double_t xi, Double_t sigma) | |
506 | { | |
507 | if (xi <= 0) return 0; | |
21676077 | 508 | if (sigma <= 0) return 0; |
c3fc8cd6 | 509 | const Double_t c = SigmaShiftC(); |
510 | const Double_t p = SigmaShiftP(); | |
511 | const Double_t u = sigma / xi; | |
512 | const Double_t q = p*u*TMath::Sqrt(u); | |
513 | if (q > 100) return 0; | |
514 | return c * sigma / TMath::Power(1+1./i, q); | |
515 | } | |
516 | //____________________________________________________________________ | |
517 | inline Double_t | |
518 | AliLandauGaus::Fl(Double_t x, Double_t delta, Double_t xi) | |
519 | { | |
520 | Double_t deltaP = delta - xi * MPShift(); | |
521 | return TMath::Landau(x, deltaP, xi, true); | |
522 | } | |
523 | //____________________________________________________________________ | |
524 | inline Double_t | |
525 | AliLandauGaus::F(Double_t x, Double_t delta, Double_t xi, | |
526 | Double_t sigma, Double_t sigmaN) | |
527 | { | |
528 | if (xi <= 0) return 0; | |
529 | ||
530 | const Int_t nSteps = NSteps(); | |
531 | const Double_t nSigma = NSigma(); | |
532 | const Double_t deltaP = delta; // - sigma * sigmaShift; // + sigma * mpshift; | |
533 | const Double_t sigma2 = sigmaN*sigmaN + sigma*sigma; | |
534 | const Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2); | |
535 | const Double_t xlow = x - nSigma * sigma1; | |
536 | const Double_t xhigh = x + nSigma * sigma1; | |
537 | const Double_t step = (xhigh - xlow) / nSteps; | |
538 | Double_t sum = 0; | |
539 | ||
540 | for (Int_t i = 0; i <= nSteps/2; i++) { | |
541 | const Double_t x1 = xlow + (i - .5) * step; | |
542 | const Double_t x2 = xhigh - (i - .5) * step; | |
543 | sum += Fl(x1, deltaP, xi) * TMath::Gaus(x, x1, sigma1); | |
544 | sum += Fl(x2, deltaP, xi) * TMath::Gaus(x, x2, sigma1); | |
545 | } | |
546 | return step * sum * InvSq2Pi() / sigma1; | |
547 | } | |
548 | ||
549 | //____________________________________________________________________ | |
550 | inline Double_t | |
551 | AliLandauGaus::Fi(Double_t x, Double_t delta, Double_t xi, | |
552 | Double_t sigma, Double_t sigmaN, Int_t i) | |
553 | { | |
554 | Double_t deltaI = delta; | |
555 | Double_t xiI = xi; | |
556 | Double_t sigmaI = sigma; | |
557 | IPars(i, deltaI, xiI, sigmaI); | |
558 | if (sigmaI < 1e-10) | |
559 | // Fall back to landau | |
560 | return Fl(x, deltaI, xiI); | |
561 | ||
562 | return F(x, deltaI, xiI, sigmaI, sigmaN); | |
563 | } | |
564 | //____________________________________________________________________ | |
565 | inline Double_t | |
566 | AliLandauGaus::Fn(Double_t x, Double_t delta, Double_t xi, | |
567 | Double_t sigma, Double_t sigmaN, Int_t n, | |
568 | const Double_t* a) | |
569 | { | |
570 | Double_t result = Fi(x, delta, xi, sigma, sigmaN, 1); | |
571 | for (Int_t i = 2; i <= n; i++) | |
572 | result += a[i-2] * Fi(x,delta,xi,sigma,sigmaN,i); | |
573 | return result; | |
574 | } | |
575 | ||
576 | //____________________________________________________________________ | |
577 | inline Double_t | |
578 | AliLandauGaus::DFidPar(Double_t x, | |
579 | UShort_t par, Double_t dPar, | |
580 | Double_t delta, Double_t xi, | |
581 | Double_t sigma, Double_t sigmaN, | |
582 | Int_t i) | |
583 | { | |
584 | if (dPar == 0) return 0; | |
585 | Double_t dp = dPar; | |
586 | Double_t d2 = dPar / 2; | |
587 | Double_t deltaI = i * (delta + xi * TMath::Log(i)); | |
588 | Double_t xiI = i * xi; | |
589 | Double_t si = TMath::Sqrt(Double_t(i)); | |
590 | Double_t sigmaI = si*sigma; | |
591 | Double_t y1 = 0; | |
592 | Double_t y2 = 0; | |
593 | Double_t y3 = 0; | |
594 | Double_t y4 = 0; | |
595 | switch (par) { | |
596 | case 0: | |
597 | y1 = Fi(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i); | |
598 | y2 = Fi(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i); | |
599 | y3 = Fi(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i); | |
600 | y4 = Fi(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i); | |
601 | break; | |
602 | case 1: | |
603 | y1 = Fi(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i); | |
604 | y2 = Fi(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i); | |
605 | y3 = Fi(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i); | |
606 | y4 = Fi(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i); | |
607 | break; | |
608 | case 2: | |
609 | y1 = Fi(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i); | |
610 | y2 = Fi(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i); | |
611 | y3 = Fi(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i); | |
612 | y4 = Fi(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i); | |
613 | break; | |
614 | case 3: | |
615 | y1 = Fi(x, deltaI, xiI, sigmaI, sigmaN+dp, i); | |
616 | y2 = Fi(x, deltaI, xiI, sigmaI, sigmaN+d2, i); | |
617 | y3 = Fi(x, deltaI, xiI, sigmaI, sigmaN-d2, i); | |
618 | y4 = Fi(x, deltaI, xiI, sigmaI, sigmaN-dp, i); | |
619 | break; | |
620 | default: | |
621 | return 0; | |
622 | } | |
623 | ||
624 | Double_t d0 = y1 - y4; | |
625 | Double_t d1 = 2 * (y2 - y3); | |
626 | ||
627 | Double_t g = 1/(2*dp) * (4*d1 - d0) / 3; | |
628 | ||
629 | return g; | |
630 | } | |
631 | ||
632 | ||
633 | //____________________________________________________________________ | |
634 | inline Color_t | |
635 | AliLandauGaus::GetIColor(Int_t i) | |
636 | { | |
637 | const Int_t kColors[] = { kRed+1, | |
638 | kPink+3, | |
639 | kMagenta+2, | |
640 | kViolet+2, | |
641 | kBlue+1, | |
642 | kAzure+3, | |
643 | kCyan+1, | |
644 | kTeal+2, | |
645 | kGreen+2, | |
646 | kSpring+3, | |
647 | kYellow+2, | |
648 | kOrange+2 }; | |
649 | return kColors[((i-1) % 12)]; | |
650 | } | |
651 | //____________________________________________________________________ | |
652 | inline Double_t | |
653 | AliLandauGaus::F1Func(Double_t* xp, Double_t* pp) | |
654 | { | |
655 | Double_t x = xp[0]; | |
656 | Double_t constant = pp[kC]; | |
657 | Double_t delta = pp[kDelta]; | |
658 | Double_t xi = pp[kXi]; | |
659 | Double_t sigma = pp[kSigma]; | |
660 | Double_t sigmaN = pp[kSigmaN]; | |
661 | ||
662 | return constant * F(x, delta, xi, sigma, sigmaN); | |
663 | } | |
664 | //____________________________________________________________________ | |
665 | inline Double_t | |
666 | AliLandauGaus::FiFunc(Double_t* xp, Double_t* pp) | |
667 | { | |
668 | Double_t x = xp[0]; | |
669 | Double_t constant = pp[kC]; | |
670 | Double_t delta = pp[kDelta]; | |
671 | Double_t xi = pp[kXi]; | |
672 | Double_t sigma = pp[kSigma]; | |
673 | Double_t sigmaN = pp[kSigmaN]; | |
674 | Int_t i = Int_t(pp[kN]); | |
675 | ||
676 | return constant * Fi(x, delta, xi, sigma, sigmaN, i); | |
677 | } | |
678 | //____________________________________________________________________ | |
679 | inline Double_t | |
680 | AliLandauGaus::FnFunc(Double_t* xp, Double_t* pp) | |
681 | { | |
682 | Double_t x = xp[0]; | |
683 | Double_t constant = pp[kC]; | |
684 | Double_t delta = pp[kDelta]; | |
685 | Double_t xi = pp[kXi]; | |
686 | Double_t sigma = pp[kSigma]; | |
687 | Double_t sigmaN = pp[kSigmaN]; | |
688 | Int_t n = Int_t(pp[kN]); | |
689 | Double_t* a = &(pp[kA]); | |
690 | ||
691 | return constant * Fn(x, delta, xi, sigma, sigmaN, n, a); | |
692 | } | |
693 | //____________________________________________________________________ | |
694 | inline Double_t | |
695 | AliLandauGaus::CompFunc(Double_t* xp, Double_t* pp) | |
696 | { | |
697 | Double_t x = xp[0]; | |
698 | Double_t cP = pp[kC]; | |
699 | Double_t deltaP = pp[kDelta]; | |
700 | Double_t xiP = pp[kXi]; | |
701 | Double_t sigmaP = pp[kSigma]; | |
702 | Double_t cS = pp[kSigma+1]; | |
703 | Double_t deltaS = deltaP; // pp[kSigma+2]; | |
704 | Double_t xiS = pp[kSigma+2/*3*/]; | |
705 | Double_t sigmaS = sigmaP; // pp[kSigma+4]; | |
706 | ||
707 | return (cP * F(x,deltaP,xiP,sigmaP,0) + | |
708 | cS * F(x,deltaS,xiS,sigmaS,0)); | |
709 | ||
710 | } | |
711 | //____________________________________________________________________ | |
712 | inline TF1* | |
713 | AliLandauGaus::MakeF1(Double_t c, | |
714 | Double_t delta, Double_t xi, | |
715 | Double_t sigma, Double_t sigmaN, | |
716 | Double_t xmin, Double_t xmax) | |
717 | { | |
718 | // Define the function to fit | |
719 | TF1* f = new TF1("landau1", &F1Func, xmin,xmax,kSigmaN+1); | |
720 | ||
721 | // Set initial guesses, parameter names, and limits | |
722 | f->SetParameters(c,delta,xi,sigma,sigmaN); | |
723 | f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}"); | |
724 | f->SetNpx(500); | |
725 | f->SetLineColor(GetIColor(1)); | |
726 | ||
727 | return f; | |
728 | } | |
729 | ||
730 | //____________________________________________________________________ | |
731 | inline TF1* | |
732 | AliLandauGaus::MakeFn(Double_t c, | |
733 | Double_t delta, Double_t xi, | |
734 | Double_t sigma, Double_t sigmaN, Int_t n, | |
735 | const Double_t* a, | |
736 | Double_t xmin, Double_t xmax) | |
737 | { | |
738 | Int_t npar = kN+n; | |
739 | TF1* f = new TF1(Form("nlandau%d", n), &FnFunc,xmin,xmax,npar); | |
740 | f->SetLineColor(GetIColor(n)); | |
741 | f->SetLineWidth(2); | |
742 | f->SetNpx(500); | |
743 | f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N"); | |
744 | f->SetParameter(kC, c); | |
745 | f->SetParameter(kDelta, delta); | |
746 | f->SetParameter(kXi, xi); | |
747 | f->SetParameter(kSigma, sigma); | |
748 | f->SetParameter(kSigmaN, sigmaN); | |
749 | f->FixParameter(kN, n); | |
750 | for (UShort_t i = 2; i <= n; i++) { | |
751 | f->SetParameter(kA+i-2, a[i-2]); | |
752 | f->SetParName(kA+i-2, Form("a_{%d}", i)); | |
753 | } | |
754 | return f; | |
755 | } | |
756 | //____________________________________________________________________ | |
757 | inline TF1* | |
758 | AliLandauGaus::MakeFi(Double_t c, | |
759 | Double_t delta, Double_t xi, | |
760 | Double_t sigma, Double_t sigmaN, Int_t i, | |
761 | Double_t xmin, Double_t xmax) | |
762 | { | |
763 | Int_t npar = kN+1; | |
764 | TF1* f = new TF1(Form("ilandau%d", i), &FiFunc,xmin,xmax,npar); | |
765 | f->SetLineColor(GetIColor(i)); | |
766 | f->SetLineWidth(1); | |
767 | f->SetNpx(500); | |
768 | f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i"); | |
769 | f->SetParameter(kC, c); | |
770 | f->SetParameter(kDelta, delta); | |
771 | f->SetParameter(kXi, xi); | |
772 | f->SetParameter(kSigma, sigma); | |
773 | f->SetParameter(kSigmaN, sigmaN); | |
774 | f->FixParameter(kN, i); | |
775 | ||
776 | return f; | |
777 | } | |
778 | //____________________________________________________________________ | |
779 | inline TF1* | |
780 | AliLandauGaus::MakeComposite(Double_t c1, | |
781 | Double_t delta, | |
782 | Double_t xi1, | |
783 | Double_t sigma, | |
784 | Double_t c2, | |
785 | Double_t xi2, | |
786 | Double_t xmin, | |
787 | Double_t xmax) | |
788 | { | |
789 | TF1* comp = new TF1("composite", &CompFunc, xmin, xmax, kSigma+1+2); | |
790 | comp->SetParNames("C", "#Delta_{p}", "#xi", "#sigma", | |
791 | "C#prime", "#xi#prime"); | |
792 | comp->SetParameters(c1, // 0 Primary weight | |
793 | delta, // 1 Primary Delta | |
794 | xi1, // 2 primary Xi | |
795 | sigma, // 3 primary sigma | |
796 | c2, // 5 Secondary weight | |
797 | xi2); // 6 secondary Xi | |
798 | return comp; | |
799 | } | |
800 | #endif | |
801 | // Local Variables: | |
802 | // mode: C++ | |
803 | // End: | |
804 | ||
805 | ||
806 |