Line data Source code
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 : /* $Id$ */
17 :
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 :
24 : #include <TClonesArray.h>
25 : #include <TMath.h>
26 : #include <TPDGCode.h>
27 : #include <TParticle.h>
28 : #include <TLorentzVector.h>
29 : #include <TVector3.h>
30 :
31 : #include "AliGenMC.h"
32 : #include "AliGenEventHeader.h"
33 : #include "AliRun.h"
34 : #include "AliGeometry.h"
35 : #include "AliDecayer.h"
36 :
37 6 : ClassImp(AliGenMC)
38 :
39 : AliGenMC::AliGenMC()
40 0 : :AliGenerator(),
41 0 : fParticles("TParticle", 1000),
42 0 : fParentSelect(8),
43 0 : fChildSelect(8),
44 0 : fCutOnChild(0),
45 0 : fChildPtMin(0.),
46 0 : fChildPtMax(1.e10),
47 0 : fChildPMin(0.),
48 0 : fChildPMax(1.e10),
49 0 : fChildPhiMin(0.),
50 0 : fChildPhiMax(2. * TMath::Pi()),
51 0 : fChildThetaMin(0.),
52 0 : fChildThetaMax(TMath::Pi()),
53 0 : fChildYMin(-12.),
54 0 : fChildYMax(12.),
55 0 : fXingAngleX(0.),
56 0 : fXingAngleY(0.),
57 0 : fForceDecay(kAll),
58 0 : fMaxLifeTime(1.e-15),
59 0 : fDyBoost(0.),
60 0 : fGeometryAcceptance(0),
61 0 : fPdgCodeParticleforAcceptanceCut(0),
62 0 : fNumberOfAcceptedParticles(0),
63 0 : fNprimaries(0)
64 0 : {
65 : // Default Constructor
66 0 : fAProjectile = 1;
67 0 : fZProjectile = 1;
68 0 : fATarget = 1;
69 0 : fZTarget = 1;
70 0 : fProjectile = "P";
71 0 : fTarget = "P";
72 0 : }
73 :
74 : AliGenMC::AliGenMC(Int_t npart)
75 0 : :AliGenerator(npart),
76 0 : fParticles("TParticle", 1000),
77 0 : fParentSelect(8),
78 0 : fChildSelect(8),
79 0 : fCutOnChild(0),
80 0 : fChildPtMin(0.),
81 0 : fChildPtMax(1.e10),
82 0 : fChildPMin(0.),
83 0 : fChildPMax(1.e10),
84 0 : fChildPhiMin(0.),
85 0 : fChildPhiMax(2. * TMath::Pi()),
86 0 : fChildThetaMin(0.),
87 0 : fChildThetaMax(TMath::Pi()),
88 0 : fChildYMin(-12.),
89 0 : fChildYMax(12.),
90 0 : fXingAngleX(0.),
91 0 : fXingAngleY(0.),
92 0 : fForceDecay(kAll),
93 0 : fMaxLifeTime(1.e-15),
94 0 : fDyBoost(0.),
95 0 : fGeometryAcceptance(0),
96 0 : fPdgCodeParticleforAcceptanceCut(0),
97 0 : fNumberOfAcceptedParticles(0),
98 0 : fNprimaries(0)
99 0 : {
100 : // Constructor
101 : //
102 0 : fAProjectile = 1;
103 0 : fZProjectile = 1;
104 0 : fATarget = 1;
105 0 : fZTarget = 1;
106 0 : fProjectile = "P";
107 0 : fTarget = "P";
108 0 : for (Int_t i=0; i<8; i++) fParentSelect[i]=fChildSelect[i]=0;
109 0 : }
110 :
111 0 : AliGenMC::~AliGenMC()
112 0 : {
113 : // Destructor
114 0 : }
115 :
116 : void AliGenMC::Init()
117 : {
118 : //
119 : // Initialization
120 0 : switch (fForceDecay) {
121 : case kBSemiElectronic:
122 : case kSemiElectronic:
123 : case kDiElectron:
124 : case kBJpsiDiElectron:
125 : case kBPsiPrimeDiElectron:
126 : case kElectronEM:
127 : case kDiElectronEM:
128 0 : fChildSelect[0] = kElectron;
129 0 : break;
130 : case kHardMuons:
131 : case kBSemiMuonic:
132 : case kDSemiMuonic:
133 : case kSemiMuonic:
134 : case kDiMuon:
135 : case kJpsiDiMuon:
136 : case kBJpsiDiMuon:
137 : case kBPsiPrimeDiMuon:
138 : case kPiToMu:
139 : case kKaToMu:
140 : case kWToMuon:
141 : case kWToCharmToMuon:
142 : case kZDiMuon:
143 : case kZDiElectron:
144 0 : fChildSelect[0]=kMuonMinus;
145 0 : break;
146 : case kWToCharm:
147 : break;
148 : case kHadronicD:
149 : case kHadronicDWithout4Bodies:
150 : case kHadronicDWithV0:
151 : case kHadronicDWithout4BodiesWithV0:
152 0 : fChildSelect[0]=kPiPlus;
153 0 : fChildSelect[1]=kKPlus;
154 0 : break;
155 : case kPhiKK:
156 0 : fChildSelect[0]=kKPlus;
157 0 : break;
158 : case kBJpsi:
159 : case kBJpsiUndecayed:
160 0 : fChildSelect[0]= 443;
161 0 : break;
162 : case kChiToJpsiGammaToMuonMuon:
163 0 : fChildSelect[0]= 22;
164 0 : fChildSelect[1]= 13;
165 0 : break;
166 : case kChiToJpsiGammaToElectronElectron:
167 0 : fChildSelect[0]= 22;
168 0 : fChildSelect[1]= 11;
169 0 : break;
170 : case kLambda:
171 0 : fChildSelect[0]= kProton;
172 0 : fChildSelect[1]= 211;
173 0 : break;
174 : case kPsiPrimeJpsiDiElectron:
175 0 : fChildSelect[0]= 211;
176 0 : fChildSelect[1]= 11;
177 0 : break;
178 : case kGammaEM:
179 0 : fChildSelect[0] = kGamma;
180 0 : break;
181 : case kOmega:
182 : case kAll:
183 : case kAllMuonic:
184 : case kNoDecay:
185 : case kNoDecayHeavy:
186 : case kNoDecayBeauty:
187 : case kNeutralPion:
188 : case kBeautyUpgrade:
189 : break;
190 : }
191 :
192 0 : if (fZTarget != 0 && fAProjectile != 0)
193 : {
194 0 : fDyBoost = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) /
195 0 : (Double_t(fZTarget) * Double_t(fAProjectile)));
196 0 : }
197 0 : }
198 :
199 :
200 : Bool_t AliGenMC::ParentSelected(Int_t ip) const
201 : {
202 : // True if particle is in list of parent particles to be selected
203 0 : for (Int_t i=0; i<8; i++)
204 : {
205 0 : if (fParentSelect.At(i) == ip) return kTRUE;
206 : }
207 0 : return kFALSE;
208 0 : }
209 :
210 : Bool_t AliGenMC::ChildSelected(Int_t ip) const
211 : {
212 : // True if particle is in list of decay products to be selected
213 0 : for (Int_t i=0; i<5; i++)
214 : {
215 0 : if (fChildSelect.At(i) == ip) return kTRUE;
216 : }
217 0 : return kFALSE;
218 0 : }
219 :
220 : Bool_t AliGenMC::KinematicSelection(const TParticle *particle, Int_t flag) const
221 : {
222 : // Perform kinematic selection
223 0 : Double_t pz = particle->Pz();
224 0 : Double_t pt = particle->Pt();
225 0 : Double_t p = particle->P();
226 0 : Double_t theta = particle->Theta();
227 0 : Double_t mass = particle->GetCalcMass();
228 0 : Double_t mt2 = pt * pt + mass * mass;
229 0 : Double_t phi = particle->Phi();
230 0 : Double_t e = particle->Energy();
231 :
232 0 : if (e == 0.)
233 0 : e = TMath::Sqrt(p * p + mass * mass);
234 :
235 :
236 : Double_t y, y0;
237 :
238 0 : if (TMath::Abs(pz) < e) {
239 0 : y = 0.5*TMath::Log((e+pz)/(e-pz));
240 0 : } else {
241 : y = 1.e10;
242 : }
243 :
244 0 : if (mt2) {
245 0 : y0 = 0.5*TMath::Log((e+TMath::Abs(pz))*(e+TMath::Abs(pz))/mt2);
246 0 : } else {
247 0 : if (TMath::Abs(y) < 1.e10) {
248 : y0 = y;
249 0 : } else {
250 : y0 = 1.e10;
251 : }
252 : }
253 :
254 0 : y = (pz < 0) ? -y0 : y0;
255 :
256 0 : if (flag == 0) {
257 : //
258 : // Primary particle cuts
259 : //
260 : // transverse momentum cut
261 0 : if (pt > fPtMax || pt < fPtMin) {
262 : // printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
263 0 : return kFALSE;
264 : }
265 : //
266 : // momentum cut
267 0 : if (p > fPMax || p < fPMin) {
268 : // printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
269 0 : return kFALSE;
270 : }
271 : //
272 : // theta cut
273 0 : if (theta > fThetaMax || theta < fThetaMin) {
274 : // printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
275 0 : return kFALSE;
276 : }
277 : //
278 : // rapidity cut
279 0 : if (y > fYMax || y < fYMin) {
280 : // printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
281 0 : return kFALSE;
282 : }
283 : //
284 : // phi cut
285 0 : if (phi > fPhiMax || phi < fPhiMin) {
286 : // printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
287 0 : return kFALSE;
288 : }
289 : } else {
290 : //
291 : // Decay product cuts
292 : //
293 : // transverse momentum cut
294 0 : if (pt > fChildPtMax || pt < fChildPtMin) {
295 : // printf("\n failed pt cut %f %f %f \n",pt,fChildPtMin,fChildPtMax);
296 0 : return kFALSE;
297 : }
298 : //
299 : // momentum cut
300 0 : if (p > fChildPMax || p < fChildPMin) {
301 : // printf("\n failed p cut %f %f %f \n",p,fChildPMin,fChildPMax);
302 0 : return kFALSE;
303 : }
304 : //
305 : // theta cut
306 0 : if (theta > fChildThetaMax || theta < fChildThetaMin) {
307 : // printf("\n failed theta cut %f %f %f \n",theta,fChildThetaMin,fChildThetaMax);
308 0 : return kFALSE;
309 : }
310 : //
311 : // rapidity cut
312 0 : if (y > fChildYMax || y < fChildYMin) {
313 : // printf("\n failed y cut %f %f %f \n",y,fChildYMin,fChildYMax);
314 0 : return kFALSE;
315 : }
316 : //
317 : // phi cut
318 0 : if (phi > fChildPhiMax || phi < fChildPhiMin) {
319 : // printf("\n failed phi cut %f %f %f \n",phi,fChildPhiMin,fChildPhiMax);
320 0 : return kFALSE;
321 : }
322 : }
323 :
324 0 : return kTRUE;
325 0 : }
326 :
327 : Bool_t AliGenMC::CheckAcceptanceGeometry(Int_t np, TClonesArray* particles)
328 : {
329 : // Check the geometrical acceptance for particle.
330 :
331 : Bool_t check ;
332 : Int_t numberOfPdgCodeParticleforAcceptanceCut = 0;
333 : Int_t numberOfAcceptedPdgCodeParticleforAcceptanceCut = 0;
334 : TParticle * particle;
335 : Int_t i;
336 0 : for (i = 0; i < np; i++) {
337 0 : particle = (TParticle *) particles->At(i);
338 0 : if( TMath::Abs( particle->GetPdgCode() ) == TMath::Abs( fPdgCodeParticleforAcceptanceCut ) ) {
339 0 : numberOfPdgCodeParticleforAcceptanceCut++;
340 0 : if (fGeometryAcceptance->Impact(particle)) numberOfAcceptedPdgCodeParticleforAcceptanceCut++;
341 : }
342 : }
343 0 : if ( numberOfAcceptedPdgCodeParticleforAcceptanceCut > (fNumberOfAcceptedParticles-1) )
344 0 : check = kTRUE;
345 : else
346 : check = kFALSE;
347 :
348 0 : return check;
349 : }
350 :
351 : Int_t AliGenMC::CheckPDGCode(Int_t pdgcode) const
352 : {
353 : //
354 : // If the particle is in a diffractive state, then take action accordingly
355 0 : switch (pdgcode) {
356 : case 91:
357 0 : return 92;
358 : case 110:
359 : //rho_diff0 -- difficult to translate, return rho0
360 0 : return 113;
361 : case 210:
362 : //pi_diffr+ -- change to pi+
363 0 : return 211;
364 : case 220:
365 : //omega_di0 -- change to omega0
366 0 : return 223;
367 : case 330:
368 : //phi_diff0 -- return phi0
369 0 : return 333;
370 : case 440:
371 : //J/psi_di0 -- return J/psi
372 0 : return 443;
373 : case 2110:
374 : //n_diffr -- return neutron
375 0 : return 2112;
376 : case 2210:
377 : //p_diffr+ -- return proton
378 0 : return 2212;
379 : }
380 : //non diffractive state -- return code unchanged
381 0 : return pdgcode;
382 0 : }
383 :
384 : void AliGenMC::Boost()
385 : {
386 : //
387 : // Boost cms into LHC lab frame
388 : //
389 :
390 0 : Double_t beta = TMath::TanH(fDyBoost);
391 0 : Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta));
392 0 : Double_t gb = gamma * beta;
393 :
394 : // printf("\n Boosting particles to lab frame %f %f %f", fDyBoost, beta, gamma);
395 :
396 : Int_t i;
397 0 : Int_t np = fParticles.GetEntriesFast();
398 0 : for (i = 0; i < np; i++)
399 : {
400 0 : TParticle* iparticle = (TParticle*) fParticles.At(i);
401 :
402 0 : Double_t e = iparticle->Energy();
403 0 : Double_t px = iparticle->Px();
404 0 : Double_t py = iparticle->Py();
405 0 : Double_t pz = iparticle->Pz();
406 :
407 0 : Double_t eb = gamma * e - gb * pz;
408 0 : Double_t pzb = -gb * e + gamma * pz;
409 :
410 0 : iparticle->SetMomentum(px, py, pzb, eb);
411 : }
412 0 : }
413 :
414 : void AliGenMC::BeamCrossAngle()
415 : {
416 : // Applies a boost in the y-direction in order to take into account the
417 : // beam crossing angle
418 :
419 : Double_t thetaPr0, pyPr2, pzPr2;
420 0 : TVector3 beta;
421 :
422 0 : thetaPr0 = fXingAngleY / 2.;
423 :
424 : // Momentum of the CMS system
425 0 : pyPr2 = TMath::Sqrt(fEnergyCMS * fEnergyCMS/ 4 - 0.938 * 0.938) * TMath::Sin(thetaPr0);
426 0 : pzPr2 = TMath::Sqrt(fEnergyCMS * fEnergyCMS/ 4 - 0.938 * 0.938) * TMath::Cos(thetaPr0);
427 :
428 0 : TLorentzVector proj1Vect, proj2Vect, projVect;
429 0 : proj1Vect.SetPxPyPzE(0., pyPr2, pzPr2, fEnergyCMS/2);
430 0 : proj2Vect.SetPxPyPzE(0., pyPr2,-pzPr2, fEnergyCMS/2);
431 0 : projVect = proj1Vect + proj2Vect;
432 0 : beta=(1. / projVect.E()) * (projVect.Vect());
433 :
434 : Int_t i;
435 0 : Int_t np = fParticles.GetEntriesFast();
436 0 : for (i = 0; i < np; i++)
437 : {
438 0 : TParticle* iparticle = (TParticle*) fParticles.At(i);
439 :
440 0 : Double_t e = iparticle->Energy();
441 0 : Double_t px = iparticle->Px();
442 0 : Double_t py = iparticle->Py();
443 0 : Double_t pz = iparticle->Pz();
444 :
445 0 : TLorentzVector partIn;
446 0 : partIn.SetPxPyPzE(px,py,pz,e);
447 0 : partIn.Boost(beta);
448 0 : iparticle->SetMomentum(partIn.Px(),partIn.Py(),partIn.Pz(),partIn.E());
449 0 : }
450 0 : }
451 :
452 :
453 : void AliGenMC::AddHeader(AliGenEventHeader* header)
454 : {
455 : // Passes header either to the container or to gAlice
456 0 : if (fContainer) {
457 0 : header->SetName(fName);
458 0 : fContainer->AddHeader(header);
459 0 : } else {
460 0 : if (gAlice) gAlice->SetGenEventHeader(header);
461 : }
462 0 : }
|
