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