newly implemented model for slow nucleon production
[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,2);  //  2==> 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,4);  //  4==> 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,3);  //  3==> 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,10);  //  10==> 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,1);  //  1==> 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,5);  //  5==> 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,6);  //  6==> 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   0=>PI, 1=>K, 2=>ETA, 3=>OMEGA, 4=>ETA', 5=>PHI 
973 //        6=>BARYON-BARYONBAR, 10==>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   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
983   Double_t kfmax2=f5/kfmax[np];
984   // PIONS
985   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
986   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
987                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
988   return fmtscal*ptpion;
989 }
990
991 //==========================================================================
992 //
993 //                     Set Getters 
994 //    
995 //==========================================================================
996
997 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
998
999 typedef Int_t (*GenFuncIp) (TRandom *);
1000
1001 GenFunc AliGenGSIlib::GetPt(Int_t param, const char * tname) const
1002 {
1003 // Return pointer to pT parameterisation
1004    GenFunc func=0;
1005    TString sname(tname);
1006
1007    switch (param) 
1008     {
1009     case kUpsilon:
1010       if (sname=="FLAT"){
1011         func= PtUpsilonFlat;
1012         break;
1013       }
1014       if (sname=="MUON"){
1015         func= PtUpsilonMUON;
1016         break;
1017       }
1018       if (sname=="RITMAN"){
1019         func=PtUpsilonRitman;
1020         break;
1021       }
1022       if (sname=="KAREL"){
1023         func=PtUpsilonKarel;
1024         break;
1025       }
1026       func=0;
1027       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1028       break;
1029
1030     case kJPsi:
1031       if (sname=="FLAT"){
1032         func= PtJpsiFlat;
1033         break;
1034       }
1035       if (sname=="MUON"){
1036         func= PtJpsiMUON;
1037         break;
1038       }
1039       //      if (sname=="SERGEI"){
1040       //        func= PtJpsi;
1041       //        break;
1042       //     }
1043       func=0;
1044       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1045       break;
1046
1047     case kCharm: 
1048       if (sname=="FLAT"){
1049         func= PtCharmFlat;
1050         break;
1051       }
1052
1053       if (sname=="MUON"){
1054         func= PtCharmMUON;
1055         break;
1056       }
1057
1058       if (sname=="GSI"){
1059         func= PtCharmGSI;
1060         break;
1061       }
1062       func=0;
1063       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1064       break;
1065
1066     case kBeauty: 
1067       if (sname=="FLAT"){
1068         func= PtBeautyFlat;
1069         break;
1070       }
1071       if (sname=="MUON"){
1072         func= PtBeautyMUON;
1073         break;
1074       }
1075       if (sname=="GSI"){
1076         func= PtBeautyGSI;
1077         break;
1078       }
1079       func=0;
1080       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1081       break;
1082
1083
1084     case kEta:
1085       func=PtEtaPHOS;
1086       break;
1087
1088     case kEtaprime:
1089       func=PtEtaprimePHOS;
1090       break;
1091
1092     case kOmega:
1093       func=PtOmega;
1094       break;
1095
1096     case kRho:
1097       func=PtRho;
1098       break;
1099
1100     case kKaon:
1101       func=PtKaonPHOS;
1102       break;
1103
1104     case kPion:
1105       func=PtPion;
1106       break;
1107
1108     case kPhi:
1109       func=PtPhiPHOS;
1110       break;
1111
1112          //    case kLambda:
1113          //         func=PtLambda;
1114          //         break;
1115
1116     case kBaryons:
1117       func=PtBaryons;
1118       break;
1119
1120     default:
1121       func=0;
1122       printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1123     }
1124    return func;
1125 }
1126
1127 GenFunc AliGenGSIlib::GetY(Int_t param, const char * tname) const
1128 {
1129 // Return pointer to y- parameterisation
1130    GenFunc func=0;
1131    TString sname(tname);
1132
1133    switch (param) 
1134     {
1135     case kUpsilon:
1136       if (sname=="FLAT"){
1137         func= YUpsilonFlat;
1138         break;
1139       }
1140       if (sname=="MUON"){
1141         func= YUpsilonMUON;
1142         break;
1143       }
1144       if (sname=="RITMAN"){
1145         func=YUpsilonRitman;
1146         break;
1147       }
1148       if (sname=="KAREL"){
1149         func=YUpsilonKarel;
1150         break;
1151       }
1152       func=0;
1153       printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
1154       break;
1155
1156     case kJPsi:
1157       if (sname=="FLAT"){
1158         func= YJpsiFlat;
1159         break;
1160       }
1161       if (sname=="MUON"){
1162         func= YJpsiMUON;
1163         break;
1164       }
1165       //      if (sname=="SERGEI"){
1166       //        func= YJpsi;
1167       //        break;
1168       //      }
1169    
1170       func=0;
1171       printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
1172       break;
1173      
1174     case kCharm: 
1175         func= YCharm;
1176         break;
1177
1178     case kBeauty: 
1179         func= YBeauty;
1180         break;
1181
1182     case kEta:
1183          func=YEtaPHOS;
1184          break;
1185
1186     case kEtaprime:
1187          func=YEtaprimePHOS;
1188          break;
1189
1190     case kOmega:
1191          func=YOmega;
1192          break;
1193
1194     case kRho:
1195          func=YRho;
1196          break;
1197
1198     case kKaon:
1199          func=YKaonPHOS;
1200          break;
1201
1202     case kPion:
1203          func=YPion;
1204          break;
1205
1206     case kPhi:
1207          func=YPhiPHOS;
1208          break;
1209
1210          //    case kLambda:
1211          //         func=YLambda;
1212          //         break;
1213
1214     case kBaryons:
1215          func=YBaryons;
1216          break;
1217
1218     default:
1219         func=0;
1220         printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
1221     }
1222     return func;
1223 }
1224
1225 GenFuncIp AliGenGSIlib::GetIp(Int_t param, const char * tname) const
1226 {
1227 // Return pointer to particle type parameterisation
1228    GenFuncIp func=0;
1229    TString sname(tname);
1230
1231    switch (param) 
1232     {
1233     case kUpsilon:
1234         func= IpUpsilon;
1235         break;
1236
1237     case kJPsi:
1238         func= IpJpsi;
1239         break;
1240
1241     case kCharm: 
1242         func= IpCharm;
1243         break;
1244
1245     case kBeauty: 
1246         func= IpBeauty;
1247         break;
1248
1249     case kEta:
1250          func=IpEta;
1251          break;
1252
1253     case kEtaprime:
1254          func=IpEtaprime;
1255          break;
1256
1257     case kOmega:
1258          func=IpOmega;
1259          break;
1260
1261     case kRho:
1262          func=IpRho;
1263          break;
1264
1265     case kKaon:
1266          func=IpKaonPHOS;
1267          break;
1268
1269     case kPion:
1270          func=IpPionPHOS;
1271          break;
1272
1273     case kPhi:
1274          func=IpPhi;
1275          break;
1276
1277          //    case kLambda:
1278          //         func=IpLambda;
1279          //         break;
1280
1281     case kBaryons:
1282          func=IpBaryons;
1283          break;
1284
1285     default:
1286         func=0;
1287         printf("<AliGenGSIlib::GetIp> unknown parametrisation\n");
1288     }
1289     return func;
1290 }
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302