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