Quarkonia for pp add. (D. Stocco)
[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::PtJpsiCDFscaled( Double_t *px, Double_t */*dummy*/)
166 {
167 // J/Psi pT
168   const Double_t kpt0 = 4.703;
169   const Double_t kxn  = 3.826;
170   Double_t x=*px;
171   //
172   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
173   return x/TMath::Power(pass1,kxn);
174 }
175
176 Double_t AliGenMUONlib::PtJpsiCDFscaledPP( Double_t *px, Double_t */*dummy*/)
177 {
178 // J/Psi pT
179 //
180 // pp 14 TeV
181 //
182 // scaled from CDF data at 2 TeV
183
184   const Double_t kpt0 = 5.355;
185   const Double_t kxn  = 3.821;
186   Double_t x=*px;
187   //
188   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
189   return x/TMath::Power(pass1,kxn);
190 }
191
192 Double_t AliGenMUONlib::PtJpsiFlat( Double_t */*px*/, Double_t */*dummy*/ )
193 {
194   return 1.;
195 }
196
197 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t */*dummy*/)
198 {
199 // J/Psi pT spectrum
200 //
201 // R. Vogt 2002
202 // PbPb 5.5 TeV
203 // MRST HO
204 // mc = 1.4 GeV, pt-kick 1 GeV
205 //
206     Float_t x = px[0];
207     Float_t c[8] = {
208         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00, 
209         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
210     };
211     
212     Double_t y;
213     if (x < 10.) {
214         Int_t j;
215         y = c[j = 7];
216         while (j > 0) y  = y * x +c[--j];
217         y = x * TMath::Exp(y);
218     } else {
219         y = 0.;
220     }
221     return y;
222 }
223
224 Double_t AliGenMUONlib::PtJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
225 {
226 // J/Psi pT spectrum
227 // B -> J/Psi X
228     Double_t x0 =   4.0384;
229     Double_t  n =   3.0288;
230     
231     Double_t x = px[0];
232     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
233     
234     return y;
235 }
236
237
238 Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t */*dummy*/)
239 {
240 // J/Psi pT spectrum
241 //
242 // R. Vogt 2002
243 // pp 14 TeV
244 // MRST HO
245 // mc = 1.4 GeV, pt-kick 1 GeV
246 //
247     Float_t x = px[0];
248     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
249  
250     Double_t y;
251     if (x < 10.) {
252         Int_t j;
253         y = c[j = 3];
254         while (j > 0) y  = y * x +c[--j];
255         y = x * TMath::Exp(y);
256     } else {
257         y = 0.;
258     }
259     return y;
260 }
261
262 //
263 //               y-distribution
264 //____________________________________________________________
265 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t */*dummy*/)
266 {
267 // J/psi y
268   const Double_t ky0 = 4.;
269   const Double_t kb=1.;
270   Double_t yj;
271   Double_t y=TMath::Abs(*py);
272   //
273   if (y < ky0)
274     yj=kb;
275   else
276     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
277   return yj;
278 }
279
280 Double_t AliGenMUONlib::YJpsiFlat( Double_t */*py*/, Double_t */*dummy*/ )
281 {
282   return 1.;
283 }
284
285
286 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t */*dummy*/)
287 {
288
289 //
290 // J/Psi y
291 //
292 //
293 // R. Vogt 2002
294 // PbPb 5.5 TeV
295 // MRST HO
296 // mc = 1.4 GeV, pt-kick 1 GeV
297 //
298     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
299     Double_t x = TMath::Abs(px[0]);
300     Double_t y;
301     
302     if (x < 4.) {
303         y = 31.754;
304     } else if (x < 6) {
305         Int_t j;
306         y = c[j = 4];
307         while (j > 0) y  = y * x + c[--j];
308     } else {
309         y =0.;
310     }
311     
312     return y;
313 }
314
315 Double_t AliGenMUONlib::YJpsiCDFscaled( Double_t *px, Double_t* dummy)
316 {
317     // J/Psi y 
318     return AliGenMUONlib::YJpsiPbPb(px, dummy);
319 }
320
321 Double_t AliGenMUONlib::YJpsiCDFscaledPP( Double_t *px, Double_t* dummy)
322 {
323     // J/Psi y 
324     return AliGenMUONlib::YJpsiPP(px, dummy);
325 }
326
327 Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t */*dummy*/)
328 {
329
330 //
331 // J/Psi y
332 //
333 //
334 // R. Vogt 2002
335 // pp 14  TeV
336 // MRST HO
337 // mc = 1.4 GeV, pt-kick 1 GeV
338 //
339
340     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
341     Double_t x = TMath::Abs(px[0]);
342     Double_t y;
343     
344     if (x < 2.5) {
345         y = 96.455 - 0.8483 * x * x;
346     } else if (x < 7.9) {
347         Int_t j;
348         y = c[j = 4];
349         while (j > 0) y  = y * x + c[--j];
350     } else {
351         y =0.;
352     }
353     
354     return y;
355 }
356
357 Double_t AliGenMUONlib::YJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
358 {
359
360 //
361 // J/Psi from B->J/Psi X
362 //
363 //
364     
365
366     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
367     
368     Double_t x = TMath::Abs(px[0]);
369     Double_t y;
370     
371     if (x > 6.) {
372         y = 0.;
373     } else {
374         Int_t j;
375         y = c[j = 6];
376         while (j > 0) y  = y * x + c[--j];
377     } 
378     
379     return y;
380 }
381
382
383
384 //                 particle composition
385 //
386 Int_t AliGenMUONlib::IpJpsi(TRandom *)
387 {
388 // J/Psi composition
389     return 443;
390 }
391 Int_t AliGenMUONlib::IpPsiP(TRandom *)
392 {
393 // Psi prime composition
394     return 100443;
395 }
396 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
397 {
398 // J/Psi composition
399   Int_t ip;
400   Float_t r = gRandom->Rndm();
401   if (r < 0.98) {
402     ip = 443;
403   } else {
404     ip = 100443;
405   }
406   return ip;
407 }
408
409
410
411 //                      Upsilon
412 //
413 //
414 //                  pt-distribution
415 //____________________________________________________________
416 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t */*dummy*/ )
417 {
418 // Upsilon pT
419   const Double_t kpt0 = 5.3;
420   const Double_t kxn  = 2.5;
421   Double_t x=*px;
422   //
423   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
424   return x/TMath::Power(pass1,kxn);
425 }
426
427 Double_t AliGenMUONlib::PtUpsilonCDFscaled( Double_t *px, Double_t */*dummy*/ )
428 {
429 // Upsilon pT
430   const Double_t kpt0 = 7.753;
431   const Double_t kxn  = 3.042;
432   Double_t x=*px;
433   //
434   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
435   return x/TMath::Power(pass1,kxn);
436 }
437
438 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( Double_t *px, Double_t */*dummy*/ )
439 {
440 // Upsilon pT
441 //
442 // pp 14 TeV
443 //
444 // scaled from CDF data at 2 TeV
445
446   const Double_t kpt0 = 8.610;
447   const Double_t kxn  = 3.051;
448   Double_t x=*px;
449   //
450   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
451   return x/TMath::Power(pass1,kxn);
452 }
453
454 Double_t AliGenMUONlib::PtUpsilonFlat( Double_t */*px*/, Double_t */*dummy*/ )
455 {
456   return 1.;
457 }
458
459 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
460 {
461
462 //
463 // Upsilon pT
464 //
465 //
466 // R. Vogt 2002
467 // PbPb 5.5 TeV
468 // MRST HO
469 // mc = 1.4 GeV, pt-kick 1 GeV
470 //
471     Float_t x = px[0];
472     Double_t c[8] = {
473         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,       
474         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
475     };
476     Double_t y;
477     if (x < 10.) {
478         Int_t j;
479         y = c[j = 7];
480         while (j > 0) y  = y * x +c[--j];
481         y = x * TMath::Exp(y);
482     } else {
483         y = 0.;
484     }
485     return y;
486 }
487
488 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t */*dummy*/)
489 {
490
491 //
492 // Upsilon pT
493 //
494 //
495 // R. Vogt 2002
496 // pp 14 TeV
497 // MRST HO
498 // mc = 1.4 GeV, pt-kick 1 GeV
499 //
500     Float_t x = px[0];
501     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,   
502                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
503     
504     Double_t y;
505     if (x < 10.) {
506         Int_t j;
507         y = c[j = 7];
508         while (j > 0) y  = y * x +c[--j];
509         y = x * TMath::Exp(y);
510     } else {
511         y = 0.;
512     }
513     return y;
514 }
515
516 //
517 //                    y-distribution
518 //
519 //____________________________________________________________
520 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t */*dummy*/)
521 {
522 // Upsilon y
523   const Double_t ky0 = 3.;
524   const Double_t kb=1.;
525   Double_t yu;
526   Double_t y=TMath::Abs(*py);
527   //
528   if (y < ky0)
529     yu=kb;
530   else
531     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
532   return yu;
533 }
534
535
536 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
537 {
538
539 //
540 // Upsilon y
541 //
542 //
543 // R. Vogt 2002
544 // PbPb 5.5 TeV
545 // MRST HO
546 // mc = 1.4 GeV, pt-kick 1 GeV
547 //
548
549     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
550                      -2.99753e-09, 1.28895e-05};
551         
552     Double_t x = px[0];
553     if (TMath::Abs(x) > 5.55) return 0.;
554     Int_t j;
555     Double_t y = c[j = 6];
556     while (j > 0) y  = y * x +c[--j];
557     return y;
558 }
559
560 Double_t AliGenMUONlib::YUpsilonCDFscaled( Double_t *px, Double_t *dummy)
561 {
562     // Upsilon y
563     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
564     
565 }
566
567 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( Double_t *px, Double_t *dummy)
568 {
569     // Upsilon y
570     return AliGenMUONlib::YUpsilonPP(px, dummy);
571     
572 }
573
574 Double_t AliGenMUONlib::YUpsilonFlat( Double_t */*px*/, Double_t */*dummy*/)
575 {
576     // Upsilon y
577     return 1.;
578     
579 }
580
581 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t */*dummy*/)
582 {
583
584 //
585 // Upsilon y
586 //
587 //
588 // R. Vogt 2002
589 // p p  14. TeV
590 // MRST HO
591 // mc = 1.4 GeV, pt-kick 1 GeV
592 //
593     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04, 
594                      -6.20539e-10, 1.29943e-05};
595                 
596     Double_t x = px[0];
597     if (TMath::Abs(x) > 6.2) return 0.;
598     Int_t j;
599     Double_t y = c[j = 6];
600     while (j > 0) y  = y * x +c[--j];
601     return y;
602 }
603
604 //                 particle composition
605 //
606 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
607 {
608 // y composition
609     return 553;
610 }
611 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
612 {
613 // y composition
614     return 100553;
615 }
616 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
617 {
618 // y composition
619     return 200553;
620 }
621 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
622 {
623 // y composition
624   Int_t ip;
625   Float_t r = gRandom->Rndm();
626   
627   if (r < 0.712) {
628     ip = 553;
629   } else if (r < 0.896) {
630     ip = 100553;
631   } else {
632     ip = 200553;
633   }
634   return ip;
635 }
636
637
638 //
639 //                        Phi
640 //
641 //
642 //    pt-distribution (by scaling of pion distribution)
643 //____________________________________________________________
644 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t */*dummy*/)
645 {
646 // Phi pT
647   return PtScal(*px,7);
648 }
649 //    y-distribution
650 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t */*dummy*/)
651 {
652 // Phi y
653     Double_t *dum=0;
654     return YJpsi(px,dum);
655 }
656 //                 particle composition
657 //
658 Int_t AliGenMUONlib::IpPhi(TRandom *)
659 {
660 // Phi composition
661     return 333;
662 }
663
664 //
665 //                        omega
666 //
667 //
668 //    pt-distribution (by scaling of pion distribution)
669 //____________________________________________________________
670 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t */*dummy*/)
671 {
672 // Omega pT
673   return PtScal(*px,5);
674 }
675 //    y-distribution
676 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t */*dummy*/)
677 {
678 // Omega y
679     Double_t *dum=0;
680     return YJpsi(px,dum);
681 }
682 //                 particle composition
683 //
684 Int_t AliGenMUONlib::IpOmega(TRandom *)
685 {
686 // Omega composition
687     return 223;
688 }
689
690
691 //
692 //                        Eta
693 //
694 //
695 //    pt-distribution (by scaling of pion distribution)
696 //____________________________________________________________
697 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t */*dummy*/)
698 {
699 // Eta pT
700   return PtScal(*px,3);
701 }
702 //    y-distribution
703 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t */*dummy*/)
704 {
705 // Eta y
706     Double_t *dum=0;
707     return YJpsi(px,dum);
708 }
709 //                 particle composition
710 //
711 Int_t AliGenMUONlib::IpEta(TRandom *)
712 {
713 // Eta composition
714     return 221;
715 }
716
717 //
718 //                        Charm
719 //
720 //
721 //                    pt-distribution
722 //____________________________________________________________
723 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t */*dummy*/)
724 {
725 // Charm pT
726   const Double_t kpt0 = 2.25;
727   const Double_t kxn  = 3.17;
728
729   Double_t x=*px;
730   //
731   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
732   return x/TMath::Power(pass1,kxn);
733 }
734
735 Double_t AliGenMUONlib::PtCharmCentral( Double_t *px, Double_t */*dummy*/)
736 {
737 // Charm pT
738   const Double_t kpt0 = 2.12;
739   const Double_t kxn  = 2.78;
740
741   Double_t x=*px;
742   //
743   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
744   return x/TMath::Power(pass1,kxn);
745 }
746 //                  y-distribution
747 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t */*dummy*/)
748 {
749 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
750 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
751 // shadowing + kt broadening 
752
753     Double_t x=px[0];
754     Double_t c[2]={-2.42985e-03,-2.31001e-04};
755     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
756     Double_t ycharm;
757     
758     if (TMath::Abs(x)>8) {
759       ycharm=0.;
760     }
761     else {
762       ycharm=TMath::Power(y,3);
763     }
764     
765     return ycharm;
766 }
767
768
769 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
770 {  
771 // Charm composition
772     Float_t random;
773     Int_t ip;
774 //    411,421,431,4122
775     random = ran->Rndm();
776 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
777 //  >>>>> cf. tab 4 p 11
778   
779     if (random < 0.30) {                       
780         ip=421;
781     } else if (random < 0.60) {
782         ip=-421;
783     } else if (random < 0.70) {
784         ip=411;
785     } else if (random < 0.80) {
786         ip=-411;
787     } else if (random < 0.86) {
788         ip=431;
789     } else if (random < 0.92) {
790         ip=-431;        
791     } else if (random < 0.96) {
792         ip=4122;
793     } else {
794         ip=-4122;
795     }
796     
797     return ip;
798 }
799
800 //
801 //                        Beauty
802 //
803 //
804 //                    pt-distribution
805 //____________________________________________________________
806 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t */*dummy*/)
807 {
808 // Beauty pT
809   const Double_t kpt0 = 6.53;
810   const Double_t kxn  = 3.59;
811   Double_t x=*px;
812   //
813   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
814   return x/TMath::Power(pass1,kxn);
815 }
816
817 Double_t AliGenMUONlib::PtBeautyCentral( Double_t *px, Double_t */*dummy*/)
818 {
819 // Beauty pT
820   const Double_t kpt0 = 6.14;
821   const Double_t kxn  = 2.93;
822   Double_t x=*px;
823   //
824   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
825   return x/TMath::Power(pass1,kxn);
826 }
827 //                     y-distribution
828 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t */*dummy*/)
829 {
830 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
831 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
832 // shadowing + kt broadening 
833
834     Double_t x=px[0];
835     Double_t c[2]={-1.27590e-02,-2.42731e-04};
836     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
837     Double_t ybeauty;
838     
839     if (TMath::Abs(x)>6) {
840       ybeauty=0.;
841     }
842     else {
843       ybeauty=TMath::Power(y,3);
844     }
845     
846     return ybeauty;
847 }
848
849
850 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
851 {  
852 // Beauty Composition
853     Float_t random;
854     Int_t ip;
855     random = ran->Rndm(); 
856     
857 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
858 //  >>>>> cf. tab 4 p 11
859     
860  if (random < 0.20) {                       
861         ip=511;
862     } else if (random < 0.40) {
863         ip=-511;
864     } else if (random < 0.605) {
865         ip=521;
866     } else if (random < 0.81) {
867         ip=-521;
868     } else if (random < 0.87) {
869         ip=531;
870     } else if (random < 0.93) {
871         ip=-531;        
872     } else if (random < 0.965) {
873         ip=5122;
874     } else {
875         ip=-5122;
876     }
877     
878  return ip;
879 }
880
881
882 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
883 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
884 {
885 // Return pointer to pT parameterisation
886     TString sname = TString(tname);
887     GenFunc func;
888     switch (param) 
889     {
890     case kPhi:
891         func=PtPhi;
892         break;
893     case kOmega:
894         func=PtOmega;
895         break;
896     case kEta:
897         func=PtEta;
898         break;
899     case kJpsiFamily:
900     case kPsiP:
901     case kJpsi:
902         if (sname == "Vogt" || sname == "Vogt PbPb") {
903             func=PtJpsiPbPb;
904         } else if (sname == "Vogt pp") {
905             func=PtJpsiPP;
906         } else if (sname == "CDF scaled") {
907             func=PtJpsiCDFscaled;
908         } else if (sname == "CDF pp") {
909             func=PtJpsiCDFscaledPP;
910         } else if (sname == "Flat") {
911             func=PtJpsiFlat;
912         } else {
913             func=PtJpsi;
914         }
915         break;
916     case kJpsiFromB:
917         func = PtJpsiBPbPb;
918         break;
919     case kUpsilonFamily:
920     case kUpsilonP:
921     case kUpsilonPP:
922     case kUpsilon:
923         if (sname == "Vogt" || sname == "Vogt PbPb") {
924             func=PtUpsilonPbPb;
925         } else if (sname == "Vogt pp") {
926             func=PtUpsilonPP;
927         } else if (sname == "CDF scaled") {
928             func=PtUpsilonCDFscaled;
929         } else if (sname == "CDF pp") {
930             func=PtUpsilonCDFscaledPP;
931         } else if (sname == "Flat") {
932             func=PtUpsilonFlat;
933         } else {
934             func=PtUpsilon;
935         }
936         break;  
937     case kCharm:
938         if (sname == "central") {
939             func=PtCharmCentral;
940         } else {
941             func=PtCharm;
942         }
943         break;
944     case kBeauty:
945         if (sname == "central") {
946             func=PtBeautyCentral;
947         } else {
948             func=PtBeauty;
949         }
950         break;
951     case kPion:
952         func=PtPion;
953         break;
954     case kKaon:
955         func=PtKaon;
956         break;
957     default:
958         func=0;
959         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
960     }
961     return func;
962 }
963
964 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
965 {
966   //    
967   // Return pointer to y- parameterisation
968   //
969     TString sname = TString(tname);
970     GenFunc func;
971     switch (param) 
972     {
973     case kPhi:
974         func=YPhi;
975         break;
976     case kEta:
977         func=YEta;
978         break;
979     case kOmega:
980         func=YOmega;
981         break;
982     case kJpsiFamily:
983     case kPsiP:
984     case kJpsi:
985         if (sname == "Vogt" || sname == "Vogt PbPb") {
986             func=YJpsiPbPb;
987         } else if (sname == "Vogt pp"){
988             func=YJpsiPP;
989         } else if (sname == "CDF scaled") {
990             func=YJpsiCDFscaled;
991         } else if (sname == "CDF pp") {
992             func=YJpsiCDFscaledPP;
993         } else if (sname == "Flat") {
994             func=YJpsiFlat;
995         } else {
996             func=YJpsi;
997         }
998         break;
999     case kJpsiFromB:
1000         func = YJpsiBPbPb;
1001         break;
1002     case kUpsilonFamily:
1003     case kUpsilonP:
1004     case kUpsilonPP:
1005     case kUpsilon:
1006         if (sname == "Vogt" || sname == "Vogt PbPb") {
1007             func=YUpsilonPbPb;
1008         } else if (sname == "Vogt pp") {
1009             func = YUpsilonPP;
1010         } else if (sname == "CDF scaled") {
1011             func=YUpsilonCDFscaled;
1012         } else if (sname == "CDF pp") {
1013             func=YUpsilonCDFscaledPP;
1014         } else if (sname == "Flat") {
1015             func=YUpsilonFlat;
1016         } else {
1017             func=YUpsilon;
1018         }
1019         break;
1020     case kCharm:
1021         func=YCharm;
1022         break;
1023     case kBeauty:
1024         func=YBeauty;
1025         break;
1026     case kPion:
1027         func=YPion;
1028         break;
1029     case kKaon:
1030         func=YKaon;
1031         break;
1032     default:
1033         func=0;
1034         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
1035     }
1036     return func;
1037 }
1038 typedef Int_t (*GenFuncIp) (TRandom *);
1039 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
1040 {
1041 // Return pointer to particle type parameterisation
1042     GenFuncIp func;
1043     switch (param) 
1044     {
1045     case kPhi:
1046         func=IpPhi;
1047         break;
1048     case kEta:
1049         func=IpEta;
1050         break;
1051     case kOmega:
1052         func=IpOmega;
1053         break;
1054     case kJpsiFamily:
1055         func=IpJpsiFamily;
1056         break;
1057     case kPsiP:
1058         func=IpPsiP;
1059         break;
1060     case kJpsi:
1061     case kJpsiFromB:
1062         func=IpJpsi;
1063         break;
1064     case kUpsilon:
1065         func=IpUpsilon;
1066         break;
1067     case kUpsilonFamily:
1068       func=IpUpsilonFamily;
1069       break;
1070     case kUpsilonP:
1071         func=IpUpsilonP;
1072         break;
1073     case kUpsilonPP:
1074         func=IpUpsilonPP;
1075         break;
1076     case kCharm:
1077         func=IpCharm;
1078         break;
1079     case kBeauty:
1080         func=IpBeauty;
1081         break;
1082     case kPion:
1083         func=IpPion;
1084         break;
1085     case kKaon:
1086         func=IpKaon;
1087         break;
1088     default:
1089         func=0;
1090         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
1091     }
1092     return func;
1093 }
1094
1095
1096
1097 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0, 
1098                                    Float_t dx,
1099                                    Int_t n, Int_t no)
1100 {
1101 //
1102 // Neville's alorithm for interpolation
1103 //
1104 // x:  x-value
1105 // y:  Input array
1106 // x0: minimum x 
1107 // dx: step size
1108 //  n: number of data points
1109 // no: order of polynom 
1110 //
1111     Float_t*  c = new Float_t[n];
1112     Float_t*  d = new Float_t[n];
1113     Int_t m, i;
1114     for (i = 0; i < n; i++) {
1115         c[i] = y[i];
1116         d[i] = y[i];
1117     }
1118     
1119     Int_t   ns  = int((x - x0)/dx);
1120     
1121     Float_t y1  = y[ns];
1122     ns--;    
1123     for (m = 0; m < no; m++) {
1124         for (i = 0; i < n-m; i++) {     
1125             Float_t ho = x0 + Float_t(i) * dx - x;
1126             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
1127             Float_t w  = c[i+1] - d[i];
1128             Float_t den = ho-hp;
1129             den = w/den;
1130             d[i] = hp * den;
1131             c[i] = ho * den;
1132         }
1133         Float_t dy;
1134         
1135         if (2*ns < (n-m-1)) {
1136             dy  = c[ns+1];
1137         } else {
1138             dy  = d[ns--];
1139         }
1140         y1 += dy;}
1141     delete[] c;
1142     delete[] d;
1143
1144     return y1;
1145 }
1146
1147