Merging the VirtualMC branch to the main development branch (HEAD)
[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.11.4.1  2002/07/24 08:56:28  alibrary
19 Updating EVGEN on TVirtulaMC
20
21 Revision 1.12  2002/06/19 06:56:34  hristov
22 Memory leak corrected
23
24 Revision 1.11  2002/03/20 10:21:13  hristov
25 Set fPtMax to 15 GeV in order to avoid some numerical problems
26
27 Revision 1.10  2001/10/15 16:44:46  morsch
28 - Possibility for vertex distribution truncation.
29 - Write mc header with vertex position.
30
31 Revision 1.9  2001/07/27 17:09:36  morsch
32 Use local SetTrack, KeepTrack and SetHighWaterMark methods
33 to delegate either to local stack or to stack owned by AliRun.
34 (Piotr Skowronski, A.M.)
35
36 Revision 1.8  2001/07/20 11:03:58  morsch
37 Issue warning message if used outside allowed eta range (-8 to 8).
38
39 Revision 1.7  2001/07/17 12:41:01  morsch
40 - Calculation of fraction of event corresponding to selected pt-range corrected
41 (R. Turrisi)
42 - Parent weight corrected.
43
44 Revision 1.6  2001/05/16 14:57:10  alibrary
45 New files for folders and Stack
46
47 Revision 1.5  2000/12/21 16:24:06  morsch
48 Coding convention clean-up
49
50 Revision 1.4  2000/11/30 07:12:50  alibrary
51 Introducing new Rndm and QA classes
52
53 Revision 1.3  2000/10/02 21:28:06  fca
54 Removal of useless dependecies via forward declarations
55
56 Revision 1.2  2000/07/11 18:24:55  fca
57 Coding convention corrections + few minor bug fixes
58
59 Revision 1.1  2000/06/09 20:20:30  morsch
60 Same class as previously in AliSimpleGen.cxx
61 All coding rule violations except RS3 corrected (AM)
62
63 */
64
65 // Parameterisation of pi and K, eta and pt distributions
66 // used for the ALICE TDRs.
67 // eta: according to HIJING (shadowing + quenching)
68 // pT : according to CDF measurement at 1.8 TeV
69 // Author: andreas.morsch@cern.ch
70
71
72 //Begin_Html
73 /*
74 <img src="picts/AliGeneratorClass.gif">
75 </pre>
76 <br clear=left>
77 <font size=+2 color=red>
78 <p>The responsible person for this module is
79 <a href="mailto:andreas.morsch@cern.ch">Andreas Morsch</a>.
80 </font>
81 <pre>
82 */
83 //End_Html
84 //                                                               //
85 ///////////////////////////////////////////////////////////////////
86
87 #include "AliGenHIJINGpara.h"
88 #include "AliGenEventHeader.h"
89 #include "AliRun.h"
90 #include "AliConst.h"
91 #include "AliPDG.h"
92
93 #include <TF1.h>
94 #include <TArrayF.h>
95
96 ClassImp(AliGenHIJINGpara)
97
98 AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para)
99 {
100 // copy constructor
101 }
102
103 //_____________________________________________________________________________
104 static Double_t ptpi(Double_t *px, Double_t *)
105 {
106   //
107   //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
108   //     POWER LAW FOR PT > 500 MEV
109   //     MT SCALING BELOW (T=160 MEV)
110   //
111   const Double_t kp0 = 1.3;
112   const Double_t kxn = 8.28;
113   const Double_t kxlim=0.5;
114   const Double_t kt=0.160;
115   const Double_t kxmpi=0.139;
116   const Double_t kb=1.;
117   Double_t y, y1, xmpi2, ynorm, a;
118   Double_t x=*px;
119   //
120   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
121   xmpi2=kxmpi*kxmpi;
122   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
123   a=ynorm/y1;
124   if (x > kxlim)
125     y=a*TMath::Power(kp0/(kp0+x),kxn);
126   else
127     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
128   return y*x;
129 }
130
131 //_____________________________________________________________________________
132 static Double_t ptscal(Double_t pt, Int_t np)
133 {
134     //    SCALING EN MASSE PAR RAPPORT A PTPI
135     //     MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
136     const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
137     //     VALUE MESON/PI AT 5 GEV
138     const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
139     np--;
140     Double_t f5=TMath::Power(((
141         sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
142     Double_t fmax2=f5/kfmax[np];
143     // PIONS
144     Double_t ptpion=100.*ptpi(&pt, (Double_t*) 0);
145     Double_t fmtscal=TMath::Power(((
146         sqrt(pt*pt+0.018215)+2.)/ (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ 
147         fmax2;
148     return fmtscal*ptpion;
149 }
150
151 //_____________________________________________________________________________
152 static Double_t ptka( Double_t *px, Double_t *)
153 {
154     //
155     // pt parametrisation for k
156     //
157     return ptscal(*px,2);
158 }
159
160
161 //_____________________________________________________________________________
162 static Double_t etapic( Double_t *py, Double_t *)
163 {
164   //
165   // eta parametrisation for pi
166   //
167     const Double_t ka1    = 4913.;
168     const Double_t ka2    = 1819.;
169     const Double_t keta1  = 0.22;
170     const Double_t keta2  = 3.66;
171     const Double_t kdeta1 = 1.47;
172     const Double_t kdeta2 = 1.51;
173     Double_t y=TMath::Abs(*py);
174     //
175     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
176     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
177     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
178 }
179
180 //_____________________________________________________________________________
181 static Double_t etakac( Double_t *py, Double_t *)
182 {
183     //
184     // eta parametrisation for ka
185     //
186     const Double_t ka1    = 497.6;
187     const Double_t ka2    = 215.6;
188     const Double_t keta1  = 0.79;
189     const Double_t keta2  = 4.09;
190     const Double_t kdeta1 = 1.54;
191     const Double_t kdeta2 = 1.40;
192     Double_t y=TMath::Abs(*py);
193     //
194     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
195     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
196     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
197 }
198
199 //_____________________________________________________________________________
200 AliGenHIJINGpara::AliGenHIJINGpara()
201   :AliGenerator()
202 {
203     //
204     // Default constructor
205     //
206     fPtpi = 0;
207     fPtka = 0;
208     fETApic = 0;
209     fETAkac = 0;
210     SetCutVertexZ();
211     SetPtRange();
212 }
213
214 //_____________________________________________________________________________
215 AliGenHIJINGpara::AliGenHIJINGpara(Int_t npart)
216   :AliGenerator(npart)
217 {
218   // 
219   // Standard constructor
220   //
221     fName="HIGINGpara";
222     fTitle="HIJING Parametrisation Particle Generator";
223     fPtpi = 0;
224     fPtka = 0;
225     fETApic = 0;
226     fETAkac = 0;
227     SetCutVertexZ();
228     SetPtRange();
229 }
230
231 //_____________________________________________________________________________
232 AliGenHIJINGpara::~AliGenHIJINGpara()
233 {
234   //
235   // Standard destructor
236   //
237     delete fPtpi;
238     delete fPtka;
239     delete fETApic;
240     delete fETAkac;
241 }
242
243 //_____________________________________________________________________________
244 void AliGenHIJINGpara::Init()
245 {
246   //
247   // Initialise the HIJING parametrisation
248   //
249     Float_t etaMin =-TMath::Log(TMath::Tan(
250         TMath::Min((Double_t)fThetaMax/2,TMath::Pi()/2-1.e-10)));
251     Float_t etaMax = -TMath::Log(TMath::Tan(
252         TMath::Max((Double_t)fThetaMin/2,1.e-10)));
253     fPtpi   = new TF1("ptpi",&ptpi,0,20,0);
254     fPtka   = new TF1("ptka",&ptka,0,20,0);
255     fETApic = new TF1("etapic",&etapic,etaMin,etaMax,0);
256     fETAkac = new TF1("etakac",&etakac,etaMin,etaMax,0);
257
258     TF1 etaPic0("etapic",&etapic,-7,7,0);
259     TF1 etaKac0("etakac",&etakac,-7,7,0);
260
261     TF1 ptPic0("ptpi",&ptpi,0.,15.,0);
262     TF1 ptKac0("ptka",&ptka,0.,15.,0);
263
264     Float_t intETApi  = etaPic0.Integral(-0.5, 0.5);
265     Float_t intETAka  = etaKac0.Integral(-0.5, 0.5);
266     Float_t scalePi   = 7316/(intETApi/1.5);
267     Float_t scaleKa   =  684/(intETAka/2.0);
268
269 //  Fraction of events corresponding to the selected pt-range    
270     Float_t intPt    = (0.877*ptPic0.Integral(0, 15)+
271                         0.123*ptKac0.Integral(0, 15));
272     Float_t intPtSel = (0.877*ptPic0.Integral(fPtMin, fPtMax)+
273                         0.123*ptKac0.Integral(fPtMin, fPtMax));
274     Float_t ptFrac   = intPtSel/intPt;
275
276 //  Fraction of events corresponding to the selected eta-range    
277     Float_t intETASel  = (scalePi*etaPic0.Integral(etaMin, etaMax)+
278                           scaleKa*etaKac0.Integral(etaMin, etaMax));
279 //  Fraction of events corresponding to the selected phi-range    
280     Float_t phiFrac    = (fPhiMax-fPhiMin)/2/TMath::Pi();
281
282     fParentWeight = Float_t(fNpart)/(intETASel*ptFrac*phiFrac);
283     
284     printf("%s: The number of particles in the selected kinematic region corresponds to %f percent of a full event\n ", 
285            ClassName(),100.*fParentWeight);
286
287 // Issue warning message if etaMin or etaMax are outside the alowed range 
288 // of the parametrization
289     if (etaMin < -8.001 || etaMax > 8.001) {
290         printf("\n \n WARNING FROM AliGenHIJINGPara !");
291         printf("\n YOU ARE USING THE PARAMETERISATION OUTSIDE ");       
292         printf("\n THE ALLOWED PSEUDORAPIDITY RANGE (-8. - 8.)");           
293         printf("\n YOUR LIMITS: %f %f \n \n ", etaMin, etaMax);
294     }
295 }
296
297 //_____________________________________________________________________________
298 void AliGenHIJINGpara::Generate()
299 {
300   //
301   // Generate one trigger
302   //
303
304   
305     const Float_t kRaKpic=0.14;
306     const Float_t kBorne=1/(1+kRaKpic);
307     Float_t polar[3]= {0,0,0};
308     //
309     const Int_t kPions[3] = {kPi0, kPiPlus, kPiMinus};
310     const Int_t kKaons[4] = {kK0Long, kK0Short, kKPlus, kKMinus};
311     //
312     Float_t origin[3];
313     Float_t pt, pl, ptot;
314     Float_t phi, theta;
315     Float_t p[3];
316     Int_t i, part, nt, j;
317     //
318     TF1 *ptf;
319     TF1 *etaf;
320     //
321     Float_t random[6];
322     //
323     for (j=0;j<3;j++) origin[j]=fOrigin[j];
324
325     if(fVertexSmear == kPerEvent) {
326         Float_t dv[3];
327         dv[2] = 1.e10;
328         while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
329             Rndm(random,6);
330             for (j=0; j < 3; j++) {
331                 dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
332                     TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
333             }
334         }
335         for (j=0; j < 3; j++) origin[j] += dv[j];
336     } // if kPerEvent
337     TArrayF eventVertex;
338     eventVertex.Set(3);
339     eventVertex[0] = origin[0];
340     eventVertex[1] = origin[1];
341     eventVertex[2] = origin[2];
342
343     for(i=0;i<fNpart;i++) {
344         while(1) {
345             Rndm(random,3);
346             if(random[0]<kBorne) {
347                 part=kPions[Int_t (random[1]*3)];
348                 ptf=fPtpi;
349               etaf=fETApic;
350             } else {
351                 part=kKaons[Int_t (random[1]*4)];
352                 ptf=fPtka;
353                 etaf=fETAkac;
354             }
355             phi=fPhiMin+random[2]*(fPhiMax-fPhiMin);
356             theta=2*TMath::ATan(TMath::Exp(-etaf->GetRandom()));
357             if(theta<fThetaMin || theta>fThetaMax) continue;
358             pt=ptf->GetRandom();
359             pl=pt/TMath::Tan(theta);
360             ptot=TMath::Sqrt(pt*pt+pl*pl);
361             if(ptot<fPMin || ptot>fPMax) continue;
362             p[0]=pt*TMath::Cos(phi);
363             p[1]=pt*TMath::Sin(phi);
364             p[2]=pl;
365             if(fVertexSmear==kPerTrack) {
366                 Rndm(random,6);
367                 for (j=0;j<3;j++) {
368                     origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
369                         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
370                 }
371             }
372             SetTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,nt,fParentWeight);
373             break;
374         }
375     }
376 // Header
377     AliGenEventHeader* header = new AliGenEventHeader("HIJINGparam");
378 // Event Vertex
379     header->SetPrimaryVertex(eventVertex);
380     gAlice->SetGenEventHeader(header); 
381 }
382
383 AliGenHIJINGpara& AliGenHIJINGpara::operator=(const  AliGenHIJINGpara& rhs)
384 {
385 // Assignment operator
386     return *this;
387 }
388
389 void AliGenHIJINGpara::SetPtRange(Float_t ptmin, Float_t ptmax) {
390   AliGenerator::SetPtRange(ptmin, ptmax);
391 }
392