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