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