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