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