Coding convention corrections + few minor bug fixes
[u/mrichter/AliRoot.git] / EVGEN / AliGenHIJINGpara.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 /*
17 $Log$
18 Revision 1.1  2000/06/09 20:20:30  morsch
19 Same class as previously in AliSimpleGen.cxx
20 All coding rule violations except RS3 corrected (AM)
21
22 */
23 ///////////////////////////////////////////////////////////////////
24 //                                                               //
25 //    Generate the final state of the interaction as the input   //
26 //    to the MonteCarlo                                          //
27 //
28 //Begin_Html
29 /*
30 <img src="picts/AliGeneratorClass.gif">
31 </pre>
32 <br clear=left>
33 <font size=+2 color=red>
34 <p>The responsible person for this module is
35 <a href="mailto:andreas.morsch@cern.ch">Andreas Morsch</a>.
36 </font>
37 <pre>
38 */
39 //End_Html
40 //                                                               //
41 ///////////////////////////////////////////////////////////////////
42
43 #include "AliGenHIJINGpara.h"
44 #include "AliRun.h"
45 #include "AliConst.h"
46 #include "AliPDG.h"
47
48 ClassImp(AliGenHIJINGpara)
49
50 AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para)
51 {
52 // copy constructor
53 }
54
55 //_____________________________________________________________________________
56 static Double_t ptpi(Double_t *px, Double_t *)
57 {
58   //
59   //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
60   //     POWER LAW FOR PT > 500 MEV
61   //     MT SCALING BELOW (T=160 MEV)
62   //
63   const Double_t kp0 = 1.3;
64   const Double_t kxn = 8.28;
65   const Double_t kxlim=0.5;
66   const Double_t kt=0.160;
67   const Double_t kxmpi=0.139;
68   const Double_t kb=1.;
69   Double_t y, y1, xmpi2, ynorm, a;
70   Double_t x=*px;
71   //
72   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
73   xmpi2=kxmpi*kxmpi;
74   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
75   a=ynorm/y1;
76   if (x > kxlim)
77     y=a*TMath::Power(kp0/(kp0+x),kxn);
78   else
79     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
80   return y*x;
81 }
82
83 //_____________________________________________________________________________
84 static Double_t ptscal(Double_t pt, Int_t np)
85 {
86     //    SCALING EN MASSE PAR RAPPORT A PTPI
87     //     MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
88     const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
89     //     VALUE MESON/PI AT 5 GEV
90     const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
91     np--;
92     Double_t f5=TMath::Power(((
93         sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
94     Double_t fmax2=f5/kfmax[np];
95     // PIONS
96     Double_t ptpion=100.*ptpi(&pt, (Double_t*) 0);
97     Double_t fmtscal=TMath::Power(((
98         sqrt(pt*pt+0.018215)+2.)/ (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ 
99         fmax2;
100     return fmtscal*ptpion;
101 }
102
103 //_____________________________________________________________________________
104 static Double_t ptka( Double_t *px, Double_t *)
105 {
106     //
107     // pt parametrisation for k
108     //
109     return ptscal(*px,2);
110 }
111
112
113 //_____________________________________________________________________________
114 static Double_t etapic( Double_t *py, Double_t *)
115 {
116   //
117   // eta parametrisation for pi
118   //
119     const Double_t ka1    = 4913.;
120     const Double_t ka2    = 1819.;
121     const Double_t keta1  = 0.22;
122     const Double_t keta2  = 3.66;
123     const Double_t kdeta1 = 1.47;
124     const Double_t kdeta2 = 1.51;
125     Double_t y=TMath::Abs(*py);
126     //
127     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
128     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
129     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
130 }
131
132 //_____________________________________________________________________________
133 static Double_t etakac( Double_t *py, Double_t *)
134 {
135     //
136     // eta parametrisation for ka
137     //
138     const Double_t ka1    = 497.6;
139     const Double_t ka2    = 215.6;
140     const Double_t keta1  = 0.79;
141     const Double_t keta2  = 4.09;
142     const Double_t kdeta1 = 1.54;
143     const Double_t kdeta2 = 1.40;
144     Double_t y=TMath::Abs(*py);
145     //
146     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
147     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
148     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
149 }
150
151 //_____________________________________________________________________________
152 AliGenHIJINGpara::AliGenHIJINGpara()
153   :AliGenerator()
154 {
155     //
156     // Default constructor
157     //
158     fPtpi = 0;
159     fPtka = 0;
160     fETApic = 0;
161     fETAkac = 0;
162 }
163
164 //_____________________________________________________________________________
165 AliGenHIJINGpara::AliGenHIJINGpara(Int_t npart)
166   :AliGenerator(npart)
167 {
168   // 
169   // Standard constructor
170   //
171     fName="HIGINGpara";
172     fTitle="HIJING Parametrisation Particle Generator";
173     fPtpi = 0;
174     fPtka = 0;
175     fETApic = 0;
176     fETAkac = 0;
177 }
178
179 //_____________________________________________________________________________
180 AliGenHIJINGpara::~AliGenHIJINGpara()
181 {
182   //
183   // Standard destructor
184   //
185     delete fPtpi;
186     delete fPtka;
187     delete fETApic;
188     delete fETAkac;
189 }
190
191 //_____________________________________________________________________________
192 void AliGenHIJINGpara::Init()
193 {
194   //
195   // Initialise the HIJING parametrisation
196   //
197     Float_t etaMin =-TMath::Log(TMath::Tan(
198         TMath::Min((Double_t)fThetaMax/2,TMath::Pi()/2-1.e-10)));
199     Float_t etaMax = -TMath::Log(TMath::Tan(
200         TMath::Max((Double_t)fThetaMin/2,1.e-10)));
201     fPtpi = new TF1("ptpi",&ptpi,0,20,0);
202     fPtka = new TF1("ptka",&ptka,0,20,0);
203     fETApic = new TF1("etapic",&etapic,etaMin,etaMax,0);
204     fETAkac = new TF1("etakac",&etakac,etaMin,etaMax,0);
205     TF1 *etaPic0 = new TF1("etapic",&etapic,-7,7,0);
206     TF1 *etaKac0 = new TF1("etakac",&etakac,-7,7,0);
207     Float_t intETApi  = etaPic0->Integral(-0.5, 0.5);
208     Float_t intETAka  = etaKac0->Integral(-0.5, 0.5);
209     Float_t scalePi=7316/(intETApi/1.5);
210     Float_t scaleKa= 684/(intETAka/2.0);
211     
212     Float_t intPt  = (0.877*etaPic0->Integral(0, 15)+
213                       0.123*etaKac0->Integral(0, 15));
214     Float_t intPtSel = (0.877*etaPic0->Integral(fPtMin, fPtMax)+
215                         0.123*etaKac0->Integral(fPtMin, fPtMax));
216     Float_t ptFrac = intPtSel/intPt;
217     
218     
219     Float_t intETASel  = (scalePi*etaPic0->Integral(etaMin, etaMax)+
220                           scaleKa*etaKac0->Integral(etaMin, etaMax));
221     Float_t phiFrac = (fPhiMax-fPhiMin)/2/TMath::Pi();
222     fParentWeight = Float_t(fNpart)/intETASel*ptFrac*phiFrac;
223     
224     printf("\n The number of particles in the selected kinematic region corresponds to %f percent of a full event\n ", 100.*fParentWeight);
225     
226 }
227
228 //_____________________________________________________________________________
229 void AliGenHIJINGpara::Generate()
230 {
231   //
232   // Generate one trigger
233   //
234
235   
236     const Float_t kRaKpic=0.14;
237     const Float_t kBorne=1/(1+kRaKpic);
238     Float_t polar[3]= {0,0,0};
239     //
240     const Int_t kPions[3] = {kPi0, kPiPlus, kPiMinus};
241     const Int_t kKaons[4] = {kK0Long, kK0Short, kKPlus, kKMinus};
242     //
243     Float_t origin[3];
244     Float_t pt, pl, ptot;
245     Float_t phi, theta;
246     Float_t p[3];
247     Int_t i, part, nt, j;
248     //
249     TF1 *ptf;
250     TF1 *etaf;
251     //
252     Float_t random[6];
253     //
254     for (j=0;j<3;j++) origin[j]=fOrigin[j];
255     if(fVertexSmear==kPerEvent) {
256         gMC->Rndm(random,6);
257         for (j=0;j<3;j++) {
258             origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
259                 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
260         }
261     }
262     for(i=0;i<fNpart;i++) {
263         while(1) {
264             gMC->Rndm(random,3);
265             if(random[0]<kBorne) {
266                 part=kPions[Int_t (random[1]*3)];
267                 ptf=fPtpi;
268               etaf=fETApic;
269             } else {
270                 part=kKaons[Int_t (random[1]*4)];
271                 ptf=fPtka;
272                 etaf=fETAkac;
273             }
274             phi=fPhiMin+random[2]*(fPhiMax-fPhiMin);
275             theta=2*TMath::ATan(TMath::Exp(-etaf->GetRandom()));
276             if(theta<fThetaMin || theta>fThetaMax) continue;
277             pt=ptf->GetRandom();
278             pl=pt/TMath::Tan(theta);
279             ptot=TMath::Sqrt(pt*pt+pl*pl);
280             if(ptot<fPMin || ptot>fPMax) continue;
281             p[0]=pt*TMath::Cos(phi);
282             p[1]=pt*TMath::Sin(phi);
283             p[2]=pl;
284             if(fVertexSmear==kPerTrack) {
285                 gMC->Rndm(random,6);
286                 for (j=0;j<3;j++) {
287                     origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
288                         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
289                 }
290             }
291             gAlice->SetTrack(fTrackIt,-1,part,p,origin,polar,0,"Primary",nt,fParentWeight);
292             break;
293         }
294     }
295 }
296
297 AliGenHIJINGpara& AliGenHIJINGpara::operator=(const  AliGenHIJINGpara& rhs)
298 {
299 // Assignment operator
300     return *this;
301 }
302
303