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