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