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