Coding Rule violations corrected.
[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
88cb7938 16/* $Id$ */
e36044d6 17
62d0ae06 18// Base class for generators using external MC generators.
19// For example AliGenPythia using Pythia.
20// Provides basic functionality: setting of kinematic cuts on
21// decay products and particle selection.
22// andreas.morsch@cern.ch
23
65de9f86 24#include <TMath.h>
116cbefd 25#include <TPDGCode.h>
e36044d6 26#include <TParticle.h>
27
116cbefd 28#include "AliGenMC.h"
29
198bb1c7 30ClassImp(AliGenMC)
e36044d6 31
32AliGenMC::AliGenMC()
198bb1c7 33 :AliGenerator()
e36044d6 34{
35// Default Constructor
36 SetCutOnChild();
37 SetChildMomentumRange();
38 SetChildPtRange();
39 SetChildPhiRange();
40 SetChildThetaRange();
41 SetChildYRange();
47fc6bd5 42 SetMaximumLifetime();
65de9f86 43 SetGeometryAcceptance();
44 SetPdgCodeParticleforAcceptanceCut();
45 SetNumberOfAcceptedParticles();
3b945a60 46 SetTarget();
47 SetProjectile();
e36044d6 48}
49
50AliGenMC::AliGenMC(Int_t npart)
198bb1c7 51 :AliGenerator(npart)
e36044d6 52{
53// Constructor
54 SetCutOnChild();
55 SetChildMomentumRange();
56 SetChildPtRange();
57 SetChildPhiRange();
58 SetChildThetaRange();
59 SetChildYRange();
60//
61 fParentSelect.Set(8);
62 fChildSelect.Set(8);
63 for (Int_t i=0; i<8; i++) fParentSelect[i]=fChildSelect[i]=0;
47fc6bd5 64 SetMaximumLifetime();
65de9f86 65 SetGeometryAcceptance();
66 SetPdgCodeParticleforAcceptanceCut();
67 SetNumberOfAcceptedParticles();
3b945a60 68 SetTarget();
69 SetProjectile();
e36044d6 70}
71
198bb1c7 72AliGenMC::AliGenMC(const AliGenMC & mc):
73 AliGenerator(mc)
e36044d6 74{
198bb1c7 75// Copy constructor
76 mc.Copy(*this);
e36044d6 77}
78
79AliGenMC::~AliGenMC()
80{
81// Destructor
82}
83
84void AliGenMC::Init()
85{
86//
87// Initialization
88 switch (fForceDecay) {
89 case kSemiElectronic:
90 case kDiElectron:
91 case kBJpsiDiElectron:
92 case kBPsiPrimeDiElectron:
47fc6bd5 93 fChildSelect[0] = kElectron;
e36044d6 94 break;
2dcb5874 95 case kHardMuons:
e36044d6 96 case kSemiMuonic:
97 case kDiMuon:
98 case kBJpsiDiMuon:
99 case kBPsiPrimeDiMuon:
100 case kPiToMu:
101 case kKaToMu:
102 fChildSelect[0]=kMuonMinus;
103 break;
104 case kHadronicD:
105 fChildSelect[0]=kPiPlus;
106 fChildSelect[1]=kKPlus;
107 break;
00cb2f27 108 case kPhiKK:
027677ae 109 fChildSelect[0]=kKPlus;
62d0ae06 110 case kOmega:
e36044d6 111 case kAll:
112 case kNoDecay:
d5f86442 113 case kNoDecayHeavy:
e36044d6 114 break;
115 }
36d770d9 116
117 if (fZTarget != 0 && fAProjectile != 0)
118 {
119 fDyBoost = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) /
120 (Double_t(fZTarget) * Double_t(fAProjectile)));
121 }
e36044d6 122}
123
124
62d0ae06 125Bool_t AliGenMC::ParentSelected(Int_t ip) const
e36044d6 126{
127// True if particle is in list of parent particles to be selected
128 for (Int_t i=0; i<8; i++)
129 {
62d0ae06 130 if (fParentSelect.At(i) == ip) return kTRUE;
e36044d6 131 }
132 return kFALSE;
133}
134
62d0ae06 135Bool_t AliGenMC::ChildSelected(Int_t ip) const
e36044d6 136{
137// True if particle is in list of decay products to be selected
138 for (Int_t i=0; i<5; i++)
139 {
62d0ae06 140 if (fChildSelect.At(i) == ip) return kTRUE;
e36044d6 141 }
142 return kFALSE;
143}
144
62d0ae06 145Bool_t AliGenMC::KinematicSelection(TParticle *particle, Int_t flag) const
e36044d6 146{
147// Perform kinematic selection
e36044d6 148 Float_t pz = particle->Pz();
149 Float_t e = particle->Energy();
150 Float_t pt = particle->Pt();
151 Float_t p = particle->P();
152 Float_t theta = particle->Theta();
cf01ee36 153 Float_t mass = particle->GetCalcMass();
95450c6d 154 Float_t mt2 = pt * pt + mass * mass;
0b6d2cb3 155 Float_t phi = particle->Phi();
95450c6d 156
95450c6d 157 Double_t y, y0;
158
fef61e14 159 if (TMath::Abs(pz) < e) {
95450c6d 160 y = 0.5*TMath::Log((e+pz)/(e-pz));
161 } else {
162 y = 1.e10;
163 }
e36044d6 164
95450c6d 165 if (mt2) {
166 y0 = 0.5*TMath::Log((e+TMath::Abs(pz))*(e+TMath::Abs(pz))/mt2);
e36044d6 167 } else {
95450c6d 168 if (TMath::Abs(y) < 1.e10) {
169 y0 = y;
170 } else {
171 y0 = 1.e10;
172 }
e36044d6 173 }
95450c6d 174
175 y = (pz < 0) ? -y0 : y0;
e36044d6 176
177 if (flag == 0) {
178//
179// Primary particle cuts
180//
181// transverse momentum cut
182 if (pt > fPtMax || pt < fPtMin) {
183// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
184 return kFALSE;
185 }
186//
187// momentum cut
188 if (p > fPMax || p < fPMin) {
189// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
190 return kFALSE;
191 }
192//
193// theta cut
194 if (theta > fThetaMax || theta < fThetaMin) {
195// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
196 return kFALSE;
197 }
198//
199// rapidity cut
200 if (y > fYMax || y < fYMin) {
201// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
202 return kFALSE;
203 }
204//
205// phi cut
206 if (phi > fPhiMax || phi < fPhiMin) {
207// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
208 return kFALSE;
209 }
210 } else {
211//
212// Decay product cuts
213//
214// transverse momentum cut
215 if (pt > fChildPtMax || pt < fChildPtMin) {
216// printf("\n failed pt cut %f %f %f \n",pt,fChildPtMin,fChildPtMax);
217 return kFALSE;
218 }
219//
220// momentum cut
221 if (p > fChildPMax || p < fChildPMin) {
222// printf("\n failed p cut %f %f %f \n",p,fChildPMin,fChildPMax);
223 return kFALSE;
224 }
225//
226// theta cut
227 if (theta > fChildThetaMax || theta < fChildThetaMin) {
228// printf("\n failed theta cut %f %f %f \n",theta,fChildThetaMin,fChildThetaMax);
229 return kFALSE;
230 }
231//
232// rapidity cut
233 if (y > fChildYMax || y < fChildYMin) {
234// printf("\n failed y cut %f %f %f \n",y,fChildYMin,fChildYMax);
235 return kFALSE;
236 }
237//
238// phi cut
239 if (phi > fChildPhiMax || phi < fChildPhiMin) {
240// printf("\n failed phi cut %f %f %f \n",phi,fChildPhiMin,fChildPhiMax);
241 return kFALSE;
242 }
243 }
244
e36044d6 245 return kTRUE;
246}
247
65de9f86 248Bool_t AliGenMC::CheckAcceptanceGeometry(Int_t np, TClonesArray* particles)
249{
250 Bool_t Check ; // All fPdgCodeParticleforAcceptanceCut particles are in in the fGeometryAcceptance acceptance
251 Int_t NumberOfPdgCodeParticleforAcceptanceCut=0;
252 Int_t NumberOfAcceptedPdgCodeParticleforAcceptanceCut=0;
253 TParticle * particle;
254 Int_t i;
255 for (i=0; i<np; i++) {
256 particle = (TParticle *) particles->At(i);
257 if( TMath::Abs( particle->GetPdgCode() ) == TMath::Abs( fPdgCodeParticleforAcceptanceCut ) ) {
258 NumberOfPdgCodeParticleforAcceptanceCut++;
259 if (fGeometryAcceptance->Impact(particle)) NumberOfAcceptedPdgCodeParticleforAcceptanceCut++;
260 }
261 }
262 if ( NumberOfAcceptedPdgCodeParticleforAcceptanceCut > (fNumberOfAcceptedParticles-1) )
263 Check = kTRUE;
264 else
265 Check = kFALSE;
266
267 return Check;
268}
269
62d0ae06 270Int_t AliGenMC::CheckPDGCode(Int_t pdgcode) const
e36044d6 271{
272//
273// If the particle is in a diffractive state, then take action accordingly
274 switch (pdgcode) {
275 case 91:
276 return 92;
277 case 110:
278 //rho_diff0 -- difficult to translate, return rho0
279 return 113;
280 case 210:
281 //pi_diffr+ -- change to pi+
282 return 211;
283 case 220:
284 //omega_di0 -- change to omega0
285 return 223;
286 case 330:
287 //phi_diff0 -- return phi0
288 return 333;
289 case 440:
290 //J/psi_di0 -- return J/psi
291 return 443;
292 case 2110:
293 //n_diffr -- return neutron
294 return 2112;
295 case 2210:
296 //p_diffr+ -- return proton
297 return 2212;
298 }
299 //non diffractive state -- return code unchanged
300 return pdgcode;
301}
3b945a60 302
36d770d9 303void AliGenMC::Boost()
3b945a60 304{
305//
306// Boost cms into LHC lab frame
307//
308
36d770d9 309 Double_t beta = TMath::TanH(fDyBoost);
3b945a60 310 Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
311 Double_t gb = gamma * beta;
312
7cdba479 313 // printf("\n Boosting particles to lab frame %f %f %f", fDyBoost, beta, gamma);
3b945a60 314
315 Int_t i;
316 Int_t np = fParticles->GetEntriesFast();
317 for (i = 0; i < np; i++)
318 {
319 TParticle* iparticle = (TParticle*) fParticles->At(i);
320
321 Double_t e = iparticle->Energy();
322 Double_t px = iparticle->Px();
323 Double_t py = iparticle->Py();
324 Double_t pz = iparticle->Pz();
325
326 Double_t eb = gamma * e - gb * pz;
327 Double_t pzb = -gb * e + gamma * pz;
328
329 iparticle->SetMomentum(px, py, pzb, eb);
330 }
331}
332
333
e36044d6 334
335AliGenMC& AliGenMC::operator=(const AliGenMC& rhs)
336{
337// Assignment operator
198bb1c7 338 rhs.Copy(*this);
e36044d6 339 return *this;
340}
341
198bb1c7 342void AliGenMC::Copy(AliGenMC&) const
343{
344 //
345 // Copy
346 //
347 Fatal("Copy","Not implemented!\n");
348}
349
350