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