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