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