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