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