File indexing completed on 2026-08-22 08:26:25
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035 #include "G4Pythia6Decayer.hh"
0036
0037 #include "Pythia6.hh"
0038
0039 #include "G4DecayProducts.hh"
0040 #include "G4DecayTable.hh"
0041 #include "G4DynamicParticle.hh"
0042 #include "G4ParticleTable.hh"
0043 #include "G4SystemOfUnits.hh"
0044 #include "G4Track.hh"
0045
0046 #include <CLHEP/Vector/LorentzVector.h>
0047 #include <cmath>
0048
0049 const EDecayType G4Pythia6Decayer::fgkDefaultDecayType = kAll;
0050
0051
0052
0053 G4Pythia6Decayer::G4Pythia6Decayer()
0054 : G4VExtDecayer("G4Pythia6Decayer"),
0055 fMessenger(this),
0056 fVerboseLevel(0),
0057 fDecayType(fgkDefaultDecayType),
0058 fDecayProductsArray(0)
0059 {
0060
0061
0062 fDecayProductsArray = new ParticleVector();
0063
0064 ForceDecay(fDecayType);
0065 }
0066
0067
0068
0069 G4Pythia6Decayer::~G4Pythia6Decayer()
0070 {
0071
0072
0073 delete fDecayProductsArray;
0074 }
0075
0076
0077
0078
0079
0080
0081
0082 G4ParticleDefinition* G4Pythia6Decayer::GetParticleDefinition(const Pythia6Particle* particle,
0083 G4bool warn) const
0084 {
0085
0086
0087
0088 G4int pdgEncoding = particle->fKF;
0089 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0090 G4ParticleDefinition* particleDefinition = 0;
0091 if (pdgEncoding != 0) particleDefinition = particleTable->FindParticle(pdgEncoding);
0092
0093 if (particleDefinition == 0 && warn) {
0094 G4cerr << "G4Pythia6Decayer: GetParticleDefinition: " << std::endl
0095 << "G4ParticleTable::FindParticle() for particle with PDG = " << pdgEncoding
0096 << " failed." << std::endl;
0097 }
0098
0099 return particleDefinition;
0100 }
0101
0102
0103
0104 G4DynamicParticle* G4Pythia6Decayer::CreateDynamicParticle(const Pythia6Particle* particle) const
0105 {
0106
0107
0108
0109 const G4ParticleDefinition* particleDefinition = GetParticleDefinition(particle);
0110 if (!particleDefinition) return 0;
0111
0112 G4ThreeVector momentum = GetParticleMomentum(particle);
0113
0114
0115 G4DynamicParticle* dynamicParticle = new G4DynamicParticle(particleDefinition, momentum);
0116
0117 return dynamicParticle;
0118 }
0119
0120
0121
0122 G4ThreeVector G4Pythia6Decayer::GetParticlePosition(const Pythia6Particle* particle) const
0123 {
0124
0125
0126 G4ThreeVector position =
0127 G4ThreeVector(particle->fVx * cm, particle->fVy * cm, particle->fVz * cm);
0128 return position;
0129 }
0130
0131
0132
0133 G4ThreeVector G4Pythia6Decayer::GetParticleMomentum(const Pythia6Particle* particle) const
0134 {
0135
0136
0137 G4ThreeVector momentum =
0138 G4ThreeVector(particle->fPx * GeV, particle->fPy * GeV, particle->fPz * GeV);
0139 return momentum;
0140 }
0141
0142
0143
0144 G4int G4Pythia6Decayer::CountProducts(G4int channel, G4int particle)
0145 {
0146
0147
0148 G4int np = 0;
0149 for (G4int i = 1; i <= 5; i++)
0150 if (std::abs(Pythia6::Instance()->GetKFDP(channel, i)) == particle) np++;
0151 return np;
0152 }
0153
0154
0155
0156 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int product, G4int mult)
0157 {
0158
0159
0160 Pythia6* pythia6 = Pythia6::Instance();
0161
0162 G4int kc = pythia6->Pycomp(particle);
0163 pythia6->SetMDCY(kc, 1, 1);
0164
0165 G4int ifirst = pythia6->GetMDCY(kc, 2);
0166 G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
0167
0168
0169
0170 for (G4int channel = ifirst; channel <= ilast; channel++) {
0171 if (CountProducts(channel, product) >= mult) {
0172 pythia6->SetMDME(channel, 1, 1);
0173 }
0174 else {
0175 pythia6->SetMDME(channel, 1, 0);
0176 }
0177 }
0178 }
0179
0180
0181
0182 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int* products, G4int* mult, G4int npart)
0183 {
0184
0185
0186 Pythia6* pythia6 = Pythia6::Instance();
0187
0188 G4int kc = pythia6->Pycomp(particle);
0189 pythia6->SetMDCY(kc, 1, 1);
0190 G4int ifirst = pythia6->GetMDCY(kc, 2);
0191 G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
0192
0193
0194 for (G4int channel = ifirst; channel <= ilast; channel++) {
0195 G4int nprod = 0;
0196 for (G4int i = 0; i < npart; i++)
0197 nprod += (CountProducts(channel, products[i]) >= mult[i]);
0198 if (nprod)
0199 pythia6->SetMDME(channel, 1, 1);
0200 else {
0201 pythia6->SetMDME(channel, 1, 0);
0202 }
0203 }
0204 }
0205
0206
0207
0208 void G4Pythia6Decayer::ForceHadronicD()
0209 {
0210
0211
0212 const G4int kNHadrons = 4;
0213 G4int channel;
0214 G4int hadron[kNHadrons] = {411, 421, 431, 4112};
0215
0216
0217 G4int iKstar0 = 313;
0218 G4int iKstarbar0 = -313;
0219 G4int iKPlus = 321;
0220 G4int iKMinus = -321;
0221 G4int iPiPlus = 211;
0222 G4int iPiMinus = -211;
0223
0224 G4int products[2] = {iKPlus, iPiMinus}, mult[2] = {1, 1};
0225 ForceParticleDecay(iKstar0, products, mult, 2);
0226
0227
0228 G4int iPhi = 333;
0229 ForceParticleDecay(iPhi, iKPlus, 2);
0230
0231 G4int decayP1[kNHadrons][3] = {
0232 {iKMinus, iPiPlus, iPiPlus}, {iKMinus, iPiPlus, 0}, {iKPlus, iKstarbar0, 0}, {-1, -1, -1}};
0233 G4int decayP2[kNHadrons][3] = {
0234 {iKstarbar0, iPiPlus, 0}, {-1, -1, -1}, {iPhi, iPiPlus, 0}, {-1, -1, -1}};
0235
0236 Pythia6* pythia6 = Pythia6::Instance();
0237 for (G4int ihadron = 0; ihadron < kNHadrons; ihadron++) {
0238 G4int kc = pythia6->Pycomp(hadron[ihadron]);
0239 pythia6->SetMDCY(kc, 1, 1);
0240 G4int ifirst = pythia6->GetMDCY(kc, 2);
0241 G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
0242
0243 for (channel = ifirst; channel <= ilast; channel++) {
0244 if ((pythia6->GetKFDP(channel, 1) == decayP1[ihadron][0]
0245 && pythia6->GetKFDP(channel, 2) == decayP1[ihadron][1]
0246 && pythia6->GetKFDP(channel, 3) == decayP1[ihadron][2]
0247 && pythia6->GetKFDP(channel, 4) == 0)
0248 || (pythia6->GetKFDP(channel, 1) == decayP2[ihadron][0]
0249 && pythia6->GetKFDP(channel, 2) == decayP2[ihadron][1]
0250 && pythia6->GetKFDP(channel, 3) == decayP2[ihadron][2]
0251 && pythia6->GetKFDP(channel, 4) == 0))
0252 {
0253 pythia6->SetMDME(channel, 1, 1);
0254 }
0255 else {
0256 pythia6->SetMDME(channel, 1, 0);
0257 }
0258 }
0259 }
0260 }
0261
0262
0263
0264 void G4Pythia6Decayer::ForceOmega()
0265 {
0266
0267
0268 Pythia6* pythia6 = Pythia6::Instance();
0269
0270 G4int iLambda0 = 3122;
0271 G4int iKMinus = -321;
0272
0273 G4int kc = pythia6->Pycomp(3334);
0274 pythia6->SetMDCY(kc, 1, 1);
0275 G4int ifirst = pythia6->GetMDCY(kc, 2);
0276 G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
0277
0278 for (G4int channel = ifirst; channel <= ilast; channel++) {
0279 if (pythia6->GetKFDP(channel, 1) == iLambda0 && pythia6->GetKFDP(channel, 2) == iKMinus
0280 && pythia6->GetKFDP(channel, 3) == 0)
0281 pythia6->SetMDME(channel, 1, 1);
0282 else
0283 pythia6->SetMDME(channel, 1, 0);
0284
0285 }
0286 }
0287
0288
0289
0290 void G4Pythia6Decayer::ForceDecay(EDecayType decayType)
0291 {
0292
0293
0294 Pythia6::Instance()->SetMSTJ(21, 2);
0295
0296 if (fDecayType == kNoDecayHeavy) return;
0297
0298
0299
0300 G4int products[3];
0301 G4int mult[3];
0302
0303 switch (decayType) {
0304 case kHardMuons:
0305 products[0] = 13;
0306 products[1] = 443;
0307 products[2] = 100443;
0308 mult[0] = 1;
0309 mult[1] = 1;
0310 mult[2] = 1;
0311 ForceParticleDecay(511, products, mult, 3);
0312 ForceParticleDecay(521, products, mult, 3);
0313 ForceParticleDecay(531, products, mult, 3);
0314 ForceParticleDecay(5122, products, mult, 3);
0315 ForceParticleDecay(5132, products, mult, 3);
0316 ForceParticleDecay(5232, products, mult, 3);
0317 ForceParticleDecay(5332, products, mult, 3);
0318 ForceParticleDecay(100443, 443, 1);
0319 ForceParticleDecay(443, 13, 2);
0320
0321 ForceParticleDecay(411, 13, 1);
0322 ForceParticleDecay(421, 13, 1);
0323 ForceParticleDecay(431, 13, 1);
0324 ForceParticleDecay(4122, 13, 1);
0325 ForceParticleDecay(4132, 13, 1);
0326 ForceParticleDecay(4232, 13, 1);
0327 ForceParticleDecay(4332, 13, 1);
0328 break;
0329
0330 case kSemiMuonic:
0331 ForceParticleDecay(411, 13, 1);
0332 ForceParticleDecay(421, 13, 1);
0333 ForceParticleDecay(431, 13, 1);
0334 ForceParticleDecay(4122, 13, 1);
0335 ForceParticleDecay(4132, 13, 1);
0336 ForceParticleDecay(4232, 13, 1);
0337 ForceParticleDecay(4332, 13, 1);
0338 ForceParticleDecay(511, 13, 1);
0339 ForceParticleDecay(521, 13, 1);
0340 ForceParticleDecay(531, 13, 1);
0341 ForceParticleDecay(5122, 13, 1);
0342 ForceParticleDecay(5132, 13, 1);
0343 ForceParticleDecay(5232, 13, 1);
0344 ForceParticleDecay(5332, 13, 1);
0345 break;
0346
0347 case kDiMuon:
0348 ForceParticleDecay(113, 13, 2);
0349 ForceParticleDecay(221, 13, 2);
0350 ForceParticleDecay(223, 13, 2);
0351 ForceParticleDecay(333, 13, 2);
0352 ForceParticleDecay(443, 13, 2);
0353 ForceParticleDecay(100443, 13, 2);
0354 ForceParticleDecay(553, 13, 2);
0355 ForceParticleDecay(100553, 13, 2);
0356 ForceParticleDecay(200553, 13, 2);
0357 break;
0358
0359 case kSemiElectronic:
0360 ForceParticleDecay(411, 11, 1);
0361 ForceParticleDecay(421, 11, 1);
0362 ForceParticleDecay(431, 11, 1);
0363 ForceParticleDecay(4122, 11, 1);
0364 ForceParticleDecay(4132, 11, 1);
0365 ForceParticleDecay(4232, 11, 1);
0366 ForceParticleDecay(4332, 11, 1);
0367 ForceParticleDecay(511, 11, 1);
0368 ForceParticleDecay(521, 11, 1);
0369 ForceParticleDecay(531, 11, 1);
0370 ForceParticleDecay(5122, 11, 1);
0371 ForceParticleDecay(5132, 11, 1);
0372 ForceParticleDecay(5232, 11, 1);
0373 ForceParticleDecay(5332, 11, 1);
0374 break;
0375
0376 case kDiElectron:
0377 ForceParticleDecay(113, 11, 2);
0378 ForceParticleDecay(333, 11, 2);
0379 ForceParticleDecay(221, 11, 2);
0380 ForceParticleDecay(223, 11, 2);
0381 ForceParticleDecay(443, 11, 2);
0382 ForceParticleDecay(100443, 11, 2);
0383 ForceParticleDecay(553, 11, 2);
0384 ForceParticleDecay(100553, 11, 2);
0385 ForceParticleDecay(200553, 11, 2);
0386 break;
0387
0388 case kBJpsiDiMuon:
0389
0390 products[0] = 443;
0391 products[1] = 100443;
0392 mult[0] = 1;
0393 mult[1] = 1;
0394
0395 ForceParticleDecay(511, products, mult, 2);
0396 ForceParticleDecay(521, products, mult, 2);
0397 ForceParticleDecay(531, products, mult, 2);
0398 ForceParticleDecay(5122, products, mult, 2);
0399 ForceParticleDecay(100443, 443, 1);
0400 ForceParticleDecay(443, 13, 2);
0401 break;
0402
0403 case kBPsiPrimeDiMuon:
0404 ForceParticleDecay(511, 100443, 1);
0405 ForceParticleDecay(521, 100443, 1);
0406 ForceParticleDecay(531, 100443, 1);
0407 ForceParticleDecay(5122, 100443, 1);
0408 ForceParticleDecay(100443, 13, 2);
0409 break;
0410
0411 case kBJpsiDiElectron:
0412 ForceParticleDecay(511, 443, 1);
0413 ForceParticleDecay(521, 443, 1);
0414 ForceParticleDecay(531, 443, 1);
0415 ForceParticleDecay(5122, 443, 1);
0416 ForceParticleDecay(443, 11, 2);
0417 break;
0418
0419 case kBJpsi:
0420 ForceParticleDecay(511, 443, 1);
0421 ForceParticleDecay(521, 443, 1);
0422 ForceParticleDecay(531, 443, 1);
0423 ForceParticleDecay(5122, 443, 1);
0424 break;
0425
0426 case kBPsiPrimeDiElectron:
0427 ForceParticleDecay(511, 100443, 1);
0428 ForceParticleDecay(521, 100443, 1);
0429 ForceParticleDecay(531, 100443, 1);
0430 ForceParticleDecay(5122, 100443, 1);
0431 ForceParticleDecay(100443, 11, 2);
0432 break;
0433
0434 case kPiToMu:
0435 ForceParticleDecay(211, 13, 1);
0436 break;
0437
0438 case kKaToMu:
0439 ForceParticleDecay(321, 13, 1);
0440 break;
0441
0442 case kWToMuon:
0443 ForceParticleDecay(24, 13, 1);
0444 break;
0445
0446 case kWToCharm:
0447 ForceParticleDecay(24, 4, 1);
0448 break;
0449
0450 case kWToCharmToMuon:
0451 ForceParticleDecay(24, 4, 1);
0452 ForceParticleDecay(411, 13, 1);
0453 ForceParticleDecay(421, 13, 1);
0454 ForceParticleDecay(431, 13, 1);
0455 ForceParticleDecay(4122, 13, 1);
0456 ForceParticleDecay(4132, 13, 1);
0457 ForceParticleDecay(4232, 13, 1);
0458 ForceParticleDecay(4332, 13, 1);
0459 break;
0460
0461 case kZDiMuon:
0462 ForceParticleDecay(23, 13, 2);
0463 break;
0464
0465 case kHadronicD:
0466 ForceHadronicD();
0467 break;
0468
0469 case kPhiKK:
0470 ForceParticleDecay(333, 321, 2);
0471 break;
0472
0473 case kOmega:
0474 ForceOmega();
0475
0476 case kAll:
0477 break;
0478
0479 case kNoDecay:
0480 Pythia6::Instance()->SetMSTJ(21, 0);
0481 break;
0482
0483 case kNoDecayHeavy:
0484 break;
0485
0486 case kMaxDecay:
0487 break;
0488 }
0489 }
0490
0491
0492
0493 void G4Pythia6Decayer::Decay(G4int pdg, const CLHEP::HepLorentzVector& p)
0494 {
0495
0496
0497 Pythia6::Instance()->Py1ent(0, pdg, p.e(), p.theta(), p.phi());
0498 }
0499
0500
0501
0502 G4int G4Pythia6Decayer::ImportParticles(ParticleVector* particles)
0503 {
0504
0505
0506 return Pythia6::Instance()->ImportParticles(particles, "All");
0507 }
0508
0509
0510
0511
0512
0513
0514
0515 G4DecayProducts* G4Pythia6Decayer::ImportDecayProducts(const G4Track& track)
0516 {
0517
0518
0519
0520 G4ThreeVector momentum = track.GetMomentum();
0521 G4double etot = track.GetDynamicParticle()->GetTotalEnergy();
0522 ;
0523 CLHEP::HepLorentzVector p;
0524 p[0] = momentum.x() / GeV;
0525 p[1] = momentum.y() / GeV;
0526 p[2] = momentum.z() / GeV;
0527 p[3] = etot / GeV;
0528
0529
0530
0531
0532
0533 G4ParticleDefinition* particleDef = track.GetDefinition();
0534 G4int pdgEncoding = particleDef->GetPDGEncoding();
0535
0536
0537
0538 Decay(pdgEncoding, p);
0539 G4int nofParticles = ImportParticles(fDecayProductsArray);
0540
0541 if (fVerboseLevel > 0) {
0542 G4cout << "nofParticles: " << nofParticles << G4endl;
0543 }
0544
0545
0546
0547 G4DecayProducts* decayProducts = new G4DecayProducts(*(track.GetDynamicParticle()));
0548
0549 G4int counter = 0;
0550 for (G4int i = 0; i < nofParticles; i++) {
0551
0552 Pythia6Particle* particle = (*fDecayProductsArray)[i];
0553
0554 G4int status = particle->fKS;
0555 G4int pdg = particle->fKF;
0556 if (status > 0 && status < 11 && std::abs(pdg) != 12 && std::abs(pdg) != 14
0557 && std::abs(pdg) != 16)
0558 {
0559
0560
0561
0562 if (fVerboseLevel > 0) {
0563 G4cout << " " << i << "th particle PDG: " << pdg << " ";
0564 }
0565
0566
0567 G4DynamicParticle* dynamicParticle = CreateDynamicParticle(particle);
0568
0569 if (dynamicParticle) {
0570 if (fVerboseLevel > 0) {
0571 G4cout << " G4 particle name: " << dynamicParticle->GetDefinition()->GetParticleName()
0572 << G4endl;
0573 }
0574
0575
0576 decayProducts->PushProducts(dynamicParticle);
0577
0578 counter++;
0579 }
0580 }
0581 }
0582 if (fVerboseLevel > 0) {
0583 G4cout << "nofParticles for tracking: " << counter << G4endl;
0584 }
0585
0586 return decayProducts;
0587 }
0588
0589
0590
0591 void G4Pythia6Decayer::ForceDecayType(EDecayType decayType)
0592 {
0593
0594
0595
0596 if (decayType == fDecayType) return;
0597
0598 fDecayType = decayType;
0599 ForceDecay(fDecayType);
0600 }