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