Psi and Upsilon family generated.
[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.21  2003/03/25 12:01:11  morsch
19 Quarkonia for pp @ 14 TeV added.
20
21 Revision 1.20  2003/03/13 11:54:39  morsch
22 Limited pT range for parameterized Upsilon and J/Psi pT distributions.
23
24 Revision 1.19  2003/02/24 16:46:11  morsch
25 New parameterisation for Psi and Upsilon (PbPb)
26
27 Revision 1.18  2002/11/07 09:13:42  morsch
28 Use "Vogt" to label new distributions.
29
30 Revision 1.17  2002/11/07 09:06:10  morsch
31 J/Psi and Upsilon pt and y distributions from R. Vogt 2002 added.
32
33 Revision 1.16  2002/10/14 14:55:35  hristov
34 Merging the VirtualMC branch to the main development branch (HEAD)
35
36 Revision 1.14.6.1  2002/06/10 14:57:41  hristov
37 Merged with v3-08-02
38
39 Revision 1.15  2002/04/17 10:11:51  morsch
40 Coding Rule violations corrected.
41
42 Revision 1.14  2002/02/22 17:26:43  morsch
43 Eta and omega added.
44
45 Revision 1.13  2001/03/27 11:01:04  morsch
46 Charm pt-distribution corrected. More realistic y-distribution for pi and K.
47
48 Revision 1.12  2001/03/09 13:01:41  morsch
49 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
50 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
51
52 Revision 1.11  2000/11/30 07:12:50  alibrary
53 Introducing new Rndm and QA classes
54
55 Revision 1.10  2000/06/29 21:08:27  morsch
56 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
57 This allows to pass a pointer to a library directly to AliGenParam and avoids the
58 use of function pointers in Config.C.
59
60 Revision 1.9  2000/06/14 15:20:56  morsch
61 Include clean-up (IH)
62
63 Revision 1.8  2000/06/09 20:32:11  morsch
64 All coding rule violations except RS3 corrected
65
66 Revision 1.7  2000/05/02 08:12:13  morsch
67 Coding rule violations corrected.
68
69 Revision 1.6  1999/09/29 09:24:14  fca
70 Introduction of the Copyright and cvs Log
71
72 */
73
74 // Library class for particle pt and y distributions used for 
75 // muon spectrometer simulations.
76 // To be used with AliGenParam.
77 // The following particle typed can be simulated:
78 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons. 
79 //
80 // andreas.morsch@cern.ch
81 //
82
83 #include "TMath.h"
84 #include "TRandom.h"
85
86 #include "AliGenMUONlib.h"
87
88 ClassImp(AliGenMUONlib)
89 //
90 //  Pions
91 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
92 {
93 //
94 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
95 //     POWER LAW FOR PT > 500 MEV
96 //     MT SCALING BELOW (T=160 MEV)
97 //
98   const Double_t kp0 = 1.3;
99   const Double_t kxn = 8.28;
100   const Double_t kxlim=0.5;
101   const Double_t kt=0.160;
102   const Double_t kxmpi=0.139;
103   const Double_t kb=1.;
104   Double_t y, y1, xmpi2, ynorm, a;
105   Double_t x=*px;
106   //
107   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
108   xmpi2=kxmpi*kxmpi;
109   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
110   a=ynorm/y1;
111   if (x > kxlim)
112     y=a*TMath::Power(kp0/(kp0+x),kxn);
113   else
114     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
115   return y*x;
116 }
117 //
118 // y-distribution
119 //
120 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
121 {
122 // Pion y
123   Double_t y=TMath::Abs(*py);
124 /*
125   const Double_t ka    = 7000.;
126   const Double_t kdy   = 4.;
127   Double_t ex = y*y/(2*kdy*kdy);
128   return ka*TMath::Exp(-ex);
129 */
130   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
131   
132 }
133 //                 particle composition
134 //
135 Int_t AliGenMUONlib::IpPion(TRandom *ran)
136 {
137 // Pion composition 
138     if (ran->Rndm() < 0.5) {
139         return  211;
140     } else {
141         return -211;
142     }
143 }
144
145 //____________________________________________________________
146 //
147 // Mt-scaling
148
149 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
150 {
151   //    SCALING EN MASSE PAR RAPPORT A PTPI
152   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
153   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
154   //     VALUE MESON/PI AT 5 GEV
155   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
156   np--;
157   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
158   Double_t fmax2=f5/kfmax[np];
159   // PIONS
160   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
161   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
162                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
163   return fmtscal*ptpion;
164 }
165 //
166 // kaon
167 //
168 //                pt-distribution
169 //____________________________________________________________
170 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
171 {
172 // Kaon pT
173   return PtScal(*px,2);
174 }
175
176 // y-distribution
177 //____________________________________________________________
178 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
179 {
180 // Kaon y
181   Double_t y=TMath::Abs(*py);
182 /*
183   const Double_t ka    = 1000.;
184   const Double_t kdy   = 4.;
185   //
186   Double_t ex = y*y/(2*kdy*kdy);
187   return ka*TMath::Exp(-ex);
188 */
189
190   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
191 }
192
193 //                 particle composition
194 //
195 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
196 {
197 // Kaon composition
198     if (ran->Rndm() < 0.5) {
199         return  321;
200     } else {
201         return -321;
202     }
203 }
204
205 //                    J/Psi 
206 //
207 //
208 //                pt-distribution
209 //____________________________________________________________
210 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
211 {
212 // J/Psi pT
213   const Double_t kpt0 = 4.;
214   const Double_t kxn  = 3.6;
215   Double_t x=*px;
216   //
217   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
218   return x/TMath::Power(pass1,kxn);
219 }
220
221 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
222 {
223 // J/Psi pT spectrum
224 //
225 // R. Vogt 2002
226 // PbPb 5.5 TeV
227 // MRST HO
228 // mc = 1.4 GeV, pt-kick 1 GeV
229 //
230     Float_t x = px[0];
231     Float_t c[8] = {
232         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00, 
233         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
234     };
235     
236     Double_t y;
237     if (x < 10.) {
238         Int_t j;
239         y = c[j = 7];
240         while (j > 0) y  = y * x +c[--j];
241         y = x * TMath::Exp(y);
242     } else {
243         y = 0.;
244     }
245     return y;
246 }
247 Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t *dummy)
248 {
249 // J/Psi pT spectrum
250 //
251 // R. Vogt 2002
252 // pp 14 TeV
253 // MRST HO
254 // mc = 1.4 GeV, pt-kick 1 GeV
255 //
256     Float_t x = px[0];
257     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
258  
259     Double_t y;
260     if (x < 10.) {
261         Int_t j;
262         y = c[j = 3];
263         while (j > 0) y  = y * x +c[--j];
264         y = x * TMath::Exp(y);
265     } else {
266         y = 0.;
267     }
268     return y;
269 }
270
271 //
272 //               y-distribution
273 //____________________________________________________________
274 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
275 {
276 // J/psi y
277   const Double_t ky0 = 4.;
278   const Double_t kb=1.;
279   Double_t yj;
280   Double_t y=TMath::Abs(*py);
281   //
282   if (y < ky0)
283     yj=kb;
284   else
285     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
286   return yj;
287 }
288
289
290 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t *dummy)
291 {
292
293 //
294 // J/Psi y
295 //
296 //
297 // R. Vogt 2002
298 // PbPb 5.5 TeV
299 // MRST HO
300 // mc = 1.4 GeV, pt-kick 1 GeV
301 //
302     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
303     Double_t x = TMath::Abs(px[0]);
304     Double_t y;
305     
306     if (x < 4.) {
307         y = 31.754;
308     } else if (x < 6) {
309         Int_t j;
310         y = c[j = 4];
311         while (j > 0) y  = y * x + c[--j];
312     } else {
313         y =0.;
314     }
315     
316     return y;
317 }
318
319 Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t *dummy)
320 {
321
322 //
323 // J/Psi y
324 //
325 //
326 // R. Vogt 2002
327 // pp 14  TeV
328 // MRST HO
329 // mc = 1.4 GeV, pt-kick 1 GeV
330 //
331
332     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
333     Double_t x = TMath::Abs(px[0]);
334     Double_t y;
335     
336     if (x < 2.5) {
337         y = 96.455 - 0.8483 * x * x;
338     } else if (x < 7.9) {
339         Int_t j;
340         y = c[j = 4];
341         while (j > 0) y  = y * x + c[--j];
342     } else {
343         y =0.;
344     }
345     
346     return y;
347 }
348
349 //                 particle composition
350 //
351 Int_t AliGenMUONlib::IpJpsi(TRandom *)
352 {
353 // J/Psi composition
354     Int_t ip;
355     Float_t r = gRandom->Rndm();
356     if (r < 0.98) {
357         ip = 443;
358     } else {
359         ip = 100443;
360     }
361     return ip;
362 }
363
364 //                      Upsilon
365 //
366 //
367 //                  pt-distribution
368 //____________________________________________________________
369 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
370 {
371 // Upsilon pT
372   const Double_t kpt0 = 5.3;
373   const Double_t kxn  = 2.5;
374   Double_t x=*px;
375   //
376   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
377   return x/TMath::Power(pass1,kxn);
378 }
379
380 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t *dummy)
381 {
382
383 //
384 // Upsilon pT
385 //
386 //
387 // R. Vogt 2002
388 // PbPb 5.5 TeV
389 // MRST HO
390 // mc = 1.4 GeV, pt-kick 1 GeV
391 //
392     Float_t x = px[0];
393     Double_t c[8] = {
394         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,       
395         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
396     };
397     Double_t y;
398     if (x < 10.) {
399         Int_t j;
400         y = c[j = 7];
401         while (j > 0) y  = y * x +c[--j];
402         y = x * TMath::Exp(y);
403     } else {
404         y = 0.;
405     }
406     return y;
407 }
408
409 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t *dummy)
410 {
411
412 //
413 // Upsilon pT
414 //
415 //
416 // R. Vogt 2002
417 // pp 14 TeV
418 // MRST HO
419 // mc = 1.4 GeV, pt-kick 1 GeV
420 //
421     Float_t x = px[0];
422     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,   
423                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
424     
425     Double_t y;
426     if (x < 10.) {
427         Int_t j;
428         y = c[j = 7];
429         while (j > 0) y  = y * x +c[--j];
430         y = x * TMath::Exp(y);
431     } else {
432         y = 0.;
433     }
434     return y;
435 }
436
437 //
438 //                    y-distribution
439 //
440 //____________________________________________________________
441 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
442 {
443 // Upsilon y
444   const Double_t ky0 = 3.;
445   const Double_t kb=1.;
446   Double_t yu;
447   Double_t y=TMath::Abs(*py);
448   //
449   if (y < ky0)
450     yu=kb;
451   else
452     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
453   return yu;
454 }
455
456
457 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t *dummy)
458 {
459
460 //
461 // Upsilon y
462 //
463 //
464 // R. Vogt 2002
465 // PbPb 5.5 TeV
466 // MRST HO
467 // mc = 1.4 GeV, pt-kick 1 GeV
468 //
469
470     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
471                      -2.99753e-09, 1.28895e-05};
472         
473     Double_t x = px[0];
474     if (TMath::Abs(x) > 5.55) return 0.;
475     Int_t j;
476     Double_t y = c[j = 6];
477     while (j > 0) y  = y * x +c[--j];
478     return y;
479 }
480
481 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t *dummy)
482 {
483
484 //
485 // Upsilon y
486 //
487 //
488 // R. Vogt 2002
489 // p p  14. TeV
490 // MRST HO
491 // mc = 1.4 GeV, pt-kick 1 GeV
492 //
493     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04, 
494                      -6.20539e-10, 1.29943e-05};
495                 
496     Double_t x = px[0];
497     if (TMath::Abs(x) > 6.2) return 0.;
498     Int_t j;
499     Double_t y = c[j = 6];
500     while (j > 0) y  = y * x +c[--j];
501     return y;
502 }
503
504 //                 particle composition
505 //
506 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
507 {
508 // y composition
509     Int_t ip;
510     Float_t r = gRandom->Rndm();
511     
512     if (r < 0.712) {
513         ip = 553;
514     } else if (r < 0.896) {
515         ip = 100553;
516     } else {
517         ip = 200553;
518     }
519
520     return ip;
521 }
522
523 //
524 //                        Phi
525 //
526 //
527 //    pt-distribution (by scaling of pion distribution)
528 //____________________________________________________________
529 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
530 {
531 // Phi pT
532   return PtScal(*px,7);
533 }
534 //    y-distribution
535 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
536 {
537 // Phi y
538     Double_t *dum=0;
539     return YJpsi(px,dum);
540 }
541 //                 particle composition
542 //
543 Int_t AliGenMUONlib::IpPhi(TRandom *)
544 {
545 // Phi composition
546     return 333;
547 }
548
549 //
550 //                        omega
551 //
552 //
553 //    pt-distribution (by scaling of pion distribution)
554 //____________________________________________________________
555 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
556 {
557 // Omega pT
558   return PtScal(*px,5);
559 }
560 //    y-distribution
561 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
562 {
563 // Omega y
564     Double_t *dum=0;
565     return YJpsi(px,dum);
566 }
567 //                 particle composition
568 //
569 Int_t AliGenMUONlib::IpOmega(TRandom *)
570 {
571 // Omega composition
572     return 223;
573 }
574
575
576 //
577 //                        Eta
578 //
579 //
580 //    pt-distribution (by scaling of pion distribution)
581 //____________________________________________________________
582 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
583 {
584 // Eta pT
585   return PtScal(*px,3);
586 }
587 //    y-distribution
588 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
589 {
590 // Eta y
591     Double_t *dum=0;
592     return YJpsi(px,dum);
593 }
594 //                 particle composition
595 //
596 Int_t AliGenMUONlib::IpEta(TRandom *)
597 {
598 // Eta composition
599     return 221;
600 }
601
602 //
603 //                        Charm
604 //
605 //
606 //                    pt-distribution
607 //____________________________________________________________
608 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
609 {
610 // Charm pT
611   const Double_t kpt0 = 4.08;
612   const Double_t kxn  = 9.40;
613
614   Double_t x=*px;
615   //
616   Double_t pass1 = 1.+(x/kpt0);
617   return x/TMath::Power(pass1,kxn);
618 }
619 //                  y-distribution
620 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
621 {
622 // Charm y
623     Double_t *dum=0;
624     return YJpsi(px,dum);
625 }
626
627 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
628 {  
629 // Charm composition
630     Float_t random;
631     Int_t ip;
632 //    411,421,431,4122
633     random = ran->Rndm();
634     if (random < 0.5) {
635         ip=411;
636     } else if (random < 0.75) {
637         ip=421;
638     } else if (random < 0.90) {
639         ip=431;
640     } else {
641         ip=4122;
642     }
643     if (ran->Rndm() < 0.5) {ip=-ip;}
644     
645     return ip;
646 }
647
648
649 //
650 //                        Beauty
651 //
652 //
653 //                    pt-distribution
654 //____________________________________________________________
655 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
656 {
657 // Beauty pT
658   const Double_t kpt0 = 4.;
659   const Double_t kxn  = 3.6;
660   Double_t x=*px;
661   //
662   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
663   return x/TMath::Power(pass1,kxn);
664 }
665 //                     y-distribution
666 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
667 {
668 // Beauty y
669     Double_t *dum=0;
670     return YJpsi(px,dum);
671 }
672
673 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
674 {  
675 // Beauty Composition
676     Float_t random;
677     Int_t ip;
678     random = ran->Rndm();
679     if (random < 0.5) {
680         ip=511;
681     } else if (random < 0.75) {
682         ip=521;
683     } else if (random < 0.90) {
684         ip=531;
685     } else {
686         ip=5122;
687     }
688     if (ran->Rndm() < 0.5) {ip=-ip;}
689     
690     return ip;
691 }
692
693 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
694 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
695 {
696 // Return pointer to pT parameterisation
697     TString sname = TString(tname);
698     GenFunc func;
699     switch (param) 
700     {
701     case kPhi:
702         func=PtPhi;
703         break;
704     case kOmega:
705         func=PtOmega;
706         break;
707     case kEta:
708         func=PtEta;
709         break;
710     case kJpsi:
711         if (sname == "Vogt" || sname == "Vogt PbPb") {
712             func=PtJpsiPbPb;
713         } else if (sname == "Vogt pp") {
714             func=PtJpsiPP;
715         } else {
716             func=PtJpsi;
717         }
718         break;
719     case kUpsilon:
720         if (sname == "Vogt" || sname == "Vogt PbPb") {
721             func=PtUpsilonPbPb;
722         } else if (sname == "Vogt pp") {
723             func=PtUpsilonPP;
724         } else {
725             func=PtUpsilon;
726         }
727         break;
728     case kCharm:
729         func=PtCharm;
730         break;
731     case kBeauty:
732         func=PtBeauty;
733         break;
734     case kPion:
735         func=PtPion;
736         break;
737     case kKaon:
738         func=PtKaon;
739         break;
740     default:
741         func=0;
742         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
743     }
744     return func;
745 }
746
747 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
748 {
749     TString sname = TString(tname);
750     
751 // Return pointer to y- parameterisation
752     GenFunc func;
753     switch (param) 
754     {
755     case kPhi:
756         func=YPhi;
757         break;
758     case kEta:
759         func=YEta;
760         break;
761     case kOmega:
762         func=YOmega;
763         break;
764     case kJpsi:
765         if (sname == "Vogt" || sname == "Vogt PbPb") {
766             func=YJpsiPbPb;
767         } else if (sname == "Vogt pp"){
768             func=YJpsiPP;
769         } else {
770             func=YJpsi;
771         }
772         
773         break;
774     case kUpsilon:
775         if (sname == "Vogt" || sname == "Vogt PbPb") {
776             func=YUpsilonPbPb;
777         } else if (sname == "Vogt pp") {
778             func = YUpsilonPP;
779         } else {
780             func=YUpsilon;
781         }
782         break;
783     case kCharm:
784         func=YCharm;
785         break;
786     case kBeauty:
787         func=YBeauty;
788         break;
789     case kPion:
790         func=YPion;
791         break;
792     case kKaon:
793         func=YKaon;
794         break;
795     default:
796         func=0;
797         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
798     }
799     return func;
800 }
801 typedef Int_t (*GenFuncIp) (TRandom *);
802 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
803 {
804 // Return pointer to particle type parameterisation
805     GenFuncIp func;
806     switch (param) 
807     {
808     case kPhi:
809         func=IpPhi;
810         break;
811     case kEta:
812         func=IpEta;
813         break;
814     case kOmega:
815         func=IpOmega;
816         break;
817     case kJpsi:
818         func=IpJpsi;
819         break;
820     case kUpsilon:
821         func=IpUpsilon;
822         break;
823     case kCharm:
824         func=IpCharm;
825         break;
826     case kBeauty:
827         func=IpBeauty;
828         break;
829     case kPion:
830         func=IpPion;
831         break;
832     case kKaon:
833         func=IpKaon;
834         break;
835     default:
836         func=0;
837         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
838     }
839     return func;
840 }
841
842
843
844 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0, 
845                                    Float_t dx,
846                                    Int_t n, Int_t no)
847 {
848 //
849 // Neville's alorithm for interpolation
850 //
851 // x:  x-value
852 // y:  Input array
853 // x0: minimum x 
854 // dx: step size
855 //  n: number of data points
856 // no: order of polynom 
857 //
858     Float_t*  c = new Float_t[n];
859     Float_t*  d = new Float_t[n];
860     Int_t m, i;
861     for (i = 0; i < n; i++) {
862         c[i] = y[i];
863         d[i] = y[i];
864     }
865     
866     Int_t   ns  = int((x - x0)/dx);
867     
868     Float_t y1  = y[ns];
869     ns--;    
870     for (m = 0; m < no; m++) {
871         for (i = 0; i < n-m; i++) {     
872             Float_t ho = x0 + Float_t(i) * dx - x;
873             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
874             Float_t w  = c[i+1] - d[i];
875             Float_t den = ho-hp;
876             den = w/den;
877             d[i] = hp * den;
878             c[i] = ho * den;
879         }
880         Float_t dy;
881         
882         if (2*ns < (n-m-1)) {
883             dy  = c[ns+1];
884         } else {
885             dy  = d[ns--];
886         }
887         y1 += dy;}
888     delete[] c;
889     delete[] d;
890
891     return y1;
892 }
893
894