fRefVolumeId for reference volume identification added.
[u/mrichter/AliRoot.git] / EVGEN / AliGenMC.cxx
CommitLineData
e36044d6 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$
3b945a60 18Revision 1.14 2003/01/14 10:50:19 alibrary
19Cleanup of STEER coding conventions
20
116cbefd 21Revision 1.13 2002/10/14 14:55:35 hristov
22Merging the VirtualMC branch to the main development branch (HEAD)
23
b9d0a01d 24Revision 1.5.4.2 2002/07/24 08:56:28 alibrary
25Updating EVGEN on TVirtulaMC
26
27Revision 1.12 2002/07/19 11:42:33 morsch
28Use CalcMass()
29
cf01ee36 30Revision 1.11 2002/06/06 15:26:24 morsch
31Correct child-selection for kPhiKK
32
027677ae 33Revision 1.10 2002/06/05 14:05:46 morsch
34Decayer option kPhiKK for forced phi->K+K- decay added.
35
00cb2f27 36Revision 1.9 2002/05/30 14:58:29 morsch
37Add pointer to AliGeometry to handle geometrical acceptance. (G. MArtinez)
38
65de9f86 39Revision 1.8 2002/04/26 10:42:35 morsch
40Case kNoDecayHeavy added. (N. Carrer)
41
d5f86442 42Revision 1.7 2002/04/17 10:32:32 morsch
43Coding Rule violations corrected.
44
62d0ae06 45Revision 1.6 2002/03/26 14:19:36 morsch
46Saver calculation of rapdity.
47
fef61e14 48Revision 1.5 2002/03/12 17:02:20 morsch
49Change in calculation of rapidity, include case in which numerically e == pz.
50
95450c6d 51Revision 1.4 2001/11/27 13:13:07 morsch
52Maximum lifetime for long-lived particles to be put on the stack is parameter.
53It can be set via SetMaximumLifetime(..).
54
47fc6bd5 55Revision 1.3 2001/10/16 08:48:56 morsch
56Common vertex related code moved to base class AliGenerator.
57
b80a8e41 58Revision 1.2 2001/10/15 08:15:51 morsch
59Event vertex and vertex truncation setting moved into AliMC.
60
dba8ac91 61Revision 1.1 2001/07/13 10:56:00 morsch
62AliGenMC base class for AliGenParam and AliGenPythia commonalities.
63
e36044d6 64*/
65
62d0ae06 66// Base class for generators using external MC generators.
67// For example AliGenPythia using Pythia.
68// Provides basic functionality: setting of kinematic cuts on
69// decay products and particle selection.
70// andreas.morsch@cern.ch
71
65de9f86 72#include <TMath.h>
116cbefd 73#include <TPDGCode.h>
e36044d6 74#include <TParticle.h>
75
116cbefd 76#include "AliGenMC.h"
77
e36044d6 78 ClassImp(AliGenMC)
79
80AliGenMC::AliGenMC()
81 :AliGenerator()
82{
83// Default Constructor
84 SetCutOnChild();
85 SetChildMomentumRange();
86 SetChildPtRange();
87 SetChildPhiRange();
88 SetChildThetaRange();
89 SetChildYRange();
47fc6bd5 90 SetMaximumLifetime();
65de9f86 91 SetGeometryAcceptance();
92 SetPdgCodeParticleforAcceptanceCut();
93 SetNumberOfAcceptedParticles();
3b945a60 94 SetTarget();
95 SetProjectile();
e36044d6 96}
97
98AliGenMC::AliGenMC(Int_t npart)
99 :AliGenerator(npart)
100{
101// Constructor
102 SetCutOnChild();
103 SetChildMomentumRange();
104 SetChildPtRange();
105 SetChildPhiRange();
106 SetChildThetaRange();
107 SetChildYRange();
108//
109 fParentSelect.Set(8);
110 fChildSelect.Set(8);
111 for (Int_t i=0; i<8; i++) fParentSelect[i]=fChildSelect[i]=0;
47fc6bd5 112 SetMaximumLifetime();
65de9f86 113 SetGeometryAcceptance();
114 SetPdgCodeParticleforAcceptanceCut();
115 SetNumberOfAcceptedParticles();
3b945a60 116 SetTarget();
117 SetProjectile();
e36044d6 118}
119
120AliGenMC::AliGenMC(const AliGenMC & mc)
121{
122// copy constructor
123}
124
125AliGenMC::~AliGenMC()
126{
127// Destructor
128}
129
130void AliGenMC::Init()
131{
132//
133// Initialization
134 switch (fForceDecay) {
135 case kSemiElectronic:
136 case kDiElectron:
137 case kBJpsiDiElectron:
138 case kBPsiPrimeDiElectron:
47fc6bd5 139 fChildSelect[0] = kElectron;
e36044d6 140 break;
141 case kSemiMuonic:
142 case kDiMuon:
143 case kBJpsiDiMuon:
144 case kBPsiPrimeDiMuon:
145 case kPiToMu:
146 case kKaToMu:
147 fChildSelect[0]=kMuonMinus;
148 break;
149 case kHadronicD:
150 fChildSelect[0]=kPiPlus;
151 fChildSelect[1]=kKPlus;
152 break;
00cb2f27 153 case kPhiKK:
027677ae 154 fChildSelect[0]=kKPlus;
62d0ae06 155 case kOmega:
e36044d6 156 case kAll:
157 case kNoDecay:
d5f86442 158 case kNoDecayHeavy:
e36044d6 159 break;
160 }
161}
162
163
62d0ae06 164Bool_t AliGenMC::ParentSelected(Int_t ip) const
e36044d6 165{
166// True if particle is in list of parent particles to be selected
167 for (Int_t i=0; i<8; i++)
168 {
62d0ae06 169 if (fParentSelect.At(i) == ip) return kTRUE;
e36044d6 170 }
171 return kFALSE;
172}
173
62d0ae06 174Bool_t AliGenMC::ChildSelected(Int_t ip) const
e36044d6 175{
176// True if particle is in list of decay products to be selected
177 for (Int_t i=0; i<5; i++)
178 {
62d0ae06 179 if (fChildSelect.At(i) == ip) return kTRUE;
e36044d6 180 }
181 return kFALSE;
182}
183
62d0ae06 184Bool_t AliGenMC::KinematicSelection(TParticle *particle, Int_t flag) const
e36044d6 185{
186// Perform kinematic selection
187 Float_t px = particle->Px();
188 Float_t py = particle->Py();
189 Float_t pz = particle->Pz();
190 Float_t e = particle->Energy();
191 Float_t pt = particle->Pt();
192 Float_t p = particle->P();
193 Float_t theta = particle->Theta();
cf01ee36 194 Float_t mass = particle->GetCalcMass();
95450c6d 195 Float_t mt2 = pt * pt + mass * mass;
196
e36044d6 197 Float_t phi = Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
95450c6d 198 Double_t y, y0;
199
fef61e14 200 if (TMath::Abs(pz) < e) {
95450c6d 201 y = 0.5*TMath::Log((e+pz)/(e-pz));
202 } else {
203 y = 1.e10;
204 }
e36044d6 205
95450c6d 206 if (mt2) {
207 y0 = 0.5*TMath::Log((e+TMath::Abs(pz))*(e+TMath::Abs(pz))/mt2);
e36044d6 208 } else {
95450c6d 209 if (TMath::Abs(y) < 1.e10) {
210 y0 = y;
211 } else {
212 y0 = 1.e10;
213 }
e36044d6 214 }
95450c6d 215
216 y = (pz < 0) ? -y0 : y0;
e36044d6 217
218 if (flag == 0) {
219//
220// Primary particle cuts
221//
222// transverse momentum cut
223 if (pt > fPtMax || pt < fPtMin) {
224// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
225 return kFALSE;
226 }
227//
228// momentum cut
229 if (p > fPMax || p < fPMin) {
230// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
231 return kFALSE;
232 }
233//
234// theta cut
235 if (theta > fThetaMax || theta < fThetaMin) {
236// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
237 return kFALSE;
238 }
239//
240// rapidity cut
241 if (y > fYMax || y < fYMin) {
242// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
243 return kFALSE;
244 }
245//
246// phi cut
247 if (phi > fPhiMax || phi < fPhiMin) {
248// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
249 return kFALSE;
250 }
251 } else {
252//
253// Decay product cuts
254//
255// transverse momentum cut
256 if (pt > fChildPtMax || pt < fChildPtMin) {
257// printf("\n failed pt cut %f %f %f \n",pt,fChildPtMin,fChildPtMax);
258 return kFALSE;
259 }
260//
261// momentum cut
262 if (p > fChildPMax || p < fChildPMin) {
263// printf("\n failed p cut %f %f %f \n",p,fChildPMin,fChildPMax);
264 return kFALSE;
265 }
266//
267// theta cut
268 if (theta > fChildThetaMax || theta < fChildThetaMin) {
269// printf("\n failed theta cut %f %f %f \n",theta,fChildThetaMin,fChildThetaMax);
270 return kFALSE;
271 }
272//
273// rapidity cut
274 if (y > fChildYMax || y < fChildYMin) {
275// printf("\n failed y cut %f %f %f \n",y,fChildYMin,fChildYMax);
276 return kFALSE;
277 }
278//
279// phi cut
280 if (phi > fChildPhiMax || phi < fChildPhiMin) {
281// printf("\n failed phi cut %f %f %f \n",phi,fChildPhiMin,fChildPhiMax);
282 return kFALSE;
283 }
284 }
285
e36044d6 286 return kTRUE;
287}
288
65de9f86 289Bool_t AliGenMC::CheckAcceptanceGeometry(Int_t np, TClonesArray* particles)
290{
291 Bool_t Check ; // All fPdgCodeParticleforAcceptanceCut particles are in in the fGeometryAcceptance acceptance
292 Int_t NumberOfPdgCodeParticleforAcceptanceCut=0;
293 Int_t NumberOfAcceptedPdgCodeParticleforAcceptanceCut=0;
294 TParticle * particle;
295 Int_t i;
296 for (i=0; i<np; i++) {
297 particle = (TParticle *) particles->At(i);
298 if( TMath::Abs( particle->GetPdgCode() ) == TMath::Abs( fPdgCodeParticleforAcceptanceCut ) ) {
299 NumberOfPdgCodeParticleforAcceptanceCut++;
300 if (fGeometryAcceptance->Impact(particle)) NumberOfAcceptedPdgCodeParticleforAcceptanceCut++;
301 }
302 }
303 if ( NumberOfAcceptedPdgCodeParticleforAcceptanceCut > (fNumberOfAcceptedParticles-1) )
304 Check = kTRUE;
305 else
306 Check = kFALSE;
307
308 return Check;
309}
310
62d0ae06 311Int_t AliGenMC::CheckPDGCode(Int_t pdgcode) const
e36044d6 312{
313//
314// If the particle is in a diffractive state, then take action accordingly
315 switch (pdgcode) {
316 case 91:
317 return 92;
318 case 110:
319 //rho_diff0 -- difficult to translate, return rho0
320 return 113;
321 case 210:
322 //pi_diffr+ -- change to pi+
323 return 211;
324 case 220:
325 //omega_di0 -- change to omega0
326 return 223;
327 case 330:
328 //phi_diff0 -- return phi0
329 return 333;
330 case 440:
331 //J/psi_di0 -- return J/psi
332 return 443;
333 case 2110:
334 //n_diffr -- return neutron
335 return 2112;
336 case 2210:
337 //p_diffr+ -- return proton
338 return 2212;
339 }
340 //non diffractive state -- return code unchanged
341 return pdgcode;
342}
3b945a60 343
344void AliGenMC::Boost(Float_t dy)
345{
346//
347// Boost cms into LHC lab frame
348//
349
350 Double_t beta = TMath::TanH(dy);
351 Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
352 Double_t gb = gamma * beta;
353
354 printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma);
355
356 Int_t i;
357 Int_t np = fParticles->GetEntriesFast();
358 for (i = 0; i < np; i++)
359 {
360 TParticle* iparticle = (TParticle*) fParticles->At(i);
361
362 Double_t e = iparticle->Energy();
363 Double_t px = iparticle->Px();
364 Double_t py = iparticle->Py();
365 Double_t pz = iparticle->Pz();
366
367 Double_t eb = gamma * e - gb * pz;
368 Double_t pzb = -gb * e + gamma * pz;
369
370 iparticle->SetMomentum(px, py, pzb, eb);
371 }
372}
373
374
e36044d6 375
376AliGenMC& AliGenMC::operator=(const AliGenMC& rhs)
377{
378// Assignment operator
379 return *this;
380}
381