| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file eventgenerator/pythia/decayer6/src/G4Pythia6Decayer.cc
28 /// \brief Implementation of the G4Pythia6Decayer class
29
30 // ----------------------------------------------------------------------------
31 // According to TPythia6Decayer class in Root:
32 // http://root.cern.ch/
33 // see http://root.cern.ch/root/License.html
34 // ----------------------------------------------------------------------------
35
36 #include "G4Pythia6Decayer.hh"
37
38 #include "Pythia6.hh"
39
40 #include "G4DecayProducts.hh"
41 #include "G4DecayTable.hh"
42 #include "G4DynamicParticle.hh"
43 #include "G4ParticleTable.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Track.hh"
46
47 #include <CLHEP/Vector/LorentzVector.h>
48 #include <cmath>
49
50 const EDecayType G4Pythia6Decayer::fgkDefaultDecayType = kAll;
51
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53
54 G4Pythia6Decayer::G4Pythia6Decayer()
55 : G4VExtDecayer("G4Pythia6Decayer"),
56 fMessenger(this),
57 fVerboseLevel(0),
58 fDecayType(fgkDefaultDecayType),
59 fDecayProductsArray(0)
60 {
61 /// Standard constructor
62
63 fDecayProductsArray = new ParticleVector();
64
65 ForceDecay(fDecayType);
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70 G4Pythia6Decayer::~G4Pythia6Decayer()
71 {
72 /// Destructor
73
74 delete fDecayProductsArray;
75 }
76
77 //
78 // private methods
79 //
80
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82
83 G4ParticleDefinition* G4Pythia6Decayer::GetParticleDefinition(const Pythia6Particle* particle,
84 G4bool warn) const
85 {
86 /// Return G4 particle definition for given TParticle
87
88 // get particle definition from G4ParticleTable
89 G4int pdgEncoding = particle->fKF;
90 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
91 G4ParticleDefinition* particleDefinition = 0;
92 if (pdgEncoding != 0) particleDefinition = particleTable->FindParticle(pdgEncoding);
93
94 if (particleDefinition == 0 && warn) {
95 G4cerr << "G4Pythia6Decayer: GetParticleDefinition: " << std::endl
96 << "G4ParticleTable::FindParticle() for particle with PDG = " << pdgEncoding
97 << " failed." << std::endl;
98 }
99
100 return particleDefinition;
101 }
102
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104
105 G4DynamicParticle* G4Pythia6Decayer::CreateDynamicParticle(const Pythia6Particle* particle) const
106 {
107 /// Create G4DynamicParticle.
108
109 // get particle properties
110 const G4ParticleDefinition* particleDefinition = GetParticleDefinition(particle);
111 if (!particleDefinition) return 0;
112
113 G4ThreeVector momentum = GetParticleMomentum(particle);
114
115 // create G4DynamicParticle
116 G4DynamicParticle* dynamicParticle = new G4DynamicParticle(particleDefinition, momentum);
117
118 return dynamicParticle;
119 }
120
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122
123 G4ThreeVector G4Pythia6Decayer::GetParticlePosition(const Pythia6Particle* particle) const
124 {
125 /// Return particle vertex position.
126
127 G4ThreeVector position =
128 G4ThreeVector(particle->fVx * cm, particle->fVy * cm, particle->fVz * cm);
129 return position;
130 }
131
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133
134 G4ThreeVector G4Pythia6Decayer::GetParticleMomentum(const Pythia6Particle* particle) const
135 {
136 /// Return particle momentum.
137
138 G4ThreeVector momentum =
139 G4ThreeVector(particle->fPx * GeV, particle->fPy * GeV, particle->fPz * GeV);
140 return momentum;
141 }
142
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144
145 G4int G4Pythia6Decayer::CountProducts(G4int channel, G4int particle)
146 {
147 /// Count number of decay products
148
149 G4int np = 0;
150 for (G4int i = 1; i <= 5; i++)
151 if (std::abs(Pythia6::Instance()->GetKFDP(channel, i)) == particle) np++;
152 return np;
153 }
154
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156
157 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int product, G4int mult)
158 {
159 /// Force decay of particle into products with multiplicity mult
160
161 Pythia6* pythia6 = Pythia6::Instance();
162
163 G4int kc = pythia6->Pycomp(particle);
164 pythia6->SetMDCY(kc, 1, 1);
165
166 G4int ifirst = pythia6->GetMDCY(kc, 2);
167 G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
168
169 //
170 // Loop over decay channels
171 for (G4int channel = ifirst; channel <= ilast; channel++) {
172 if (CountProducts(channel, product) >= mult) {
173 pythia6->SetMDME(channel, 1, 1);
174 }
175 else {
176 pythia6->SetMDME(channel, 1, 0);
177 }
178 }
179 }
180
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182
183 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int* products, G4int* mult, G4int npart)
184 {
185 /// Force decay of particle into products with multiplicity mult
186
187 Pythia6* pythia6 = Pythia6::Instance();
188
189 G4int kc = pythia6->Pycomp(particle);
190 pythia6->SetMDCY(kc, 1, 1);
191 G4int ifirst = pythia6->GetMDCY(kc, 2);
192 G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
193 //
194 // Loop over decay channels
195 for (G4int channel = ifirst; channel <= ilast; channel++) {
196 G4int nprod = 0;
197 for (G4int i = 0; i < npart; i++)
198 nprod += (CountProducts(channel, products[i]) >= mult[i]);
199 if (nprod)
200 pythia6->SetMDME(channel, 1, 1);
201 else {
202 pythia6->SetMDME(channel, 1, 0);
203 }
204 }
205 }
206
207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
208
209 void G4Pythia6Decayer::ForceHadronicD()
210 {
211 /// Force golden D decay modes
212
213 const G4int kNHadrons = 4;
214 G4int channel;
215 G4int hadron[kNHadrons] = {411, 421, 431, 4112};
216
217 // for D+ -> K0* (-> K- pi+) pi+
218 G4int iKstar0 = 313;
219 G4int iKstarbar0 = -313;
220 G4int iKPlus = 321;
221 G4int iKMinus = -321;
222 G4int iPiPlus = 211;
223 G4int iPiMinus = -211;
224
225 G4int products[2] = {iKPlus, iPiMinus}, mult[2] = {1, 1};
226 ForceParticleDecay(iKstar0, products, mult, 2);
227
228 // for Ds -> Phi pi+
229 G4int iPhi = 333;
230 ForceParticleDecay(iPhi, iKPlus, 2); // Phi->K+K-
231
232 G4int decayP1[kNHadrons][3] = {
233 {iKMinus, iPiPlus, iPiPlus}, {iKMinus, iPiPlus, 0}, {iKPlus, iKstarbar0, 0}, {-1, -1, -1}};
234 G4int decayP2[kNHadrons][3] = {
235 {iKstarbar0, iPiPlus, 0}, {-1, -1, -1}, {iPhi, iPiPlus, 0}, {-1, -1, -1}};
236
237 Pythia6* pythia6 = Pythia6::Instance();
238 for (G4int ihadron = 0; ihadron < kNHadrons; ihadron++) {
239 G4int kc = pythia6->Pycomp(hadron[ihadron]);
240 pythia6->SetMDCY(kc, 1, 1);
241 G4int ifirst = pythia6->GetMDCY(kc, 2);
242 G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
243
244 for (channel = ifirst; channel <= ilast; channel++) {
245 if ((pythia6->GetKFDP(channel, 1) == decayP1[ihadron][0]
246 && pythia6->GetKFDP(channel, 2) == decayP1[ihadron][1]
247 && pythia6->GetKFDP(channel, 3) == decayP1[ihadron][2]
248 && pythia6->GetKFDP(channel, 4) == 0)
249 || (pythia6->GetKFDP(channel, 1) == decayP2[ihadron][0]
250 && pythia6->GetKFDP(channel, 2) == decayP2[ihadron][1]
251 && pythia6->GetKFDP(channel, 3) == decayP2[ihadron][2]
252 && pythia6->GetKFDP(channel, 4) == 0))
253 {
254 pythia6->SetMDME(channel, 1, 1);
255 }
256 else {
257 pythia6->SetMDME(channel, 1, 0);
258 } // selected channel ?
259 } // decay channels
260 } // hadrons
261 }
262
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
264
265 void G4Pythia6Decayer::ForceOmega()
266 {
267 /// Force Omega -> Lambda K- Decay
268
269 Pythia6* pythia6 = Pythia6::Instance();
270
271 G4int iLambda0 = 3122;
272 G4int iKMinus = -321;
273
274 G4int kc = pythia6->Pycomp(3334);
275 pythia6->SetMDCY(kc, 1, 1);
276 G4int ifirst = pythia6->GetMDCY(kc, 2);
277 G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
278
279 for (G4int channel = ifirst; channel <= ilast; channel++) {
280 if (pythia6->GetKFDP(channel, 1) == iLambda0 && pythia6->GetKFDP(channel, 2) == iKMinus
281 && pythia6->GetKFDP(channel, 3) == 0)
282 pythia6->SetMDME(channel, 1, 1);
283 else
284 pythia6->SetMDME(channel, 1, 0);
285 // selected channel ?
286 } // decay channels
287 }
288
289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
290
291 void G4Pythia6Decayer::ForceDecay(EDecayType decayType)
292 {
293 /// Force a particle decay mode
294
295 Pythia6::Instance()->SetMSTJ(21, 2);
296
297 if (fDecayType == kNoDecayHeavy) return;
298
299 //
300 // select mode
301 G4int products[3];
302 G4int mult[3];
303
304 switch (decayType) {
305 case kHardMuons:
306 products[0] = 13;
307 products[1] = 443;
308 products[2] = 100443;
309 mult[0] = 1;
310 mult[1] = 1;
311 mult[2] = 1;
312 ForceParticleDecay(511, products, mult, 3);
313 ForceParticleDecay(521, products, mult, 3);
314 ForceParticleDecay(531, products, mult, 3);
315 ForceParticleDecay(5122, products, mult, 3);
316 ForceParticleDecay(5132, products, mult, 3);
317 ForceParticleDecay(5232, products, mult, 3);
318 ForceParticleDecay(5332, products, mult, 3);
319 ForceParticleDecay(100443, 443, 1); // Psi' -> J/Psi X
320 ForceParticleDecay(443, 13, 2); // J/Psi -> mu+ mu-
321
322 ForceParticleDecay(411, 13, 1); // D+/-
323 ForceParticleDecay(421, 13, 1); // D0
324 ForceParticleDecay(431, 13, 1); // D_s
325 ForceParticleDecay(4122, 13, 1); // Lambda_c
326 ForceParticleDecay(4132, 13, 1); // Xsi_c
327 ForceParticleDecay(4232, 13, 1); // Sigma_c
328 ForceParticleDecay(4332, 13, 1); // Omega_c
329 break;
330
331 case kSemiMuonic:
332 ForceParticleDecay(411, 13, 1); // D+/-
333 ForceParticleDecay(421, 13, 1); // D0
334 ForceParticleDecay(431, 13, 1); // D_s
335 ForceParticleDecay(4122, 13, 1); // Lambda_c
336 ForceParticleDecay(4132, 13, 1); // Xsi_c
337 ForceParticleDecay(4232, 13, 1); // Sigma_c
338 ForceParticleDecay(4332, 13, 1); // Omega_c
339 ForceParticleDecay(511, 13, 1); // B0
340 ForceParticleDecay(521, 13, 1); // B+/-
341 ForceParticleDecay(531, 13, 1); // B_s
342 ForceParticleDecay(5122, 13, 1); // Lambda_b
343 ForceParticleDecay(5132, 13, 1); // Xsi_b
344 ForceParticleDecay(5232, 13, 1); // Sigma_b
345 ForceParticleDecay(5332, 13, 1); // Omega_b
346 break;
347
348 case kDiMuon:
349 ForceParticleDecay(113, 13, 2); // rho
350 ForceParticleDecay(221, 13, 2); // eta
351 ForceParticleDecay(223, 13, 2); // omega
352 ForceParticleDecay(333, 13, 2); // phi
353 ForceParticleDecay(443, 13, 2); // J/Psi
354 ForceParticleDecay(100443, 13, 2); // Psi'
355 ForceParticleDecay(553, 13, 2); // Upsilon
356 ForceParticleDecay(100553, 13, 2); // Upsilon'
357 ForceParticleDecay(200553, 13, 2); // Upsilon''
358 break;
359
360 case kSemiElectronic:
361 ForceParticleDecay(411, 11, 1); // D+/-
362 ForceParticleDecay(421, 11, 1); // D0
363 ForceParticleDecay(431, 11, 1); // D_s
364 ForceParticleDecay(4122, 11, 1); // Lambda_c
365 ForceParticleDecay(4132, 11, 1); // Xsi_c
366 ForceParticleDecay(4232, 11, 1); // Sigma_c
367 ForceParticleDecay(4332, 11, 1); // Omega_c
368 ForceParticleDecay(511, 11, 1); // B0
369 ForceParticleDecay(521, 11, 1); // B+/-
370 ForceParticleDecay(531, 11, 1); // B_s
371 ForceParticleDecay(5122, 11, 1); // Lambda_b
372 ForceParticleDecay(5132, 11, 1); // Xsi_b
373 ForceParticleDecay(5232, 11, 1); // Sigma_b
374 ForceParticleDecay(5332, 11, 1); // Omega_b
375 break;
376
377 case kDiElectron:
378 ForceParticleDecay(113, 11, 2); // rho
379 ForceParticleDecay(333, 11, 2); // phi
380 ForceParticleDecay(221, 11, 2); // eta
381 ForceParticleDecay(223, 11, 2); // omega
382 ForceParticleDecay(443, 11, 2); // J/Psi
383 ForceParticleDecay(100443, 11, 2); // Psi'
384 ForceParticleDecay(553, 11, 2); // Upsilon
385 ForceParticleDecay(100553, 11, 2); // Upsilon'
386 ForceParticleDecay(200553, 11, 2); // Upsilon''
387 break;
388
389 case kBJpsiDiMuon:
390
391 products[0] = 443;
392 products[1] = 100443;
393 mult[0] = 1;
394 mult[1] = 1;
395
396 ForceParticleDecay(511, products, mult, 2); // B0 -> J/Psi (Psi') X
397 ForceParticleDecay(521, products, mult, 2); // B+/- -> J/Psi (Psi') X
398 ForceParticleDecay(531, products, mult, 2); // B_s -> J/Psi (Psi') X
399 ForceParticleDecay(5122, products, mult, 2); // Lambda_b -> J/Psi (Psi')X
400 ForceParticleDecay(100443, 443, 1); // Psi' -> J/Psi X
401 ForceParticleDecay(443, 13, 2); // J/Psi -> mu+ mu-
402 break;
403
404 case kBPsiPrimeDiMuon:
405 ForceParticleDecay(511, 100443, 1); // B0
406 ForceParticleDecay(521, 100443, 1); // B+/-
407 ForceParticleDecay(531, 100443, 1); // B_s
408 ForceParticleDecay(5122, 100443, 1); // Lambda_b
409 ForceParticleDecay(100443, 13, 2); // Psi'
410 break;
411
412 case kBJpsiDiElectron:
413 ForceParticleDecay(511, 443, 1); // B0
414 ForceParticleDecay(521, 443, 1); // B+/-
415 ForceParticleDecay(531, 443, 1); // B_s
416 ForceParticleDecay(5122, 443, 1); // Lambda_b
417 ForceParticleDecay(443, 11, 2); // J/Psi
418 break;
419
420 case kBJpsi:
421 ForceParticleDecay(511, 443, 1); // B0
422 ForceParticleDecay(521, 443, 1); // B+/-
423 ForceParticleDecay(531, 443, 1); // B_s
424 ForceParticleDecay(5122, 443, 1); // Lambda_b
425 break;
426
427 case kBPsiPrimeDiElectron:
428 ForceParticleDecay(511, 100443, 1); // B0
429 ForceParticleDecay(521, 100443, 1); // B+/-
430 ForceParticleDecay(531, 100443, 1); // B_s
431 ForceParticleDecay(5122, 100443, 1); // Lambda_b
432 ForceParticleDecay(100443, 11, 2); // Psi'
433 break;
434
435 case kPiToMu:
436 ForceParticleDecay(211, 13, 1); // pi->mu
437 break;
438
439 case kKaToMu:
440 ForceParticleDecay(321, 13, 1); // K->mu
441 break;
442
443 case kWToMuon:
444 ForceParticleDecay(24, 13, 1); // W -> mu
445 break;
446
447 case kWToCharm:
448 ForceParticleDecay(24, 4, 1); // W -> c
449 break;
450
451 case kWToCharmToMuon:
452 ForceParticleDecay(24, 4, 1); // W -> c
453 ForceParticleDecay(411, 13, 1); // D+/- -> mu
454 ForceParticleDecay(421, 13, 1); // D0 -> mu
455 ForceParticleDecay(431, 13, 1); // D_s -> mu
456 ForceParticleDecay(4122, 13, 1); // Lambda_c
457 ForceParticleDecay(4132, 13, 1); // Xsi_c
458 ForceParticleDecay(4232, 13, 1); // Sigma_c
459 ForceParticleDecay(4332, 13, 1); // Omega_c
460 break;
461
462 case kZDiMuon:
463 ForceParticleDecay(23, 13, 2); // Z -> mu+ mu-
464 break;
465
466 case kHadronicD:
467 ForceHadronicD();
468 break;
469
470 case kPhiKK:
471 ForceParticleDecay(333, 321, 2); // Phi->K+K-
472 break;
473
474 case kOmega:
475 ForceOmega();
476
477 case kAll:
478 break;
479
480 case kNoDecay:
481 Pythia6::Instance()->SetMSTJ(21, 0);
482 break;
483
484 case kNoDecayHeavy:
485 break;
486
487 case kMaxDecay:
488 break;
489 }
490 }
491
492 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
493
494 void G4Pythia6Decayer::Decay(G4int pdg, const CLHEP::HepLorentzVector& p)
495 {
496 /// Decay a particle of type IDPART (PDG code) and momentum P.
497
498 Pythia6::Instance()->Py1ent(0, pdg, p.e(), p.theta(), p.phi());
499 }
500
501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
502
503 G4int G4Pythia6Decayer::ImportParticles(ParticleVector* particles)
504 {
505 /// Get the decay products into the passed PARTICLES vector
506
507 return Pythia6::Instance()->ImportParticles(particles, "All");
508 }
509
510 //
511 // public methods
512 //
513
514 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
515
516 G4DecayProducts* G4Pythia6Decayer::ImportDecayProducts(const G4Track& track)
517 {
518 /// Import decay products
519
520 // get particle momentum
521 G4ThreeVector momentum = track.GetMomentum();
522 G4double etot = track.GetDynamicParticle()->GetTotalEnergy();
523 ;
524 CLHEP::HepLorentzVector p;
525 p[0] = momentum.x() / GeV;
526 p[1] = momentum.y() / GeV;
527 p[2] = momentum.z() / GeV;
528 p[3] = etot / GeV;
529
530 // get particle PDG
531 // ask G4Pythia6Decayer to get PDG encoding
532 // (in order to get PDG from extended TDatabasePDG
533 // in case the standard PDG code is not defined)
534 G4ParticleDefinition* particleDef = track.GetDefinition();
535 G4int pdgEncoding = particleDef->GetPDGEncoding();
536
537 // let Pythia6Decayer decay the particle
538 // and import the decay products
539 Decay(pdgEncoding, p);
540 G4int nofParticles = ImportParticles(fDecayProductsArray);
541
542 if (fVerboseLevel > 0) {
543 G4cout << "nofParticles: " << nofParticles << G4endl;
544 }
545
546 // convert decay products Pythia6Particle type
547 // to G4DecayProducts
548 G4DecayProducts* decayProducts = new G4DecayProducts(*(track.GetDynamicParticle()));
549
550 G4int counter = 0;
551 for (G4int i = 0; i < nofParticles; i++) {
552 // get particle from ParticleVector
553 Pythia6Particle* particle = (*fDecayProductsArray)[i];
554
555 G4int status = particle->fKS;
556 G4int pdg = particle->fKF;
557 if (status > 0 && status < 11 && std::abs(pdg) != 12 && std::abs(pdg) != 14
558 && std::abs(pdg) != 16)
559 {
560 // pass to tracking final particles only;
561 // skip neutrinos
562
563 if (fVerboseLevel > 0) {
564 G4cout << " " << i << "th particle PDG: " << pdg << " ";
565 }
566
567 // create G4DynamicParticle
568 G4DynamicParticle* dynamicParticle = CreateDynamicParticle(particle);
569
570 if (dynamicParticle) {
571 if (fVerboseLevel > 0) {
572 G4cout << " G4 particle name: " << dynamicParticle->GetDefinition()->GetParticleName()
573 << G4endl;
574 }
575
576 // add dynamicParticle to decayProducts
577 decayProducts->PushProducts(dynamicParticle);
578
579 counter++;
580 }
581 }
582 }
583 if (fVerboseLevel > 0) {
584 G4cout << "nofParticles for tracking: " << counter << G4endl;
585 }
586
587 return decayProducts;
588 }
589
590 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
591
592 void G4Pythia6Decayer::ForceDecayType(EDecayType decayType)
593 {
594 /// Force a given decay type
595
596 // Do nothing if the decay type is not different from current one
597 if (decayType == fDecayType) return;
598
599 fDecayType = decayType;
600 ForceDecay(fDecayType);
601 }
602