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