ty is Double_t
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.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.34  2001/11/27 13:13:07  morsch
19 Maximum lifetime for long-lived particles to be put on the stack is parameter.
20 It can be set via SetMaximumLifetime(..).
21
22 Revision 1.33  2001/10/21 18:35:56  hristov
23 Several pointers were set to zero in the default constructors to avoid memory management problems
24
25 Revision 1.32  2001/07/27 17:09:36  morsch
26 Use local SetTrack, KeepTrack and SetHighWaterMark methods
27 to delegate either to local stack or to stack owned by AliRun.
28 (Piotr Skowronski, A.M.)
29
30 Revision 1.31  2001/07/13 10:58:54  morsch
31 - Some coded moved to AliGenMC
32 - Improved handling of secondary vertices.
33
34 Revision 1.30  2001/06/15 07:55:04  morsch
35 Put only first generation decay products on the stack.
36
37 Revision 1.29  2001/03/27 10:58:41  morsch
38 Initialize decayer before generation. Important if run inside cocktail.
39
40 Revision 1.28  2001/03/09 13:01:41  morsch
41 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
42 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
43
44 Revision 1.27  2001/02/02 15:21:10  morsch
45 Set high water mark after last particle.
46 Use Vertex() method for Vertex.
47
48 Revision 1.26  2000/12/21 16:24:06  morsch
49 Coding convention clean-up
50
51 Revision 1.25  2000/11/30 07:12:50  alibrary
52 Introducing new Rndm and QA classes
53
54 Revision 1.24  2000/10/18 19:11:27  hristov
55 Division by zero fixed
56
57 Revision 1.23  2000/10/02 21:28:06  fca
58 Removal of useless dependecies via forward declarations
59
60 Revision 1.22  2000/09/12 14:14:55  morsch
61 Call fDecayer->ForceDecay() at the beginning of Generate().
62
63 Revision 1.21  2000/09/08 15:39:01  morsch
64 Handle the case fForceDecay=all during the generation, i.e. select all secondaries.
65
66 Revision 1.20  2000/09/06 14:35:44  morsch
67 Use AliDecayerPythia for particle decays.
68
69 Revision 1.19  2000/07/11 18:24:56  fca
70 Coding convention corrections + few minor bug fixes
71
72 Revision 1.18  2000/06/29 21:08:27  morsch
73 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
74 This allows to pass a pointer to a library directly to AliGenParam and avoids the
75 use of function pointers in Config.C.
76
77 Revision 1.17  2000/06/09 20:33:30  morsch
78 All coding rule violations except RS3 corrected
79
80 Revision 1.16  2000/05/02 07:51:31  morsch
81 - Control precision of pT sampling TF1::SetNpx(..)
82 - Correct initialisation of child-cuts in all constructors.
83 - Most coding rule violations corrected.
84
85 Revision 1.15  2000/04/03 15:42:12  morsch
86 Cuts on primary particles are separated from those on the decay products. Methods
87 SetChildMomentumRange, SetChildPtRange, SetChildPhiRange, SetChildThetaRange added.
88
89 Revision 1.14  1999/11/09 07:38:48  fca
90 Changes for compatibility with version 2.23 of ROOT
91
92 Revision 1.13  1999/11/04 11:30:31  fca
93 Correct the logics for SetForceDecay
94
95 Revision 1.12  1999/11/03 17:43:20  fca
96 New version from G.Martinez & A.Morsch
97
98 Revision 1.11  1999/09/29 09:24:14  fca
99 Introduction of the Copyright and cvs Log
100
101 */
102
103 #include "AliGenParam.h"
104 #include "AliDecayerPythia.h"
105 #include "AliGenMUONlib.h"
106 #include "AliRun.h"
107 #include <TParticle.h>
108 #include <TParticlePDG.h>
109 #include <TDatabasePDG.h>
110 #include <TLorentzVector.h>
111
112 #include <TF1.h>
113
114 ClassImp(AliGenParam)
115
116 //------------------------------------------------------------
117
118   //Begin_Html
119   /*
120     <img src="picts/AliGenParam.gif">
121   */
122   //End_Html
123
124 //____________________________________________________________
125 //____________________________________________________________
126 AliGenParam::AliGenParam()
127 {
128 // Deafault constructor
129     fPtPara = 0;
130     fYPara  = 0;
131     fParam  = 0;
132     fAnalog = kAnalog;
133     SetDeltaPt();
134 //
135 // Set random number generator   
136     sRandom = fRandom;
137     fDecayer = 0;
138 }
139
140 AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname):AliGenMC(npart)
141 {
142 // Constructor using number of particles parameterisation id and library
143     
144     fPtParaFunc = Library->GetPt(param, tname);
145     fYParaFunc  = Library->GetY (param, tname);
146     fIpParaFunc = Library->GetIp(param, tname);
147     
148     fPtPara = 0;
149     fYPara  = 0;
150     fParam  = param;
151     fAnalog = kAnalog;
152     SetForceDecay();
153     SetDeltaPt(); 
154 //
155 // Set random number generator   
156     sRandom = fRandom;
157 }
158
159 //____________________________________________________________
160
161 AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenMC(npart)
162 {
163 // Constructor using parameterisation id and number of particles
164 //      
165     AliGenLib* pLibrary = new AliGenMUONlib();
166  
167     fPtParaFunc = pLibrary->GetPt(param, tname);
168     fYParaFunc  = pLibrary->GetY (param, tname);
169     fIpParaFunc = pLibrary->GetIp(param, tname);
170     
171     fPtPara = 0;
172     fYPara  = 0;
173     fParam  = param;
174     fAnalog = kAnalog;
175     fChildSelect.Set(5);
176     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
177     SetForceDecay();
178     SetCutOnChild();
179     SetChildMomentumRange();
180     SetChildPtRange();
181     SetChildPhiRange();
182     SetChildThetaRange(); 
183     SetDeltaPt(); 
184 }
185
186 AliGenParam::AliGenParam(Int_t npart, Int_t param,
187                          Double_t (*PtPara) (Double_t*, Double_t*),
188                          Double_t (*YPara ) (Double_t* ,Double_t*),
189                          Int_t    (*IpPara) (TRandom *))                 
190     :AliGenMC(npart)
191 {
192 // Constructor
193 // Gines Martinez 1/10/99 
194     fPtParaFunc = PtPara; 
195     fYParaFunc  = YPara;  
196     fIpParaFunc = IpPara;
197 //  
198     fPtPara = 0;
199     fYPara  = 0;
200     fParam  = param;
201     fAnalog = kAnalog;
202     fChildSelect.Set(5);
203     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
204     SetForceDecay();
205     SetCutOnChild();
206     SetChildMomentumRange();
207     SetChildPtRange();
208     SetChildPhiRange();
209     SetChildThetaRange();  
210     SetDeltaPt();
211 }
212
213
214 AliGenParam::AliGenParam(const AliGenParam & Paramd)
215 {
216 // copy constructor
217 }
218
219 //____________________________________________________________
220 AliGenParam::~AliGenParam()
221 {
222 // Destructor
223     delete  fPtPara;
224     delete  fYPara;
225 }
226
227 //____________________________________________________________
228 void AliGenParam::Init()
229 {
230 // Initialisation
231
232     fDecayer = new AliDecayerPythia();
233   //Begin_Html
234   /*
235     <img src="picts/AliGenParam.gif">
236   */
237   //End_Html
238  
239     fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,fPtMin,fPtMax,0);
240 //  Set representation precision to 10 MeV
241     Int_t npx= Int_t((fPtMax-fPtMin)/fDeltaPt);
242     
243     fPtPara->SetNpx(npx);
244     
245     fYPara  = new TF1("Y -Parametrization",fYParaFunc,fYMin,fYMax,0);
246     TF1* ptPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
247     TF1* yPara  = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);
248
249 //
250 // dN/dy| y=0
251     Double_t y1=0;
252     Double_t y2=0;
253     
254     fdNdy0=fYParaFunc(&y1,&y2);
255 //
256 // Integral over generation region
257     Float_t intYS  = yPara ->Integral(fYMin, fYMax);
258     Float_t intPt0 = ptPara->Integral(0,15);
259     Float_t intPtS = ptPara->Integral(fPtMin,fPtMax);
260     Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
261     if (fAnalog == kAnalog) {
262         fYWgt  = intYS/fdNdy0;
263         fPtWgt = intPtS/intPt0;
264         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
265     } else {
266         fYWgt = intYS/fdNdy0;
267         fPtWgt = (fPtMax-fPtMin)/intPt0;
268         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
269     }
270 //
271 // particle decay related initialization
272     fDecayer->SetForceDecay(fForceDecay);
273     fDecayer->Init();
274
275 //
276     AliGenMC::Init();
277 }
278
279 //____________________________________________________________
280 void AliGenParam::Generate()
281 {
282 //
283 // Generate 'npart' of light and heavy mesons (J/Psi, upsilon or phi, Pion,
284 // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and 
285 // antineutrons in the the desired theta, phi and momentum windows; 
286 // Gaussian smearing on the vertex is done if selected. 
287 // The decay of heavy mesons is done using lujet, 
288 //    and the childern particle are tracked by GEANT
289 // However, light mesons are directly tracked by GEANT 
290 // setting fForceDecay = nodecay (SetForceDecay(nodecay)) 
291 //
292 //
293 //  Reinitialize decayer
294     fDecayer->Init();
295 //
296   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
297   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
298   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
299   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
300   Float_t p[3], pc[3], 
301           och[3];             // Momentum, polarisation and origin of the children particles from lujet
302   Double_t ty, xmt;
303   Int_t nt, i, j;
304   Float_t  wgtp, wgtch;
305   Double_t dummy;
306   static TClonesArray *particles;
307   //
308   if(!particles) particles = new TClonesArray("TParticle",1000);
309
310   static TDatabasePDG *pDataBase = new TDatabasePDG();
311   if(!pDataBase) pDataBase = new TDatabasePDG();
312   //
313   Float_t random[6];
314  
315 // Calculating vertex position per event
316   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
317   if(fVertexSmear==kPerEvent) {
318       Vertex();
319       for (j=0;j<3;j++) origin0[j]=fVertex[j];
320   }
321   
322   Int_t ipa=0;
323   
324 // Generating fNpart particles
325   while (ipa<fNpart) {
326       while(1) {
327 //
328 // particle type
329           Int_t iPart = fIpParaFunc(fRandom);
330           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;          
331           TParticlePDG *particle = pDataBase->GetParticle(iPart);
332           Float_t am = particle->Mass();
333
334           Rndm(random,2);
335 //
336 // phi
337           phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
338 //
339 // y
340           ty = TMath::TanH(fYPara->GetRandom());
341 //
342 // pT
343           if (fAnalog == kAnalog) {
344               pt=fPtPara->GetRandom();
345               wgtp=fParentWeight;
346               wgtch=fChildWeight;
347           } else {
348               pt=fPtMin+random[1]*(fPtMax-fPtMin);
349               Double_t ptd=pt;
350               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
351               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
352           }
353           xmt=sqrt(pt*pt+am*am);
354           if (TMath::Abs(ty)==1.) {
355               ty=0.;
356               Fatal("AliGenParam", 
357                     "Division by 0: Please check you rapidity range !");
358           }
359           
360           pl=xmt*ty/sqrt(1.-ty*ty);
361           theta=TMath::ATan2(pt,pl);
362 // Cut on theta
363           if(theta<fThetaMin || theta>fThetaMax) continue;
364           ptot=TMath::Sqrt(pt*pt+pl*pl);
365 // Cut on momentum
366           if(ptot<fPMin || ptot>fPMax) continue;
367 //
368           p[0]=pt*TMath::Cos(phi);
369           p[1]=pt*TMath::Sin(phi);
370           p[2]=pl;
371           if(fVertexSmear==kPerTrack) {
372               Rndm(random,6);
373               for (j=0;j<3;j++) {
374                   origin0[j]=
375                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
376                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
377               }
378           }
379           
380 // Looking at fForceDecay : 
381 // if fForceDecay != none Primary particle decays using 
382 // AliPythia and children are tracked by GEANT
383 //
384 // if fForceDecay == none Primary particle is tracked by GEANT 
385 // (In the latest, make sure that GEANT actually does all the decays you want)    
386 //
387
388           if (fForceDecay != kNoDecay) {
389 // Using lujet to decay particle
390               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
391               TLorentzVector pmom(p[0], p[1], p[2], energy);
392               fDecayer->Decay(iPart,&pmom);
393 //
394 // select decay particles
395               Int_t np=fDecayer->ImportParticles(particles);
396               Int_t ncsel=0;
397               Int_t* pFlag      = new Int_t[np];
398               Int_t* pParent    = new Int_t[np];
399               Int_t* pSelected  = new Int_t[np];
400               Int_t* trackIt    = new Int_t[np];
401
402               for (i=0; i<np; i++) {
403                   pFlag[i]     =  0;
404                   pSelected[i] =  0;
405                   pParent[i]   = -1;
406               }
407               
408               if (np >1) {
409                   TParticle* iparticle =  (TParticle *) particles->At(0);
410                   Int_t ipF, ipL;
411                   for (i = 1; i<np ; i++) {
412                       trackIt[i] = 1;
413                       iparticle = (TParticle *) particles->At(i);
414                       Int_t kf = iparticle->GetPdgCode();
415                       Int_t ks = iparticle->GetStatusCode();
416 // flagged particle
417
418                       if (pFlag[i] == 1) {
419                           ipF = iparticle->GetFirstDaughter();
420                           ipL = iparticle->GetLastDaughter();   
421                           if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
422                           continue;
423                       }
424
425 // flag decay products of particles with long life-time (c tau > .3 mum)                      
426                       
427                       if (ks != 1) { 
428 //                        TParticlePDG *particle = pDataBase->GetParticle(kf);
429                           
430                           Double_t lifeTime = fDecayer->GetLifetime(kf);
431 //                        Double_t mass     = particle->Mass();
432 //                        Double_t width    = particle->Width();
433                           if (lifeTime > (Double_t) fMaxLifeTime) {
434                               ipF = iparticle->GetFirstDaughter();
435                               ipL = iparticle->GetLastDaughter();       
436                               if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
437                           } else{
438                               trackIt[i]     = 0;
439                               pSelected[i]   = 1;
440                           }
441                       } // ks==1 ?
442 //
443 // children
444                       
445                       if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll && trackIt[i])
446                       {
447                           if (fCutOnChild) {
448                               pc[0]=iparticle->Px();
449                               pc[1]=iparticle->Py();
450                               pc[2]=iparticle->Pz();
451                               Bool_t  childok = KinematicSelection(iparticle, 1);
452                               if(childok) {
453                                   pSelected[i]  = 1;
454                                   ncsel++;
455                               } else {
456                                   ncsel=-1;
457                                   break;
458                               } // child kine cuts
459                           } else {
460                               pSelected[i]  = 1;
461                               ncsel++;
462                           } // if child selection
463                       } // select muon
464                   } // decay particle loop
465               } // if decay products
466               
467               Int_t iparent;
468               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
469                   ipa++;
470 //
471 // Parent
472                   SetTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
473                   pParent[0] = nt;
474                   KeepTrack(nt); 
475 //
476 // Decay Products
477 //                
478                   for (i = 1; i < np; i++) {
479                       if (pSelected[i]) {
480                           TParticle* iparticle = (TParticle *) particles->At(i);
481                           Int_t kf  = iparticle->GetPdgCode();
482                           Int_t ipa = iparticle->GetFirstMother()-1;
483                           
484                           och[0] = origin0[0]+iparticle->Vx()/10;
485                           och[1] = origin0[1]+iparticle->Vy()/10;
486                           och[2] = origin0[2]+iparticle->Vz()/10;
487                           pc[0]  = iparticle->Px();
488                           pc[1]  = iparticle->Py();
489                           pc[2]  = iparticle->Pz();
490                           
491                           if (ipa > -1) {
492                               iparent = pParent[ipa];
493                           } else {
494                               iparent = -1;
495                           }
496                           SetTrack(fTrackIt*trackIt[i], iparent, kf,
497                                            pc, och, polar,
498                                            0, kPDecay, nt, wgtch);
499                           pParent[i] = nt;
500                           KeepTrack(nt); 
501                       }
502                   }
503               }  // Decays by Lujet
504
505               if (pFlag)      delete[] pFlag;
506               if (pParent)    delete[] pParent;
507               if (pSelected)  delete[] pSelected;          
508               if (trackIt)    delete[] trackIt;     
509           } // kinematic selection
510           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
511           {
512             gAlice->
513                 SetTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
514             ipa++; 
515           }
516           break;
517     } // while
518   } // event loop
519   SetHighWaterMark(nt);
520 }
521
522 AliGenParam& AliGenParam::operator=(const  AliGenParam& rhs)
523 {
524 // Assignment operator
525     return *this;
526 }
527
528
529