03a25f9b45265abe170f2974d65aa76b24be4940
[u/mrichter/AliRoot.git] / EVGEN / AliGenPMDlib.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 #include <TMath.h>
19 #include <TPDGCode.h>
20
21 #include "AliGenPMDlib.h"
22
23 ClassImp(AliGenPMDlib)
24 //
25 //  Neutral Pions
26
27 Double_t AliGenPMDlib::PtPi0(Double_t *px, Double_t *dummy)
28 {
29 //
30 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
31 //     POWER LAW FOR PT > 500 MEV
32 //     MT SCALING BELOW (T=160 MEV)
33 //
34   const Double_t kp0 = 1.3;
35   const Double_t kxn = 8.28;
36   const Double_t kxlim=0.5;
37   const Double_t kt=0.160;
38   const Double_t kxmpi=0.139;
39   const Double_t kb=1.;
40   Double_t y, y1, xmpi2, ynorm, a;
41   Double_t x=*px;
42   //
43   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
44   xmpi2=kxmpi*kxmpi;
45   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
46   a=ynorm/y1;
47   if (x > kxlim)
48     y=a*TMath::Power(kp0/(kp0+x),kxn);
49   else
50     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
51   return y*x;
52 }
53
54 //
55 // y-distribution
56 //
57 Double_t AliGenPMDlib::YPi0( Double_t *py, Double_t *dummy)
58 {
59   //
60   // y parametrisation for pi0
61   //
62     const Double_t ka1    = 4913.;
63     const Double_t ka2    = 1819.;
64     const Double_t keta1  = 0.22;
65     const Double_t keta2  = 3.66;
66     const Double_t kdeta1 = 1.47;
67     const Double_t kdeta2 = 1.51;
68     Double_t y=TMath::Abs(*py);
69     //
70     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
71     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
72     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
73 }
74
75 //                 particle composition
76 //
77 Int_t AliGenPMDlib::IpPi0(TRandom *)
78 {
79 // Pi0
80     return kPi0;
81 }
82
83 //____________________________________________________________
84 //
85 // Mt-scaling
86
87 Double_t AliGenPMDlib::PtScal(Double_t pt, Int_t np)
88 {
89   //    SCALING EN MASSE PAR RAPPORT A PTPI
90   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
91   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
92   //     VALUE MESON/PI AT 5 GEV
93   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
94   np--;
95   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
96   Double_t fmax2=f5/kfmax[np];
97   // PIONS
98   Double_t ptpion=100.*PtPi0(&pt, (Double_t*) 0);
99   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
100                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
101   return fmtscal*ptpion;
102 }
103 //
104 // kaon
105 //
106 //                pt-distribution
107 //____________________________________________________________
108
109 Double_t AliGenPMDlib::PtEta( Double_t *px, Double_t *dummy)
110 {
111 // Kaon pT
112   return PtScal(*px,3);
113 }
114
115 // y-distribution
116 //____________________________________________________________
117 Double_t AliGenPMDlib::YEta( Double_t *py, Double_t *dummy)
118 {
119     //
120     // y parametrisation for etas
121     //
122     const Double_t ka1    = 4913.;
123     const Double_t ka2    = 1819.;
124     const Double_t keta1  = 0.22;
125     const Double_t keta2  = 3.66;
126     const Double_t kdeta1 = 1.47;
127     const Double_t kdeta2 = 1.51;
128     Double_t y=TMath::Abs(*py);
129     //
130     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
131     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
132     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
133 }
134
135 //                 particle composition
136 //
137 Int_t AliGenPMDlib::IpEta(TRandom *)
138 {
139     return 221;
140 }
141
142
143 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
144 GenFunc AliGenPMDlib::GetPt(Int_t param,  const char* tname) const
145 {
146 // Return pointer to pT parameterisation
147     GenFunc func=NULL;
148     switch (param) 
149     {
150     case kPion:
151         func=PtPi0;
152         break;
153     case kEta:
154         func=PtEta;
155         break;
156     default:
157         func=0;
158         printf("<AliGenPMDlib::GetPt> unknown parametrisation\n");
159     }
160     return func;
161 }
162
163 GenFunc AliGenPMDlib::GetY(Int_t param, const char* tname) const
164 {
165 // Return pointer to y- parameterisation
166     GenFunc func=NULL;
167     switch (param) 
168     {
169     case kPion:
170         func=YPi0;
171         break;
172     case kEta:
173         func=YEta;
174         break;
175     default:
176         func=0;
177         printf("<AliGenPMDlib::GetY> unknown parametrisation\n");
178     }
179     return func;
180
181 }
182 typedef Int_t (*GenFuncIp) (TRandom *);
183 GenFuncIp AliGenPMDlib::GetIp(Int_t param,  const char* tname) const
184 {
185 // Return pointer to particle type parameterisation
186     GenFuncIp func=NULL;
187     switch (param) 
188     {
189     case kPion:
190         func=IpPi0;
191         break;
192     case kEta:
193         func=IpEta;
194         break;
195     default:
196         printf("<AliGenPMDlib::GetIp> unknown parametrisation\n");
197     }
198     return func;
199 }
200
201
202
203