J/Psi and Upsilon pt and y distributions from R. Vogt 2002 added.
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONlib.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.16  2002/10/14 14:55:35  hristov
19 Merging the VirtualMC branch to the main development branch (HEAD)
20
21 Revision 1.14.6.1  2002/06/10 14:57:41  hristov
22 Merged with v3-08-02
23
24 Revision 1.15  2002/04/17 10:11:51  morsch
25 Coding Rule violations corrected.
26
27 Revision 1.14  2002/02/22 17:26:43  morsch
28 Eta and omega added.
29
30 Revision 1.13  2001/03/27 11:01:04  morsch
31 Charm pt-distribution corrected. More realistic y-distribution for pi and K.
32
33 Revision 1.12  2001/03/09 13:01:41  morsch
34 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
35 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
36
37 Revision 1.11  2000/11/30 07:12:50  alibrary
38 Introducing new Rndm and QA classes
39
40 Revision 1.10  2000/06/29 21:08:27  morsch
41 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
42 This allows to pass a pointer to a library directly to AliGenParam and avoids the
43 use of function pointers in Config.C.
44
45 Revision 1.9  2000/06/14 15:20:56  morsch
46 Include clean-up (IH)
47
48 Revision 1.8  2000/06/09 20:32:11  morsch
49 All coding rule violations except RS3 corrected
50
51 Revision 1.7  2000/05/02 08:12:13  morsch
52 Coding rule violations corrected.
53
54 Revision 1.6  1999/09/29 09:24:14  fca
55 Introduction of the Copyright and cvs Log
56
57 */
58
59 // Library class for particle pt and y distributions used for 
60 // muon spectrometer simulations.
61 // To be used with AliGenParam.
62 // The following particle typed can be simulated:
63 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons. 
64 //
65 // andreas.morsch@cern.ch
66 //
67
68 #include "TMath.h"
69 #include "TRandom.h"
70
71 #include "AliGenMUONlib.h"
72
73 ClassImp(AliGenMUONlib)
74 //
75 //  Pions
76 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
77 {
78 //
79 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
80 //     POWER LAW FOR PT > 500 MEV
81 //     MT SCALING BELOW (T=160 MEV)
82 //
83   const Double_t kp0 = 1.3;
84   const Double_t kxn = 8.28;
85   const Double_t kxlim=0.5;
86   const Double_t kt=0.160;
87   const Double_t kxmpi=0.139;
88   const Double_t kb=1.;
89   Double_t y, y1, xmpi2, ynorm, a;
90   Double_t x=*px;
91   //
92   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
93   xmpi2=kxmpi*kxmpi;
94   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
95   a=ynorm/y1;
96   if (x > kxlim)
97     y=a*TMath::Power(kp0/(kp0+x),kxn);
98   else
99     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
100   return y*x;
101 }
102 //
103 // y-distribution
104 //
105 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
106 {
107 // Pion y
108   Double_t y=TMath::Abs(*py);
109 /*
110   const Double_t ka    = 7000.;
111   const Double_t kdy   = 4.;
112   Double_t ex = y*y/(2*kdy*kdy);
113   return ka*TMath::Exp(-ex);
114 */
115   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
116   
117 }
118 //                 particle composition
119 //
120 Int_t AliGenMUONlib::IpPion(TRandom *ran)
121 {
122 // Pion composition 
123     if (ran->Rndm() < 0.5) {
124         return  211;
125     } else {
126         return -211;
127     }
128 }
129
130 //____________________________________________________________
131 //
132 // Mt-scaling
133
134 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
135 {
136   //    SCALING EN MASSE PAR RAPPORT A PTPI
137   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
138   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
139   //     VALUE MESON/PI AT 5 GEV
140   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
141   np--;
142   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
143   Double_t fmax2=f5/kfmax[np];
144   // PIONS
145   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
146   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
147                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
148   return fmtscal*ptpion;
149 }
150 //
151 // kaon
152 //
153 //                pt-distribution
154 //____________________________________________________________
155 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
156 {
157 // Kaon pT
158   return PtScal(*px,2);
159 }
160
161 // y-distribution
162 //____________________________________________________________
163 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
164 {
165 // Kaon y
166   Double_t y=TMath::Abs(*py);
167 /*
168   const Double_t ka    = 1000.;
169   const Double_t kdy   = 4.;
170   //
171   Double_t ex = y*y/(2*kdy*kdy);
172   return ka*TMath::Exp(-ex);
173 */
174
175   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
176 }
177
178 //                 particle composition
179 //
180 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
181 {
182 // Kaon composition
183     if (ran->Rndm() < 0.5) {
184         return  321;
185     } else {
186         return -321;
187     }
188 }
189
190 //                    J/Psi 
191 //
192 //
193 //                pt-distribution
194 //____________________________________________________________
195 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
196 {
197 // J/Psi pT
198   const Double_t kpt0 = 4.;
199   const Double_t kxn  = 3.6;
200   Double_t x=*px;
201   //
202   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
203   return x/TMath::Power(pass1,kxn);
204 }
205
206 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
207 {
208 // J/Psi pT
209 //
210 // R. Vogt 2002
211 // PbPb 5.5 TeV
212 // MRST HO
213 // mc = 1.4 GeV, pt-kick 1 GeV
214 //
215
216     Float_t ptJpsi[100] = {
217         0.0000e-01,  4.5870e+01,  6.5200e+01,  7.1740e+01,  6.5090e+01,
218         5.5070e+01,  4.9420e+01,  3.9780e+01,  3.2390e+01,  2.8120e+01,
219         2.3870e+01,  1.9540e+01,  1.6510e+01,  1.4180e+01,  1.2050e+01,
220         1.0390e+01,  8.7970e+00,  7.8680e+00,  6.7710e+00,  5.9360e+00,
221         5.3460e+00,  4.5670e+00,  4.6500e+00,  3.9360e+00,  3.5070e+00,
222         3.2070e+00,  2.8310e+00,  2.6340e+00,  2.4900e+00,  2.2410e+00,
223         2.1090e+00,  1.9070e+00,  1.7360e+00,  1.6120e+00,  1.5450e+00,
224         1.4350e+00,  1.3890e+00,  1.2610e+00,  1.0880e+00,  1.0930e+00,
225         1.0680e+00,  9.2500e-01,  8.6790e-01,  8.1790e-01,  7.9770e-01,
226         7.4660e-01,  7.3110e-01,  6.5120e-01,  6.8140e-01,  5.7960e-01,
227         5.8210e-01,  5.4640e-01,  5.1700e-01,  5.0760e-01,  4.8280e-01,
228         4.5360e-01,  4.4910e-01,  4.2410e-01,  4.2100e-01,  3.9530e-01,
229         3.7220e-01,  3.4840e-01,  3.4550e-01,  3.3000e-01,  3.1670e-01,
230         3.1470e-01,  2.8920e-01,  2.7650e-01,  2.6860e-01,  2.5390e-01,
231         2.4190e-01,  2.5200e-01,  2.2960e-01,  2.2540e-01,  2.0950e-01,
232         2.0250e-01,  1.8720e-01,  1.8200e-01,  1.7860e-01,  1.8290e-01,
233         1.6970e-01,  1.7130e-01,  1.6310e-01,  1.5500e-01,  1.5100e-01,
234         1.5770e-01,  1.4240e-01,  1.4560e-01,  1.3330e-01,  1.4190e-01,
235         1.2010e-01,  1.2430e-01,  1.2430e-01,  1.1340e-01,  1.1840e-01,
236         1.1380e-01,  1.0330e-01,  1.0130e-01,  1.0390e-01,  9.5810e-02
237     };
238     Float_t x = px[0] * px[0];
239     if (x  < 1.5 || x > 100) {
240         return 0.0;
241     } else {
242         Float_t y =  Interpolate(x, ptJpsi, 0.5, 1., 100, 4);
243         return px[0] * y;
244     }
245 }
246 //
247 //               y-distribution
248 //____________________________________________________________
249 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
250 {
251 // J/psi y
252   const Double_t ky0 = 4.;
253   const Double_t kb=1.;
254   Double_t yj;
255   Double_t y=TMath::Abs(*py);
256   //
257   if (y < ky0)
258     yj=kb;
259   else
260     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
261   return yj;
262 }
263
264
265 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t *dummy)
266 {
267
268 //
269 // J/Psi y
270 //
271 //
272 // R. Vogt 2002
273 // PbPb 5.5 TeV
274 // MRST HO
275 // mc = 1.4 GeV, pt-kick 1 GeV
276 //
277
278     Float_t yJpsi[62] = 
279         {
280           0.3981E-03, 0.1169E-01, 0.6143E-01, 0.3554E+00, 0.1249E+01 
281         , 0.1677E+01, 0.3634E+01, 0.5414E+01, 0.8242E+01, 0.1102E+02
282         , 0.1353E+02, 0.1964E+02, 0.2357E+02, 0.2662E+02, 0.3023E+02 
283         , 0.3250E+02, 0.3137E+02, 0.3243E+02, 0.3120E+02, 0.3249E+02 
284         , 0.3166E+02, 0.3104E+02, 0.3203E+02, 0.3149E+02, 0.3117E+02 
285         , 0.3210E+02, 0.3170E+02, 0.3279E+02, 0.3079E+02, 0.3208E+02 
286         , 0.3218E+02, 0.3218E+02, 0.3208E+02, 0.3079E+02, 0.3279E+02 
287         , 0.3170E+02, 0.3210E+02, 0.3118E+02, 0.3149E+02, 0.3203E+02 
288         , 0.3104E+02, 0.3167E+02, 0.3250E+02, 0.3120E+02, 0.3243E+02 
289         , 0.3137E+02, 0.3250E+02, 0.3022E+02, 0.2662E+02, 0.2357E+02 
290         , 0.1964E+02, 0.1353E+02, 0.1102E+02, 0.8242E+01, 0.5414E+01 
291         , 0.3634E+01, 0.1677E+01, 0.1249E+01, 0.3554E+00, 0.6142E-01
292         , 0.1169E-01, 0.3981E-03};
293
294     return Interpolate(px[0], yJpsi, -7.625, 0.25, 62, 2);
295 }
296
297 //                 particle composition
298 //
299 Int_t AliGenMUONlib::IpJpsi(TRandom *)
300 {
301 // J/Psi composition
302     return 443;
303 }
304
305 //                      Upsilon
306 //
307 //
308 //                  pt-distribution
309 //____________________________________________________________
310 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
311 {
312 // Upsilon pT
313   const Double_t kpt0 = 5.3;
314   const Double_t kxn  = 2.5;
315   Double_t x=*px;
316   //
317   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
318   return x/TMath::Power(pass1,kxn);
319 }
320
321 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t *dummy)
322 {
323
324 //
325 // Upsilon pT
326 //
327 //
328 // R. Vogt 2002
329 // PbPb 5.5 TeV
330 // MRST HO
331 // mc = 1.4 GeV, pt-kick 1 GeV
332 //
333
334     Float_t ptUps[100] = 
335         {   
336             0.0000e-01, -1.5290e-02,  2.6020e-01,  3.4220e-01,  3.3710e-01,
337             3.1880e-01,  3.3420e-01,  2.7740e-01,  2.3730e-01,  2.0640e-01,
338             1.7690e-01,  1.6190e-01,  1.4500e-01,  1.3310e-01,  1.1440e-01,
339             1.0800e-01,  1.0210e-01,  8.4690e-02,  8.0050e-02,  7.0710e-02,
340             6.4160e-02,  6.5200e-02,  6.6890e-02,  6.0600e-02,  5.4030e-02,
341             5.1140e-02,  4.6120e-02,  4.4800e-02,  4.2490e-02,  4.1440e-02,
342             4.0310e-02,  3.7110e-02,  3.5890e-02,  3.5420e-02,  3.0370e-02,
343             2.9970e-02,  3.0770e-02,  2.6380e-02,  2.7740e-02,  2.6690e-02,
344             2.4210e-02,  2.5200e-02,  2.3760e-02,  2.1370e-02,  2.2290e-02,
345             2.2700e-02,  2.0110e-02,  1.9320e-02,  1.8830e-02,  1.9910e-02,
346             1.9740e-02,  1.8460e-02,  1.8240e-02,  1.6740e-02,  1.6140e-02,
347             1.7340e-02,  1.5950e-02,  1.5430e-02,  1.4780e-02,  1.2750e-02,
348             1.4370e-02,  1.2810e-02,  1.2900e-02,  1.1070e-02,  1.1830e-02,
349             1.1150e-02,  1.1260e-02,  1.1610e-02,  1.0700e-02,  1.1600e-02,
350             1.0390e-02,  1.0280e-02,  1.0180e-02,  1.0030e-02,  9.6050e-03,
351             8.8050e-03,  8.9680e-03,  9.0120e-03,  8.4110e-03,  8.6660e-03,
352             8.3060e-03,  8.5850e-03,  8.2600e-03,  8.3800e-03,  8.4200e-03,
353             7.5690e-03,  7.2100e-03,  7.1230e-03,  7.3350e-03,  7.1980e-03,
354             6.7500e-03,  6.6190e-03,  6.3370e-03,  6.6270e-03,  6.8290e-03,
355             6.0880e-03,  6.6310e-03,  6.0490e-03,  5.8900e-03,  5.6100e-03
356         };
357     Float_t x = px[0] * px[0];
358     if (x  < 1.5 || x > 100) {
359         return 0.0;
360     } else {
361         Float_t y =  Interpolate(x, ptUps, 0.5, 1., 100, 4);
362         return px[0] * y;
363     }
364 }
365
366 //
367 //                    y-distribution
368 //
369 //____________________________________________________________
370 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
371 {
372 // Upsilon y
373   const Double_t ky0 = 3.;
374   const Double_t kb=1.;
375   Double_t yu;
376   Double_t y=TMath::Abs(*py);
377   //
378   if (y < ky0)
379     yu=kb;
380   else
381     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
382   return yu;
383 }
384
385
386 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t *dummy)
387 {
388
389 //
390 // Upsilon y
391 //
392 //
393 // R. Vogt 2002
394 // PbPb 5.5 TeV
395 // MRST HO
396 // mc = 1.4 GeV, pt-kick 1 GeV
397 //
398
399     Float_t yUps[52] = 
400         {
401             0.000000, 0.000065, 0.000767, 0.004332, 0.014790,
402             0.029370, 0.052060, 0.077930, 0.122500, 0.135800,
403             0.184000, 0.207900, 0.228300, 0.258500, 0.269500,
404             0.288500, 0.316200, 0.304100, 0.315800, 0.323300,
405             0.322400, 0.322600, 0.345500, 0.338100, 0.331900,
406             0.343700, 0.343700, 0.331900, 0.338100, 0.345500,
407             0.322600, 0.322400, 0.323300, 0.315800, 0.304100,
408             0.316200, 0.288500, 0.269500, 0.258500, 0.228300,
409             0.207900, 0.184000, 0.135800, 0.122500, 0.077930,
410             0.052060, 0.029380, 0.014780, 0.004332, 0.000767,
411             0.6479E-04, 0.1013E-06       
412         };
413     
414
415
416     return Interpolate(px[0], yUps, -6.5, 0.25, 52, 2);
417 }
418
419 //                 particle composition
420 //
421 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
422 {
423 // y composition
424     return 553;
425 }
426
427 //
428 //                        Phi
429 //
430 //
431 //    pt-distribution (by scaling of pion distribution)
432 //____________________________________________________________
433 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
434 {
435 // Phi pT
436   return PtScal(*px,7);
437 }
438 //    y-distribution
439 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
440 {
441 // Phi y
442     Double_t *dum=0;
443     return YJpsi(px,dum);
444 }
445 //                 particle composition
446 //
447 Int_t AliGenMUONlib::IpPhi(TRandom *)
448 {
449 // Phi composition
450     return 333;
451 }
452
453 //
454 //                        omega
455 //
456 //
457 //    pt-distribution (by scaling of pion distribution)
458 //____________________________________________________________
459 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
460 {
461 // Omega pT
462   return PtScal(*px,5);
463 }
464 //    y-distribution
465 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
466 {
467 // Omega y
468     Double_t *dum=0;
469     return YJpsi(px,dum);
470 }
471 //                 particle composition
472 //
473 Int_t AliGenMUONlib::IpOmega(TRandom *)
474 {
475 // Omega composition
476     return 223;
477 }
478
479
480 //
481 //                        Eta
482 //
483 //
484 //    pt-distribution (by scaling of pion distribution)
485 //____________________________________________________________
486 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
487 {
488 // Eta pT
489   return PtScal(*px,3);
490 }
491 //    y-distribution
492 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
493 {
494 // Eta y
495     Double_t *dum=0;
496     return YJpsi(px,dum);
497 }
498 //                 particle composition
499 //
500 Int_t AliGenMUONlib::IpEta(TRandom *)
501 {
502 // Eta composition
503     return 221;
504 }
505
506 //
507 //                        Charm
508 //
509 //
510 //                    pt-distribution
511 //____________________________________________________________
512 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
513 {
514 // Charm pT
515   const Double_t kpt0 = 4.08;
516   const Double_t kxn  = 9.40;
517
518   Double_t x=*px;
519   //
520   Double_t pass1 = 1.+(x/kpt0);
521   return x/TMath::Power(pass1,kxn);
522 }
523 //                  y-distribution
524 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
525 {
526 // Charm y
527     Double_t *dum=0;
528     return YJpsi(px,dum);
529 }
530
531 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
532 {  
533 // Charm composition
534     Float_t random;
535     Int_t ip;
536 //    411,421,431,4122
537     random = ran->Rndm();
538     if (random < 0.5) {
539         ip=411;
540     } else if (random < 0.75) {
541         ip=421;
542     } else if (random < 0.90) {
543         ip=431;
544     } else {
545         ip=4122;
546     }
547     if (ran->Rndm() < 0.5) {ip=-ip;}
548     
549     return ip;
550 }
551
552
553 //
554 //                        Beauty
555 //
556 //
557 //                    pt-distribution
558 //____________________________________________________________
559 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
560 {
561 // Beauty pT
562   const Double_t kpt0 = 4.;
563   const Double_t kxn  = 3.6;
564   Double_t x=*px;
565   //
566   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
567   return x/TMath::Power(pass1,kxn);
568 }
569 //                     y-distribution
570 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
571 {
572 // Beauty y
573     Double_t *dum=0;
574     return YJpsi(px,dum);
575 }
576
577 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
578 {  
579 // Beauty Composition
580     Float_t random;
581     Int_t ip;
582     random = ran->Rndm();
583     if (random < 0.5) {
584         ip=511;
585     } else if (random < 0.75) {
586         ip=521;
587     } else if (random < 0.90) {
588         ip=531;
589     } else {
590         ip=5122;
591     }
592     if (ran->Rndm() < 0.5) {ip=-ip;}
593     
594     return ip;
595 }
596
597 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
598 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
599 {
600 // Return pointer to pT parameterisation
601     TString sname = TString(tname);
602     GenFunc func;
603     switch (param) 
604     {
605     case kPhi:
606         func=PtPhi;
607         break;
608     case kOmega:
609         func=PtOmega;
610         break;
611     case kEta:
612         func=PtEta;
613         break;
614     case kJpsi:
615         if (sname == "PbPb") {
616             func=PtJpsiPbPb;
617         } else {
618             func=PtJpsi;
619         }
620         break;
621     case kUpsilon:
622         if (sname == "PbPb") {
623             func=PtUpsilonPbPb;
624         } else {
625             func=PtUpsilon;
626         }
627         break;
628     case kCharm:
629         func=PtCharm;
630         break;
631     case kBeauty:
632         func=PtBeauty;
633         break;
634     case kPion:
635         func=PtPion;
636         break;
637     case kKaon:
638         func=PtKaon;
639         break;
640     default:
641         func=0;
642         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
643     }
644     return func;
645 }
646
647 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
648 {
649     TString sname = TString(tname);
650     
651 // Return pointer to y- parameterisation
652     GenFunc func;
653     switch (param) 
654     {
655     case kPhi:
656         func=YPhi;
657         break;
658     case kEta:
659         func=YEta;
660         break;
661     case kOmega:
662         func=YOmega;
663         break;
664     case kJpsi:
665         if (sname == "PbPb") {
666             func=YJpsiPbPb;
667         } else {
668             func=YJpsi;
669         }
670         break;
671     case kUpsilon:
672         if (sname == "PbPb") {
673             func=YUpsilonPbPb;
674         } else {
675             func=YUpsilon;
676         }
677         break;
678     case kCharm:
679         func=YCharm;
680         break;
681     case kBeauty:
682         func=YBeauty;
683         break;
684     case kPion:
685         func=YPion;
686         break;
687     case kKaon:
688         func=YKaon;
689         break;
690     default:
691         func=0;
692         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
693     }
694     return func;
695 }
696 typedef Int_t (*GenFuncIp) (TRandom *);
697 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
698 {
699 // Return pointer to particle type parameterisation
700     GenFuncIp func;
701     switch (param) 
702     {
703     case kPhi:
704         func=IpPhi;
705         break;
706     case kEta:
707         func=IpEta;
708         break;
709     case kOmega:
710         func=IpOmega;
711         break;
712     case kJpsi:
713         func=IpJpsi;
714         break;
715     case kUpsilon:
716         func=IpUpsilon;
717         break;
718     case kCharm:
719         func=IpCharm;
720         break;
721     case kBeauty:
722         func=IpBeauty;
723         break;
724     case kPion:
725         func=IpPion;
726         break;
727     case kKaon:
728         func=IpKaon;
729         break;
730     default:
731         func=0;
732         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
733     }
734     return func;
735 }
736
737
738
739 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0, 
740                                    Float_t dx,
741                                    Int_t n, Int_t no)
742 {
743 //
744 // Neville's alorithm for interpolation
745 //
746 // x:  x-value
747 // y:  Input array
748 // x0: minimum x 
749 // dx: step size
750 //  n: number of data points
751 // no: order of polynom 
752 //
753     Float_t*  c = new Float_t[n];
754     Float_t*  d = new Float_t[n];
755     Int_t m, i;
756     for (i = 0; i < n; i++) {
757         c[i] = y[i];
758         d[i] = y[i];
759     }
760     
761     Int_t   ns  = int((x - x0)/dx);
762     
763     Float_t y1  = y[ns];
764     ns--;    
765     for (m = 0; m < no; m++) {
766         for (i = 0; i < n-m; i++) {     
767             Float_t ho = x0 + Float_t(i) * dx - x;
768             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
769             Float_t w  = c[i+1] - d[i];
770             Float_t den = ho-hp;
771             den = w/den;
772             d[i] = hp * den;
773             c[i] = ho * den;
774         }
775         Float_t dy;
776         
777         if (2*ns < (n-m-1)) {
778             dy  = c[ns+1];
779         } else {
780             dy  = d[ns--];
781         }
782         y1 += dy;}
783     delete[] c;
784     delete[] d;
785
786     return y1;
787 }
788
789