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