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