Updated version.
[u/mrichter/AliRoot.git] / EVGEN / AliGenGSIlib.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 /////////////////////////////////////////////////////////////////////////////
19 //                                                                         //
20 // Implementation of AliGenlib to collect parametrisations used for        //
21 // GSI simulations.                                                        //
22 // It is an extension of AliMUONLib providing in addition the option       //
23 // for different parametrisations of pt, y and ip for every particle type  //
24 //                                                                         //
25 // Responsible: Andres.Sandoval@cern.ch                                    //
26 //                                                                         //
27 /////////////////////////////////////////////////////////////////////////////
28
29 #include "TMath.h"
30 #include "TRandom.h"
31 #include "TString.h"
32 #include "AliGenGSIlib.h"
33
34
35 ClassImp(AliGenGSIlib)
36
37 //==========================================================================
38 //
39 //              Definition of Particle Distributions
40 //                    
41 //==========================================================================
42 //
43 //                         Upsilon 
44 //
45 //--------------------------------------------------------------------------
46 //
47 //                  upsilon particle composition
48 //
49 //--------------------------------------------------------------------------
50 Int_t AliGenGSIlib::IpUpsilon(TRandom *)
51 {
52 // Return upsilon pdg code
53
54   return 553;     
55
56 }
57 Double_t AliGenGSIlib::PtUpsilonFlat( const Double_t *px, const Double_t */*dummy*/ )
58 {
59 //--------------------------------------------------------------------------
60 //
61 //                   upsilon pt-distribution FLAT
62 //
63 //____________________________________________________________--------------
64     
65   const Double_t kptmin = 0.0;
66   const Double_t kptmax  = 15.0;
67   Double_t x=*px;
68   Double_t weight = 0.;
69
70   if ((x > kptmin) &&  (x < kptmax)) weight = 1.;
71
72   return weight;
73    
74 }
75 Double_t AliGenGSIlib::YUpsilonFlat(const Double_t *py, const Double_t */*dummy*/)
76 {
77 //--------------------------------------------------------------------------
78 //
79 //                    upsilon y-distribution FLAT
80 //
81 //--------------------------------------------------------------------------
82
83   const Double_t ky0 = 1.5;
84   const Double_t kb=1.;
85   Double_t yu;
86   Double_t y=TMath::Abs(*py);
87
88   if (y < ky0)
89     yu=kb;
90   else
91     yu = 0.;
92
93   return yu;
94
95 }
96 Double_t AliGenGSIlib::PtUpsilonRitman( const Double_t *px, const Double_t */*dummy*/ )
97 {
98 //--------------------------------------------------------------------------
99 //
100 //                 upsilon pt-distribution RITMAN
101 //
102 //--------------------------------------------------------------------------
103
104   const Double_t kpt0 = 4.7;
105   const Double_t kxn  = 3.5;
106   Double_t x=*px;
107
108   Double_t pass1 = 1.+((x*x)/(kpt0*kpt0));
109
110   return x/TMath::Power(pass1,kxn);
111    
112 }
113 Double_t AliGenGSIlib::YUpsilonRitman(const Double_t *py, const Double_t */*dummy*/)
114 {
115 //--------------------------------------------------------------------------
116 //
117 //                  upsilon y-distribution RITMAN
118 //
119 //--------------------------------------------------------------------------
120
121   const Double_t ky0 = 3.;
122   const Double_t kb=1.;
123   Double_t yu;
124   Double_t y=TMath::Abs(*py);
125
126   if (y < ky0)
127     yu=kb;
128   else
129     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
130
131   return yu;
132    
133 }
134 Double_t AliGenGSIlib::PtUpsilonKarel( const Double_t */*px*/, const Double_t */*dummy*/ )
135 {
136 //--------------------------------------------------------------------------
137 //
138 //                 upsilon pt-distribution kAREL
139 //
140 //--------------------------------------------------------------------------
141 // to implement
142
143   return 0.1;   
144
145 }
146 Double_t AliGenGSIlib::YUpsilonKarel(const Double_t */*py*/, const Double_t */*dummy*/)
147 {
148 //--------------------------------------------------------------------------
149 //
150 //                  upsilon y-distribution KAREL
151 //
152 //--------------------------------------------------------------------------
153   
154   //to implement
155
156   return 0.2;  
157
158 }
159 Double_t AliGenGSIlib::PtUpsilonMUON( const Double_t *px, const Double_t */*dummy*/ )
160 {
161 //--------------------------------------------------------------------------
162 //
163 //                 upsilon pt-distribution MUONlib
164 //
165 //--------------------------------------------------------------------------
166
167   const Double_t kpt0 = 5.3;
168   const Double_t kxn  = 2.5;
169   Double_t x=*px;
170
171   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
172
173   return x/TMath::Power(pass1,kxn);
174
175 }
176 Double_t AliGenGSIlib::YUpsilonMUON(const Double_t *py, const Double_t */*dummy*/)
177 {
178 //--------------------------------------------------------------------------
179 //
180 //                   upsilon y-distribution MUONlib
181 //
182 //--------------------------------------------------------------------------
183
184   const Double_t ky0 = 3.;
185   const Double_t kb=1.;
186   Double_t yu;
187   Double_t y=TMath::Abs(*py);
188
189   if (y < ky0)
190     yu=kb;
191   else
192     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
193
194   return yu;
195
196 }
197 //--------------------------------------------------------------------------
198 //
199 //                             J/Psi
200 //
201 Int_t AliGenGSIlib::IpJpsi(TRandom *)
202 {
203 //--------------------------------------------------------------------------
204 //
205 //                    J/Psi particle composition
206 //
207 //--------------------------------------------------------------------------
208
209   return 443;     
210
211 }
212 Double_t AliGenGSIlib::PtJpsiFlat( const Double_t *px, const Double_t */*dummy*/ )
213 {
214 //--------------------------------------------------------------------------
215 //
216 //                   J/Psi pt-distribution FLAT
217 //
218 //--------------------------------------------------------------------------
219
220   const Double_t kptmin = 0.0;
221   const Double_t kptmax  = 15.0;
222   Double_t x=*px;
223   Double_t weight = 0.;
224
225   if ((x > kptmin) && (x < kptmax)) weight = 1.;
226
227   return weight;
228    
229 }
230 Double_t AliGenGSIlib::YJpsiFlat(const Double_t *py, const Double_t */*dummy*/)
231 {
232 //--------------------------------------------------------------------------
233 //
234 //                    J/Psi y-distribution FLAT
235 //
236 //--------------------------------------------------------------------------
237
238   const Double_t ky0 = 1.5;
239   const Double_t kb=1.;
240   Double_t yu;
241   Double_t y=TMath::Abs(*py);
242
243   if (y < ky0)
244     yu=kb;
245   else
246     yu = 0.;
247
248   return yu;
249
250 }
251 Double_t AliGenGSIlib::PtJpsiMUON( const Double_t *px, const Double_t */*dummy*/ )
252 {
253 //--------------------------------------------------------------------------
254 //
255 //                   J/Psi pt-distribution MUONlib
256 //
257 //--------------------------------------------------------------------------
258
259   const Double_t kpt0 = 4.;
260   const Double_t kxn  = 3.6;
261   Double_t x=*px;
262
263   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
264   return x/TMath::Power(pass1,kxn);
265    
266 }
267 Double_t AliGenGSIlib::PtJpsiRitman( const Double_t *px, const Double_t */*dummy*/ )
268 {
269 //--------------------------------------------------------------------------
270 //
271 //                   J/Psi pt-distribution Ritman
272 //
273 //--------------------------------------------------------------------------
274
275   const Double_t kpt0 = 2.3;
276   const Double_t kxn  = 3.5;
277   Double_t x=*px;
278
279   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
280
281   return x/TMath::Power(pass1,kxn);
282    
283 }
284 Double_t AliGenGSIlib::YJpsiMUON(const Double_t *py, const Double_t */*dummy*/)
285 {
286 //--------------------------------------------------------------------------
287 //
288 //                    J/Psi y-distribution MUONlib
289 //
290 //--------------------------------------------------------------------------
291
292   const Double_t ky0 = 4.;
293   const Double_t kb=1.;
294   Double_t yj;
295   Double_t y=TMath::Abs(*py);
296
297   if (y < ky0)
298     yj=kb;
299   else
300     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
301   return yj;
302
303 }
304 //--------------------------------------------------------------------------
305 //
306 //                  J/Psi pt-distribution by Sergei
307 //
308 //--------------------------------------------------------------------------
309 //Double_t AliGenGSIlib::PtJpsi( Double_t *px, Double_t */*dummy*/ )
310 //{
311
312 //  return x = gRandom->Rndm()*10.;
313
314 //}
315 //--------------------------------------------------------------------------
316 //
317 //                  J/Psi y-distribution by Sergei
318 //
319 //--------------------------------------------------------------------------
320 /*Double_t AliGenGSIlib::YJpsi(Double_t *py, Double_t *dummy)
321 {
322
323   const Double_t ky0 = 4.;
324   const Double_t kb=1.;
325   Double_t yj;
326   Double_t y=TMath::Abs(*py);
327   //
328   if (y < ky0)
329     yj=kb;
330   else
331     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
332   return yj;
333
334 }
335 */
336 //--------------------------------------------------------------------------
337 //
338 //                        Charm
339 //
340 //--------------------------------------------------------------------------
341 Int_t AliGenGSIlib::IpCharm(TRandom *ran)
342 {
343 //
344 //                    charm particle composition
345 //
346 //--------------------------------------------------------------------------
347   
348     Float_t random;
349     Int_t ip;
350     // 411,421,431,4122
351     random = ran->Rndm();
352     if (random < 0.5) {
353         ip=411;
354     } else if (random < 0.75) {
355         ip=421;
356     } else if (random < 0.90) {
357         ip=431;
358     } else {
359         ip=4122;
360     }
361     if (ran->Rndm() < 0.5) {ip=-ip;}
362     
363     return ip;
364
365 }
366 Double_t AliGenGSIlib::PtCharmFlat( const Double_t *px, const Double_t */*dummy*/)
367 {
368 //--------------------------------------------------------------------------
369 //
370 //                    charm pt-distribution, FLAT
371 //
372 //--------------------------------------------------------------------------
373
374   Double_t x=*px;
375
376   if (x>10.)  x = 0.;
377   else x=1.;
378   return x ;
379
380 }
381 Double_t AliGenGSIlib::PtCharmGSI( const Double_t *px, const Double_t */*dummy*/)
382 {
383 //--------------------------------------------------------------------------
384 //
385 //                    charm pt-distribution, from Dariuzs Miskowiec
386 //
387 //--------------------------------------------------------------------------
388
389   //Taken from PYTHIA with MRS D-' (3031 from PDFLIB), K=3.0
390   const Double_t kp1 = 1.3;
391   const Double_t kp2  = 0.39;
392   const Double_t kp3  = 0.018;
393   const Double_t kp4  = 0.91;
394   Double_t x=*px;
395
396   Double_t pass1 =TMath::Exp(-x/kp2)  ;
397   Double_t pass2 =TMath::Exp(-x/kp4)  ;
398   return TMath::Power(x,kp1) * (pass1 + kp3 * pass2);
399
400 }
401 Double_t AliGenGSIlib::PtCharmMUON( const Double_t *px, const Double_t */*dummy*/)
402 {
403 //--------------------------------------------------------------------------
404 //
405 //                    charm pt-distribution, from MUONlib
406 //
407 //--------------------------------------------------------------------------
408
409   const Double_t kpt0 = 4.08;
410   const Double_t kxn  = 9.40;
411   Double_t x=*px;
412
413   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
414
415   return x/TMath::Power(pass1,kxn);
416
417 }
418 Double_t AliGenGSIlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
419 {
420 //--------------------------------------------------------------------------
421 //
422 //                    charm y-distribution
423 //
424 //--------------------------------------------------------------------------
425
426     Double_t *dum=0;
427
428     return YJpsiMUON(px,dum);
429
430 }
431 //--------------------------------------------------------------------------
432 //
433 //                        Beauty
434 //
435 //--------------------------------------------------------------------------
436 Int_t AliGenGSIlib::IpBeauty(TRandom *ran)
437 {  
438 //
439 //                    beauty particle composition
440 //
441 //--------------------------------------------------------------------------
442
443     Float_t random;
444     Int_t ip;
445     random = ran->Rndm();
446     if (random < 0.5) {
447         ip=511;
448     } else if (random < 0.75) {
449         ip=521;
450     } else if (random < 0.90) {
451         ip=531;
452     } else {
453         ip=5122;
454     }
455     if (ran->Rndm() < 0.5) {ip=-ip;}
456     
457     return ip;
458 }
459 Double_t AliGenGSIlib::PtBeautyFlat( const Double_t *px, const Double_t */*dummy*/)
460 {
461 //--------------------------------------------------------------------------
462 //
463 //                    beauty pt-distribution, FLAT
464 //
465 //--------------------------------------------------------------------------
466
467   Double_t x=*px;
468
469   if (x>10.) x=0.;
470   else x = 1.;
471   return x ;
472
473 }
474 Double_t AliGenGSIlib::PtBeautyGSI( const Double_t *px, const Double_t */*dummy*/)
475 {
476 //--------------------------------------------------------------------------
477 //
478 //
479 //                    beauty pt-distribution, from D. Miskowiec
480 //
481 //--------------------------------------------------------------------------
482
483   //Taken from PYTHIA with MRS D-' (3031 from PDFLIB), K=3.0
484   const Double_t kp1 = 1.3;
485   const Double_t kp2  = 1.78;
486   const Double_t kp3  = 0.0096;
487   const Double_t kp4  = 4.16;
488   Double_t x=*px;
489
490   Double_t pass1 =TMath::Exp(-x/kp2)  ;
491   Double_t pass2 =TMath::Exp(-x/kp4)  ;
492
493   return TMath::Power(x,kp1) * (pass1 + kp3 * pass2);
494
495 }
496 Double_t AliGenGSIlib::PtBeautyMUON( const Double_t *px, const Double_t */*dummy*/)
497 {
498 //--------------------------------------------------------------------------
499 //
500 //                    beauty pt-distribution, from MUONlib
501 //
502 //--------------------------------------------------------------------------
503
504   const Double_t kpt0 = 4.;
505   const Double_t kxn  = 3.6;
506   Double_t x=*px;
507
508   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
509
510   return x/TMath::Power(pass1,kxn);
511
512 }
513 Double_t AliGenGSIlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
514 {
515 //--------------------------------------------------------------------------
516 //
517 //                    beauty y-distribution
518 //
519 //--------------------------------------------------------------------------
520
521     Double_t *dum=0;
522
523     return YJpsiMUON(px,dum);
524
525 }
526 //--------------------------------------------------------------------------
527 //
528 //                          Eta
529 //
530 //--------------------------------------------------------------------------
531 Int_t AliGenGSIlib::IpEta(TRandom *)
532 {
533 //
534 //                 eta particle composition 
535 //
536 //--------------------------------------------------------------------------
537
538   return 221;     
539
540 }
541 Double_t AliGenGSIlib::PtEtaPHOS( const Double_t *px, const Double_t */*dummy*/ )
542 {
543 //--------------------------------------------------------------------------
544 //
545 //                  eta pt-distribution
546 //
547 //____________________________________________________________--------------
548
549   return PtScal(*px,3);  //  3==> Eta in the PtScal function
550    
551 }
552 Double_t AliGenGSIlib::YEtaPHOS(const Double_t *py, const Double_t */*dummy*/)
553 {
554 //--------------------------------------------------------------------------
555 //
556 //                   eta y-distribution 
557 //
558 //--------------------------------------------------------------------------
559
560   const Double_t ka    = 1000.;
561   const Double_t kdy   = 4.;
562
563   Double_t y=TMath::Abs(*py);
564
565   Double_t ex = y*y/(2*kdy*kdy);
566
567   return ka*TMath::Exp(-ex);
568
569 }
570 //--------------------------------------------------------------------------
571 //
572 //                       Etaprime
573 //
574 //--------------------------------------------------------------------------
575 Int_t AliGenGSIlib::IpEtaprime(TRandom *)
576 {
577 //
578 //                 etaprime particle composition 
579 //
580 //--------------------------------------------------------------------------
581
582   return 331;     
583
584 }
585 Double_t AliGenGSIlib::PtEtaprimePHOS( const Double_t *px, const Double_t */*dummy*/ )
586 {
587 //--------------------------------------------------------------------------
588 //
589 //                 etaprime pt-distribution
590 //
591 //____________________________________________________________--------------
592
593   return PtScal(*px,5);  //  5==> Etaprime in the PtScal function
594    
595 }
596 Double_t AliGenGSIlib::YEtaprimePHOS(const Double_t *py, const Double_t */*dummy*/)
597 {
598 //--------------------------------------------------------------------------
599 //
600 //                  etaprime y-distribution 
601 //
602 //--------------------------------------------------------------------------
603
604   const Double_t ka    = 1000.;
605   const Double_t kdy   = 4.;
606
607   Double_t y=TMath::Abs(*py);
608
609   Double_t ex = y*y/(2*kdy*kdy);
610
611   return ka*TMath::Exp(-ex);
612
613 }
614 //--------------------------------------------------------------------------
615 //
616 //                       omega 
617 //
618 //--------------------------------------------------------------------------
619 Int_t AliGenGSIlib::IpOmega(TRandom *)
620 {
621 //
622 //                 omega particle composition 
623 //
624 //--------------------------------------------------------------------------
625
626   return 223;     
627
628 }
629 Double_t AliGenGSIlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
630 {
631 //--------------------------------------------------------------------------
632 //
633 //                  omega pt-distribution
634 //
635 //____________________________________________________________--------------
636
637   return PtScal(*px,4);  //  4==> Omega in the PtScal function
638    
639 }
640 Double_t AliGenGSIlib::YOmega(const Double_t *py, const Double_t */*dummy*/)
641 {
642 //--------------------------------------------------------------------------
643 //
644 //                   omega y-distribution 
645 //
646 //--------------------------------------------------------------------------
647
648   const Double_t ka    = 1000.;
649   const Double_t kdy   = 4.;
650
651
652   Double_t y=TMath::Abs(*py);
653
654   Double_t ex = y*y/(2*kdy*kdy);
655
656   return ka*TMath::Exp(-ex);
657
658 }
659 //--------------------------------------------------------------------------
660 //
661 //                       Rho 
662 //
663 //--------------------------------------------------------------------------
664
665 Int_t AliGenGSIlib::IpRho(TRandom *)
666 {
667 //
668 //                 rho particle composition 
669 //
670 //--------------------------------------------------------------------------
671
672   return 113;     
673
674 }
675 Double_t AliGenGSIlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
676 {
677 //--------------------------------------------------------------------------
678 //
679 //                  rho pt-distribution
680 //
681 //____________________________________________________________--------------
682
683   return PtScal(*px,11);  //  11==> Rho in the PtScal function
684    
685 }
686 Double_t AliGenGSIlib::YRho(const Double_t *py, const Double_t */*dummy*/)
687 {
688 //--------------------------------------------------------------------------
689 //
690 //                   rho y-distribution 
691 //
692 //--------------------------------------------------------------------------
693
694   const Double_t ka    = 1000.;
695   const Double_t kdy   = 4.;
696
697
698   Double_t y=TMath::Abs(*py);
699
700   Double_t ex = y*y/(2*kdy*kdy);
701
702   return ka*TMath::Exp(-ex);
703
704 }
705 //--------------------------------------------------------------------------
706 //
707 //                              Pion
708 //
709 //--------------------------------------------------------------------------
710 Int_t AliGenGSIlib::IpPionPHOS(TRandom *ran)
711 {
712 //
713 //                 particle composition  pi+, pi0, pi-
714 //
715 //--------------------------------------------------------------------------
716
717     Float_t random = ran->Rndm();
718
719     if ( (3.*random)  < 1. ) 
720     {
721           return 211 ;
722     } 
723     else
724     {  
725       if ( (3.*random) >= 2.)
726       {
727          return -211 ;
728       }
729       else 
730       {
731         return 111  ;
732       }
733     }
734 }
735 Double_t AliGenGSIlib::PtPion( const Double_t *px, const Double_t */*dummy*/ )
736 {
737 //--------------------------------------------------------------------------
738 //
739 //                  pion pt-distribution
740 //
741 //       Pion transverse momentum distribtuion as in AliGenMUONlib class, 
742 //       version 3.01 of aliroot:
743 //       PT-PARAMETERIZATION CDF, PRL 61(88) 1819
744 //       POWER LAW FOR PT > 500 MEV
745 //       MT SCALING BELOW (T=160 MEV)
746 //
747 //____________________________________________________________--------------
748
749   const Double_t kp0 = 1.3;
750   const Double_t kxn = 8.28;
751   const Double_t kxlim=0.5;
752   const Double_t kt=0.160;
753   const Double_t kxmpi=0.139;
754   const Double_t kb=1.;
755   Double_t y, y1, kxmpi2, ynorm, a;
756   Double_t x=*px;
757
758   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
759   kxmpi2=kxmpi*kxmpi;
760   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+kxmpi2)/kt));
761   a=ynorm/y1;
762   if (x > kxlim)
763     y=a*TMath::Power(kp0/(kp0+x),kxn);
764   else
765     y=kb*TMath::Exp(-sqrt(x*x+kxmpi2)/kt);
766   return y*x;
767    
768 }
769 Double_t AliGenGSIlib::YPion(const Double_t *py, const Double_t */*dummy*/)
770 {
771 //--------------------------------------------------------------------------
772 //
773 //                    pion y-distribution 
774 //
775 //--------------------------------------------------------------------------
776
777   const Double_t ka    = 7000.;   
778   const Double_t kdy   = 4.;
779
780   Double_t y=TMath::Abs(*py);
781
782   Double_t ex = y*y/(2*kdy*kdy);
783
784   return ka*TMath::Exp(-ex);
785
786 }
787 Int_t AliGenGSIlib::IpKaonPHOS(TRandom *ran)
788 {
789 //--------------------------------------------------------------------------
790 //
791 //
792 //                        Kaon
793 //--------------------------------------------------------------------------
794 //
795 //                kaon particle composition K+, K-, Ko_short, Ko_long
796 //
797 //--------------------------------------------------------------------------
798
799     Float_t random = ran->Rndm();
800     Float_t random2 = ran->Rndm();
801     if (random2 < 0.5) 
802     {
803       if (random < 0.5) {       
804         return  321;   //   K+
805       } else {
806         return -321;   // K-
807       }
808     }
809     else
810     {  
811       if (random < 0.5) {       
812         return  130;   // K^0 short
813       } else {  
814         return  310;   // K^0 long
815       }
816     }
817 }
818 Double_t AliGenGSIlib::PtKaonPHOS( const Double_t *px, const Double_t */*dummy*/ )
819 {
820 //--------------------------------------------------------------------------
821 //
822 //                   kaon pt-distribution
823 //
824 //____________________________________________________________--------------
825
826   return PtScal(*px,2);  //  2==> Kaon in the PtScal function
827    
828 }
829 Double_t AliGenGSIlib::YKaonPHOS(const Double_t *py, const Double_t */*dummy*/)
830 {
831 //--------------------------------------------------------------------------
832 //
833 //                    kaon y-distribution 
834 //
835 //--------------------------------------------------------------------------
836
837   const Double_t ka    = 1000.;
838   const Double_t kdy   = 4.;
839
840   Double_t y=TMath::Abs(*py);
841
842   Double_t ex = y*y/(2*kdy*kdy);
843
844   return ka*TMath::Exp(-ex);
845
846 }
847 //--------------------------------------------------------------------------
848 //
849 //                        Phi  
850 //
851 Int_t AliGenGSIlib::IpPhi(TRandom *)
852 {
853 //--------------------------------------------------------------------------
854 //
855 //                 particle composition 
856 //
857 //--------------------------------------------------------------------------
858
859   return 333;     
860
861 }
862 Double_t AliGenGSIlib::PtPhiPHOS( const Double_t *px, const Double_t */*dummy*/ )
863 {
864 //--------------------------------------------------------------------------
865 //
866 //                   phi pt-distribution
867 //
868 //____________________________________________________________--------------
869
870   return PtScal(*px,6);  //  6==> Phi in the PtScal function
871    
872 }
873 Double_t AliGenGSIlib::YPhiPHOS(const Double_t *py, const Double_t */*dummy*/)
874 {
875 //--------------------------------------------------------------------------
876 //
877 //                    phi y-distribution 
878 //
879 //--------------------------------------------------------------------------
880
881   const Double_t ka    = 1000.;
882   const Double_t kdy   = 4.;
883
884
885   Double_t y=TMath::Abs(*py);
886  
887   Double_t ex = y*y/(2*kdy*kdy);
888
889   return ka*TMath::Exp(-ex);
890
891 }
892 Int_t AliGenGSIlib::IpBaryons(TRandom *ran)
893 {
894 //--------------------------------------------------------------------------
895 //
896 //                          Baryons
897 //
898 //--------------------------------------------------------------------------
899 //
900 //                 baryons particle composition p, pbar, n, nbar
901 //
902 //--------------------------------------------------------------------------
903
904     Float_t random = ran->Rndm();
905     Float_t random2 = ran->Rndm();
906     if (random2 < 0.5) 
907     {
908       if (random < 0.5) {       
909         return  2212;   //   p
910       } else {
911         return -2212;   // pbar
912       }
913     }
914     else
915     {  
916       if (random < 0.5) {       
917         return  2112;   // n
918       } else {  
919         return -2112;   // n bar
920       }
921     }
922 }
923 Double_t AliGenGSIlib::PtBaryons( const Double_t *px, const Double_t */*dummy*/ )
924 {
925 //--------------------------------------------------------------------------
926 //
927 //                  baryons pt-distribution
928 //
929 //____________________________________________________________--------------
930
931   return PtScal(*px,7);  //  7==> Baryon in the PtScal function
932    
933 }
934 Double_t AliGenGSIlib::YBaryons(const Double_t *py, const Double_t */*dummy*/)
935 {
936 //--------------------------------------------------------------------------
937 //
938 //                   baryons y-distribution 
939 //
940 //--------------------------------------------------------------------------
941
942   const Double_t ka    = 1000.;
943   const Double_t kdy   = 4.;
944
945   Double_t y=TMath::Abs(*py);
946
947   Double_t ex = y*y/(2*kdy*kdy);
948
949   return ka*TMath::Exp(-ex);
950
951 }
952 //============================================================= 
953 //
954 //                    Mt-scaling as in AliGenPHOSlib  
955 //
956 //============================================================= 
957 //
958  Double_t AliGenGSIlib::PtScal(Double_t pt, Int_t np)
959 {
960 // Function for the calculation of the Pt distribution for a 
961 // given particle np, from the pion Pt distribution using the 
962 // mt scaling. 
963
964 // It was taken from AliGenPHOSlib aliroot version 3.04, which 
965 // is an update of the one in AliGenMUONlib aliroot version 3.01
966 // with an extension for Baryons but supressing Rhos
967 // np = 1=>Pions 2=>Kaons 3=>Etas 4=>Omegas 5=>ETA' 6=>PHI
968 //      7=>BARYONS-BARYONBARS
969
970 // The present adds the Rhos
971
972 // MASS   1=>PI, 2=>K, 3=>ETA, 4=>OMEGA, 5=>ETA', 6=>PHI 
973 //        7=>BARYON-BARYONBAR, 11==>RHO
974
975   const Double_t khm[11] = {0.1396, 0.494,  0.547,    0.782,   0.957,   1.02, 
976                                          0.938, 0. , 0., 0., 0.769};
977
978   //     VALUE MESON/PI AT 5 GEV
979
980   const Double_t kfmax[11]={1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
981
982   np--;
983   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
984   Double_t kfmax2=f5/kfmax[np];
985   // PIONS
986   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
987   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
988                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
989   return fmtscal*ptpion;
990 }
991
992 //==========================================================================
993 //
994 //                     Set Getters 
995 //    
996 //==========================================================================
997
998 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
999
1000 typedef Int_t (*GenFuncIp) (TRandom *);
1001
1002 GenFunc AliGenGSIlib::GetPt(Int_t param, const char * tname) const
1003 {
1004 // Return pointer to pT parameterisation
1005    GenFunc func=0;
1006    TString sname(tname);
1007
1008    switch (param) 
1009     {
1010     case kUpsilon:
1011       if (sname=="FLAT"){
1012         func= PtUpsilonFlat;
1013         break;
1014       }
1015       if (sname=="MUON"){
1016         func= PtUpsilonMUON;
1017         break;
1018       }
1019       if (sname=="RITMAN"){
1020         func=PtUpsilonRitman;
1021         break;
1022       }
1023       if (sname=="KAREL"){
1024         func=PtUpsilonKarel;
1025         break;
1026       }
1027       func=0;
1028       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1029       break;
1030
1031     case kJPsi:
1032       if (sname=="FLAT"){
1033         func= PtJpsiFlat;
1034         break;
1035       }
1036       if (sname=="MUON"){
1037         func= PtJpsiMUON;
1038         break;
1039       }
1040       //      if (sname=="SERGEI"){
1041       //        func= PtJpsi;
1042       //        break;
1043       //     }
1044       func=0;
1045       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1046       break;
1047
1048     case kCharm: 
1049       if (sname=="FLAT"){
1050         func= PtCharmFlat;
1051         break;
1052       }
1053
1054       if (sname=="MUON"){
1055         func= PtCharmMUON;
1056         break;
1057       }
1058
1059       if (sname=="GSI"){
1060         func= PtCharmGSI;
1061         break;
1062       }
1063       func=0;
1064       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1065       break;
1066
1067     case kBeauty: 
1068       if (sname=="FLAT"){
1069         func= PtBeautyFlat;
1070         break;
1071       }
1072       if (sname=="MUON"){
1073         func= PtBeautyMUON;
1074         break;
1075       }
1076       if (sname=="GSI"){
1077         func= PtBeautyGSI;
1078         break;
1079       }
1080       func=0;
1081       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1082       break;
1083
1084
1085     case kEta:
1086       func=PtEtaPHOS;
1087       break;
1088
1089     case kEtaprime:
1090       func=PtEtaprimePHOS;
1091       break;
1092
1093     case kOmega:
1094       func=PtOmega;
1095       break;
1096
1097     case kRho:
1098       func=PtRho;
1099       break;
1100
1101     case kKaon:
1102       func=PtKaonPHOS;
1103       break;
1104
1105     case kPion:
1106       func=PtPion;
1107       break;
1108
1109     case kPhi:
1110       func=PtPhiPHOS;
1111       break;
1112
1113          //    case kLambda:
1114          //         func=PtLambda;
1115          //         break;
1116
1117     case kBaryons:
1118       func=PtBaryons;
1119       break;
1120
1121     default:
1122       func=0;
1123       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1124     }
1125    return func;
1126 }
1127
1128 GenFunc AliGenGSIlib::GetY(Int_t param, const char * tname) const
1129 {
1130 // Return pointer to y- parameterisation
1131    GenFunc func=0;
1132    TString sname(tname);
1133
1134    switch (param) 
1135     {
1136     case kUpsilon:
1137       if (sname=="FLAT"){
1138         func= YUpsilonFlat;
1139         break;
1140       }
1141       if (sname=="MUON"){
1142         func= YUpsilonMUON;
1143         break;
1144       }
1145       if (sname=="RITMAN"){
1146         func=YUpsilonRitman;
1147         break;
1148       }
1149       if (sname=="KAREL"){
1150         func=YUpsilonKarel;
1151         break;
1152       }
1153       func=0;
1154       printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
1155       break;
1156
1157     case kJPsi:
1158       if (sname=="FLAT"){
1159         func= YJpsiFlat;
1160         break;
1161       }
1162       if (sname=="MUON"){
1163         func= YJpsiMUON;
1164         break;
1165       }
1166       //      if (sname=="SERGEI"){
1167       //        func= YJpsi;
1168       //        break;
1169       //      }
1170    
1171       func=0;
1172       printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
1173       break;
1174      
1175     case kCharm: 
1176         func= YCharm;
1177         break;
1178
1179     case kBeauty: 
1180         func= YBeauty;
1181         break;
1182
1183     case kEta:
1184          func=YEtaPHOS;
1185          break;
1186
1187     case kEtaprime:
1188          func=YEtaprimePHOS;
1189          break;
1190
1191     case kOmega:
1192          func=YOmega;
1193          break;
1194
1195     case kRho:
1196          func=YRho;
1197          break;
1198
1199     case kKaon:
1200          func=YKaonPHOS;
1201          break;
1202
1203     case kPion:
1204          func=YPion;
1205          break;
1206
1207     case kPhi:
1208          func=YPhiPHOS;
1209          break;
1210
1211          //    case kLambda:
1212          //         func=YLambda;
1213          //         break;
1214
1215     case kBaryons:
1216          func=YBaryons;
1217          break;
1218
1219     default:
1220         func=0;
1221         printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
1222     }
1223     return func;
1224 }
1225
1226 GenFuncIp AliGenGSIlib::GetIp(Int_t param, const char * tname) const
1227 {
1228 // Return pointer to particle type parameterisation
1229    GenFuncIp func=0;
1230    TString sname(tname);
1231
1232    switch (param) 
1233     {
1234     case kUpsilon:
1235         func= IpUpsilon;
1236         break;
1237
1238     case kJPsi:
1239         func= IpJpsi;
1240         break;
1241
1242     case kCharm: 
1243         func= IpCharm;
1244         break;
1245
1246     case kBeauty: 
1247         func= IpBeauty;
1248         break;
1249
1250     case kEta:
1251          func=IpEta;
1252          break;
1253
1254     case kEtaprime:
1255          func=IpEtaprime;
1256          break;
1257
1258     case kOmega:
1259          func=IpOmega;
1260          break;
1261
1262     case kRho:
1263          func=IpRho;
1264          break;
1265
1266     case kKaon:
1267          func=IpKaonPHOS;
1268          break;
1269
1270     case kPion:
1271          func=IpPionPHOS;
1272          break;
1273
1274     case kPhi:
1275          func=IpPhi;
1276          break;
1277
1278          //    case kLambda:
1279          //         func=IpLambda;
1280          //         break;
1281
1282     case kBaryons:
1283          func=IpBaryons;
1284          break;
1285
1286     default:
1287         func=0;
1288         printf("<AliGenGSIlib::GetIp> unknown parametrisation\n");
1289     }
1290     return func;
1291 }
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303