40 #include <Geant4/G4DynamicParticle.hh>
41 #include <Geant4/G4DecayProducts.hh>
42 #include <Geant4/G4DecayTable.hh>
43 #include <Geant4/G4ParticleTable.hh>
44 #include <Geant4/G4Track.hh>
45 #include <Geant4/G4SystemOfUnits.hh>
47 #include <CLHEP/Vector/LorentzVector.h>
56 : G4VExtDecayer(
"G4Pythia6Decayer"),
59 fDecayType(fgkDefaultDecayType),
60 fDecayProductsArray(0)
66 ForceDecay(fDecayType);
75 delete fDecayProductsArray;
85 G4ParticleDefinition* G4Pythia6Decayer::
86 GetParticleDefinition(
const Pythia6Particle* particle, G4bool warn)
const
91 G4int pdgEncoding = particle->
fKF;
92 G4ParticleTable* particleTable
93 = G4ParticleTable::GetParticleTable();
94 G4ParticleDefinition* particleDefinition = 0;
96 particleDefinition = particleTable->FindParticle(pdgEncoding);
98 if ( particleDefinition == 0 && warn) {
100 <<
"G4Pythia6Decayer: GetParticleDefinition: " << std::endl
101 <<
"G4ParticleTable::FindParticle() for particle with PDG = "
103 <<
" failed." << std::endl;
106 return particleDefinition;
112 G4Pythia6Decayer::CreateDynamicParticle(
const Pythia6Particle* particle)
const
117 const G4ParticleDefinition* particleDefinition
118 = GetParticleDefinition(particle);
119 if ( ! particleDefinition )
return 0;
121 G4ThreeVector momentum = GetParticleMomentum(particle);
124 G4DynamicParticle* dynamicParticle
125 =
new G4DynamicParticle(particleDefinition, momentum);
127 return dynamicParticle;
132 G4ThreeVector G4Pythia6Decayer::GetParticlePosition(
137 G4ThreeVector position
138 = G4ThreeVector(particle->
fVx * cm,
146 G4ThreeVector G4Pythia6Decayer::GetParticleMomentum(
151 G4ThreeVector momentum
152 = G4ThreeVector(particle->
fPx * GeV,
154 particle->
fPz * GeV);
160 G4int G4Pythia6Decayer::CountProducts(G4int channel, G4int particle)
165 for ( G4int i=1; i<=5; i++ )
174 G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int product, G4int mult)
180 G4int kc = pythia6->
Pycomp(particle);
183 G4int ifirst = pythia6->
GetMDCY(kc,2);
184 G4int ilast = ifirst + pythia6->
GetMDCY(kc,3)-1;
188 for (G4int channel= ifirst; channel <= ilast; channel++) {
189 if (CountProducts(channel,product) >= mult) {
199 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int* products,
200 G4int* mult, G4int npart)
206 G4int kc = pythia6->
Pycomp(particle);
208 G4int ifirst = pythia6->
GetMDCY(kc,2);
209 G4int ilast = ifirst+pythia6->
GetMDCY(kc,3)-1;
212 for (G4int channel = ifirst; channel <= ilast; channel++) {
214 for (G4int i = 0; i < npart; i++)
215 nprod += (CountProducts(channel, products[i]) >= mult[i]);
226 void G4Pythia6Decayer::ForceHadronicD()
230 const G4int kNHadrons = 4;
232 G4int hadron[kNHadrons] = {411, 421, 431, 4112};
236 G4int iKstarbar0 = -313;
238 G4int iKMinus = -321;
240 G4int iPiMinus = -211;
242 G4int products[2] = {iKPlus, iPiMinus}, mult[2] = {1, 1};
243 ForceParticleDecay(iKstar0, products, mult, 2);
247 ForceParticleDecay(iPhi,iKPlus,2);
249 G4int decayP1[kNHadrons][3] = {
250 {iKMinus, iPiPlus, iPiPlus},
251 {iKMinus, iPiPlus, 0 },
252 {iKPlus , iKstarbar0, 0 },
255 G4int decayP2[kNHadrons][3] = {
256 {iKstarbar0, iPiPlus, 0 },
258 {iPhi , iPiPlus, 0 },
263 for ( G4int ihadron = 0; ihadron < kNHadrons; ihadron++ ) {
264 G4int kc = pythia6->
Pycomp(hadron[ihadron]);
266 G4int ifirst = pythia6->
GetMDCY(kc,2);
267 G4int ilast = ifirst + pythia6->
GetMDCY(kc,3)-1;
269 for (channel = ifirst; channel <= ilast; channel++) {
270 if ((pythia6->
GetKFDP(channel,1) == decayP1[ihadron][0] &&
271 pythia6->
GetKFDP(channel,2) == decayP1[ihadron][1] &&
272 pythia6->
GetKFDP(channel,3) == decayP1[ihadron][2] &&
273 pythia6->
GetKFDP(channel,4) == 0) ||
274 (pythia6->
GetKFDP(channel,1) == decayP2[ihadron][0] &&
275 pythia6->
GetKFDP(channel,2) == decayP2[ihadron][1] &&
276 pythia6->
GetKFDP(channel,3) == decayP2[ihadron][2] &&
277 pythia6->
GetKFDP(channel,4) == 0)) {
288 void G4Pythia6Decayer::ForceOmega()
294 G4int iLambda0 = 3122;
295 G4int iKMinus = -321;
297 G4int kc = pythia6->
Pycomp(3334);
299 G4int ifirst = pythia6->
GetMDCY(kc,2);
300 G4int ilast = ifirst + pythia6->
GetMDCY(kc,3)-1;
302 for (G4int channel = ifirst; channel <= ilast; channel++) {
303 if (pythia6->
GetKFDP(channel,1) == iLambda0 &&
304 pythia6->
GetKFDP(channel,2) == iKMinus &&
305 pythia6->
GetKFDP(channel,3) == 0)
315 void G4Pythia6Decayer::ForceDecay(
EDecayType decayType)
328 switch ( decayType ) {
333 products[2] = 100443;
337 ForceParticleDecay( 511, products, mult, 3);
338 ForceParticleDecay( 521, products, mult, 3);
339 ForceParticleDecay( 531, products, mult, 3);
340 ForceParticleDecay( 5122, products, mult, 3);
341 ForceParticleDecay( 5132, products, mult, 3);
342 ForceParticleDecay( 5232, products, mult, 3);
343 ForceParticleDecay( 5332, products, mult, 3);
344 ForceParticleDecay( 100443, 443, 1);
345 ForceParticleDecay( 443, 13, 2);
347 ForceParticleDecay( 411,13,1);
348 ForceParticleDecay( 421,13,1);
349 ForceParticleDecay( 431,13,1);
350 ForceParticleDecay( 4122,13,1);
351 ForceParticleDecay( 4132,13,1);
352 ForceParticleDecay( 4232,13,1);
353 ForceParticleDecay( 4332,13,1);
357 ForceParticleDecay( 411,13,1);
358 ForceParticleDecay( 421,13,1);
359 ForceParticleDecay( 431,13,1);
360 ForceParticleDecay( 4122,13,1);
361 ForceParticleDecay( 4132,13,1);
362 ForceParticleDecay( 4232,13,1);
363 ForceParticleDecay( 4332,13,1);
364 ForceParticleDecay( 511,13,1);
365 ForceParticleDecay( 521,13,1);
366 ForceParticleDecay( 531,13,1);
367 ForceParticleDecay( 5122,13,1);
368 ForceParticleDecay( 5132,13,1);
369 ForceParticleDecay( 5232,13,1);
370 ForceParticleDecay( 5332,13,1);
374 ForceParticleDecay( 113,13,2);
375 ForceParticleDecay( 221,13,2);
376 ForceParticleDecay( 223,13,2);
377 ForceParticleDecay( 333,13,2);
378 ForceParticleDecay( 443,13,2);
379 ForceParticleDecay(100443,13,2);
380 ForceParticleDecay( 553,13,2);
381 ForceParticleDecay(100553,13,2);
382 ForceParticleDecay(200553,13,2);
386 ForceParticleDecay( 411,11,1);
387 ForceParticleDecay( 421,11,1);
388 ForceParticleDecay( 431,11,1);
389 ForceParticleDecay( 4122,11,1);
390 ForceParticleDecay( 4132,11,1);
391 ForceParticleDecay( 4232,11,1);
392 ForceParticleDecay( 4332,11,1);
393 ForceParticleDecay( 511,11,1);
394 ForceParticleDecay( 521,11,1);
395 ForceParticleDecay( 531,11,1);
396 ForceParticleDecay( 5122,11,1);
397 ForceParticleDecay( 5132,11,1);
398 ForceParticleDecay( 5232,11,1);
399 ForceParticleDecay( 5332,11,1);
403 ForceParticleDecay( 113,11,2);
404 ForceParticleDecay( 333,11,2);
405 ForceParticleDecay( 221,11,2);
406 ForceParticleDecay( 223,11,2);
407 ForceParticleDecay( 443,11,2);
408 ForceParticleDecay(100443,11,2);
409 ForceParticleDecay( 553,11,2);
410 ForceParticleDecay(100553,11,2);
411 ForceParticleDecay(200553,11,2);
417 products[1] = 100443;
421 ForceParticleDecay( 511, products, mult, 2);
422 ForceParticleDecay( 521, products, mult, 2);
423 ForceParticleDecay( 531, products, mult, 2);
424 ForceParticleDecay( 5122, products, mult, 2);
425 ForceParticleDecay( 100443, 443, 1);
426 ForceParticleDecay( 443,13,2);
430 ForceParticleDecay( 511,100443,1);
431 ForceParticleDecay( 521,100443,1);
432 ForceParticleDecay( 531,100443,1);
433 ForceParticleDecay( 5122,100443,1);
434 ForceParticleDecay(100443,13,2);
438 ForceParticleDecay( 511,443,1);
439 ForceParticleDecay( 521,443,1);
440 ForceParticleDecay( 531,443,1);
441 ForceParticleDecay( 5122,443,1);
442 ForceParticleDecay( 443,11,2);
446 ForceParticleDecay( 511,443,1);
447 ForceParticleDecay( 521,443,1);
448 ForceParticleDecay( 531,443,1);
449 ForceParticleDecay( 5122,443,1);
453 ForceParticleDecay( 511,100443,1);
454 ForceParticleDecay( 521,100443,1);
455 ForceParticleDecay( 531,100443,1);
456 ForceParticleDecay( 5122,100443,1);
457 ForceParticleDecay(100443,11,2);
461 ForceParticleDecay(211,13,1);
465 ForceParticleDecay(321,13,1);
469 ForceParticleDecay( 24, 13,1);
473 ForceParticleDecay( 24, 4,1);
477 ForceParticleDecay( 24, 4,1);
478 ForceParticleDecay( 411,13,1);
479 ForceParticleDecay( 421,13,1);
480 ForceParticleDecay( 431,13,1);
481 ForceParticleDecay( 4122,13,1);
482 ForceParticleDecay( 4132,13,1);
483 ForceParticleDecay( 4232,13,1);
484 ForceParticleDecay( 4332,13,1);
488 ForceParticleDecay( 23, 13,2);
496 ForceParticleDecay(333,321,2);
517 void G4Pythia6Decayer::Decay(G4int pdg,
const CLHEP::HepLorentzVector& p)
526 G4int G4Pythia6Decayer::ImportParticles(
ParticleVector* particles)
544 G4ThreeVector momentum = track.GetMomentum();
545 G4double etot = track.GetDynamicParticle()->GetTotalEnergy();;
546 CLHEP::HepLorentzVector p;
547 p[0] = momentum.x() / GeV;
548 p[1] = momentum.y() / GeV;
549 p[2] = momentum.z() / GeV;
556 G4ParticleDefinition* particleDef = track.GetDefinition();
557 G4int pdgEncoding = particleDef->GetPDGEncoding();
561 Decay(pdgEncoding, p);
562 G4int nofParticles = ImportParticles(fDecayProductsArray);
564 if ( fVerboseLevel > 0 ) {
565 G4cout <<
"nofParticles: " << nofParticles << G4endl;
570 G4DecayProducts* decayProducts
571 =
new G4DecayProducts(*(track.GetDynamicParticle()));
574 for (G4int i=0; i<nofParticles; i++) {
579 G4int status = particle->
fKS;
580 G4int pdg = particle->
fKF;
581 if ( status>0 && status<11 &&
582 std::abs(pdg)!=12 && std::abs(pdg)!=14 && std::abs(pdg)!=16 ) {
586 if ( fVerboseLevel > 0 ) {
587 G4cout <<
" " << i <<
"th particle PDG: " << pdg <<
" ";
591 G4DynamicParticle* dynamicParticle
592 = CreateDynamicParticle(particle);
594 if (dynamicParticle) {
596 if ( fVerboseLevel > 0 ) {
597 G4cout <<
" G4 particle name: "
598 << dynamicParticle->GetDefinition()->GetParticleName()
603 decayProducts->PushProducts(dynamicParticle);
609 if ( fVerboseLevel > 0 ) {
610 G4cout <<
"nofParticles for tracking: " << counter << G4endl;
613 return decayProducts;
623 if ( decayType == fDecayType )
return;
625 fDecayType = decayType;
626 ForceDecay(fDecayType);
void ForceDecayType(EDecayType decayType)
virtual ~G4Pythia6Decayer()
void SetMDME(int i, int j, int m)
Structure for Pythia6 particle properties.
virtual G4DecayProducts * ImportDecayProducts(const G4Track &track)
void SetMDCY(int i, int j, int m)
static Pythia6 * Instance()
void Py1ent(int line, int kf, double pe, double theta, double phi)
int GetMDCY(int i, int j)
std::vector< Pythia6Particle * > ParticleVector
int GetKFDP(int i, int j)
ParticleVector * ImportParticles()
void SetMSTJ(int i, int m)