Overlaps corrected, new shape of sectors
[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 //
19 // Library of generators for PMD
20 // providing y and pt parametrisation
21 // for generated tracks
22 // Specific for PMD needs
23 // Author: PMD Offline Group
24 //
25
26 #include <TMath.h>
27 #include <TPDGCode.h>
28
29 #include "AliGenPMDlib.h"
30
31 ClassImp(AliGenPMDlib)
32 //
33 //  Neutral Pions
34
35 Double_t AliGenPMDlib::PtPi0(const Double_t *px, const 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 //
63 // y-distribution
64 //
65 Double_t AliGenPMDlib::YPi0( const Double_t *py, const Double_t */*dummy*/)
66 {
67   //
68   // y parametrisation for pi0
69   //
70     const Double_t ka1    = 4913.;
71     const Double_t ka2    = 1819.;
72     const Double_t keta1  = 0.22;
73     const Double_t keta2  = 3.66;
74     const Double_t kdeta1 = 1.47;
75     const Double_t kdeta2 = 1.51;
76     Double_t y=TMath::Abs(*py);
77     //
78     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
79     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
80     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
81 }
82
83 //                 particle composition
84 //
85 Int_t AliGenPMDlib::IpPi0(TRandom *)
86 {
87 // Pi0
88     return kPi0;
89 }
90
91 //____________________________________________________________
92 //
93 // Mt-scaling
94
95 Double_t AliGenPMDlib::PtScal(Double_t pt, Int_t np)
96 {
97   //    SCALING EN MASSE PAR RAPPORT A PTPI
98   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
99   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
100   //     VALUE MESON/PI AT 5 GEV
101   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
102   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
103   Double_t fmax2=f5/kfmax[np];
104   // PIONS
105   Double_t ptpion=100.*PtPi0(&pt, (Double_t*) 0);
106   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
107                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
108   return fmtscal*ptpion;
109 }
110 //
111 // kaon
112 //
113 //                pt-distribution
114 //____________________________________________________________
115
116 Double_t AliGenPMDlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
117 {
118 // Kaon pT
119   return PtScal(*px, 1);
120 }
121
122 // y-distribution
123 //____________________________________________________________
124 Double_t AliGenPMDlib::YEta( const Double_t *py, const Double_t */*dummy*/)
125 {
126     //
127     // y parametrisation for etas
128     //
129     const Double_t ka1    = 4913.;
130     const Double_t ka2    = 1819.;
131     const Double_t keta1  = 0.22;
132     const Double_t keta2  = 3.66;
133     const Double_t kdeta1 = 1.47;
134     const Double_t kdeta2 = 1.51;
135     Double_t y=TMath::Abs(*py);
136     //
137     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
138     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
139     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
140 }
141
142 //                 particle composition
143 //
144 Int_t AliGenPMDlib::IpEta(TRandom *)
145 {
146     return 221;
147 }
148
149
150 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
151 GenFunc AliGenPMDlib::GetPt(Int_t param,  const char* /*tname*/) const
152 {
153 // Return pointer to pT parameterisation
154     GenFunc func=NULL;
155     switch (param) 
156     {
157     case kPion:
158         func=PtPi0;
159         break;
160     case kEta:
161         func=PtEta;
162         break;
163     default:
164         func=0;
165         printf("<AliGenPMDlib::GetPt> unknown parametrisation\n");
166     }
167     return func;
168 }
169
170 GenFunc AliGenPMDlib::GetY(Int_t param, const char* /*tname*/) const
171 {
172 // Return pointer to y- parameterisation
173     GenFunc func=NULL;
174     switch (param) 
175     {
176     case kPion:
177         func=YPi0;
178         break;
179     case kEta:
180         func=YEta;
181         break;
182     default:
183         func=0;
184         printf("<AliGenPMDlib::GetY> unknown parametrisation\n");
185     }
186     return func;
187
188 }
189 typedef Int_t (*GenFuncIp) (TRandom *);
190 GenFuncIp AliGenPMDlib::GetIp(Int_t param,  const char* /*tname*/) const
191 {
192 // Return pointer to particle type parameterisation
193     GenFuncIp func=NULL;
194     switch (param) 
195     {
196     case kPion:
197         func=IpPi0;
198         break;
199     case kEta:
200         func=IpEta;
201         break;
202     default:
203         printf("<AliGenPMDlib::GetIp> unknown parametrisation\n");
204     }
205     return func;
206 }
207
208
209
210