Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-08-22 08:26:25

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file G4Pythia6Decayer.cc
0027 /// \brief Implementation of the G4Pythia6Decayer class
0028 
0029 // ----------------------------------------------------------------------------
0030 // According to TPythia6Decayer class in Root:
0031 // http://root.cern/
0032 // see http://root.cern/root/License.html
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0052 
0053 G4Pythia6Decayer::G4Pythia6Decayer()
0054   : G4VExtDecayer("G4Pythia6Decayer"),
0055     fMessenger(this),
0056     fVerboseLevel(0),
0057     fDecayType(fgkDefaultDecayType),
0058     fDecayProductsArray(0)
0059 {
0060   /// Standard constructor
0061 
0062   fDecayProductsArray = new ParticleVector();
0063 
0064   ForceDecay(fDecayType);
0065 }
0066 
0067 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0068 
0069 G4Pythia6Decayer::~G4Pythia6Decayer()
0070 {
0071   /// Destructor
0072 
0073   delete fDecayProductsArray;
0074 }
0075 
0076 //
0077 // private methods
0078 //
0079 
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0081 
0082 G4ParticleDefinition* G4Pythia6Decayer::GetParticleDefinition(const Pythia6Particle* particle,
0083                                                               G4bool warn) const
0084 {
0085   /// Return G4 particle definition for given TParticle
0086 
0087   // get particle definition from G4ParticleTable
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0103 
0104 G4DynamicParticle* G4Pythia6Decayer::CreateDynamicParticle(const Pythia6Particle* particle) const
0105 {
0106   /// Create G4DynamicParticle.
0107 
0108   // get particle properties
0109   const G4ParticleDefinition* particleDefinition = GetParticleDefinition(particle);
0110   if (!particleDefinition) return 0;
0111 
0112   G4ThreeVector momentum = GetParticleMomentum(particle);
0113 
0114   // create G4DynamicParticle
0115   G4DynamicParticle* dynamicParticle = new G4DynamicParticle(particleDefinition, momentum);
0116 
0117   return dynamicParticle;
0118 }
0119 
0120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0121 
0122 G4ThreeVector G4Pythia6Decayer::GetParticlePosition(const Pythia6Particle* particle) const
0123 {
0124   /// Return particle vertex position.
0125 
0126   G4ThreeVector position =
0127     G4ThreeVector(particle->fVx * cm, particle->fVy * cm, particle->fVz * cm);
0128   return position;
0129 }
0130 
0131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0132 
0133 G4ThreeVector G4Pythia6Decayer::GetParticleMomentum(const Pythia6Particle* particle) const
0134 {
0135   /// Return particle momentum.
0136 
0137   G4ThreeVector momentum =
0138     G4ThreeVector(particle->fPx * GeV, particle->fPy * GeV, particle->fPz * GeV);
0139   return momentum;
0140 }
0141 
0142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0143 
0144 G4int G4Pythia6Decayer::CountProducts(G4int channel, G4int particle)
0145 {
0146   /// Count number of decay products
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0155 
0156 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int product, G4int mult)
0157 {
0158   /// Force decay of particle into products with multiplicity mult
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   //  Loop over decay channels
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0181 
0182 void G4Pythia6Decayer::ForceParticleDecay(G4int particle, G4int* products, G4int* mult, G4int npart)
0183 {
0184   /// Force decay of particle into products with multiplicity mult
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   //  Loop over decay channels
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0207 
0208 void G4Pythia6Decayer::ForceHadronicD()
0209 {
0210   /// Force golden D decay modes
0211 
0212   const G4int kNHadrons = 4;
0213   G4int channel;
0214   G4int hadron[kNHadrons] = {411, 421, 431, 4112};
0215 
0216   // for D+ -> K0* (-> K- pi+) pi+
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   // for Ds -> Phi pi+
0228   G4int iPhi = 333;
0229   ForceParticleDecay(iPhi, iKPlus, 2);  // Phi->K+K-
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       }  // selected channel ?
0258     }  // decay channels
0259   }  // hadrons
0260 }
0261 
0262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0263 
0264 void G4Pythia6Decayer::ForceOmega()
0265 {
0266   /// Force Omega -> Lambda K- Decay
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     // selected channel ?
0285   }  // decay channels
0286 }
0287 
0288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0289 
0290 void G4Pythia6Decayer::ForceDecay(EDecayType decayType)
0291 {
0292   /// Force a particle decay mode
0293 
0294   Pythia6::Instance()->SetMSTJ(21, 2);
0295 
0296   if (fDecayType == kNoDecayHeavy) return;
0297 
0298   //
0299   // select mode
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);  // Psi'  -> J/Psi X
0319       ForceParticleDecay(443, 13, 2);  // J/Psi -> mu+ mu-
0320 
0321       ForceParticleDecay(411, 13, 1);  // D+/-
0322       ForceParticleDecay(421, 13, 1);  // D0
0323       ForceParticleDecay(431, 13, 1);  // D_s
0324       ForceParticleDecay(4122, 13, 1);  // Lambda_c
0325       ForceParticleDecay(4132, 13, 1);  // Xsi_c
0326       ForceParticleDecay(4232, 13, 1);  // Sigma_c
0327       ForceParticleDecay(4332, 13, 1);  // Omega_c
0328       break;
0329 
0330     case kSemiMuonic:
0331       ForceParticleDecay(411, 13, 1);  // D+/-
0332       ForceParticleDecay(421, 13, 1);  // D0
0333       ForceParticleDecay(431, 13, 1);  // D_s
0334       ForceParticleDecay(4122, 13, 1);  // Lambda_c
0335       ForceParticleDecay(4132, 13, 1);  // Xsi_c
0336       ForceParticleDecay(4232, 13, 1);  // Sigma_c
0337       ForceParticleDecay(4332, 13, 1);  // Omega_c
0338       ForceParticleDecay(511, 13, 1);  // B0
0339       ForceParticleDecay(521, 13, 1);  // B+/-
0340       ForceParticleDecay(531, 13, 1);  // B_s
0341       ForceParticleDecay(5122, 13, 1);  // Lambda_b
0342       ForceParticleDecay(5132, 13, 1);  // Xsi_b
0343       ForceParticleDecay(5232, 13, 1);  // Sigma_b
0344       ForceParticleDecay(5332, 13, 1);  // Omega_b
0345       break;
0346 
0347     case kDiMuon:
0348       ForceParticleDecay(113, 13, 2);  // rho
0349       ForceParticleDecay(221, 13, 2);  // eta
0350       ForceParticleDecay(223, 13, 2);  // omega
0351       ForceParticleDecay(333, 13, 2);  // phi
0352       ForceParticleDecay(443, 13, 2);  // J/Psi
0353       ForceParticleDecay(100443, 13, 2);  // Psi'
0354       ForceParticleDecay(553, 13, 2);  // Upsilon
0355       ForceParticleDecay(100553, 13, 2);  // Upsilon'
0356       ForceParticleDecay(200553, 13, 2);  // Upsilon''
0357       break;
0358 
0359     case kSemiElectronic:
0360       ForceParticleDecay(411, 11, 1);  // D+/-
0361       ForceParticleDecay(421, 11, 1);  // D0
0362       ForceParticleDecay(431, 11, 1);  // D_s
0363       ForceParticleDecay(4122, 11, 1);  // Lambda_c
0364       ForceParticleDecay(4132, 11, 1);  // Xsi_c
0365       ForceParticleDecay(4232, 11, 1);  // Sigma_c
0366       ForceParticleDecay(4332, 11, 1);  // Omega_c
0367       ForceParticleDecay(511, 11, 1);  // B0
0368       ForceParticleDecay(521, 11, 1);  // B+/-
0369       ForceParticleDecay(531, 11, 1);  // B_s
0370       ForceParticleDecay(5122, 11, 1);  // Lambda_b
0371       ForceParticleDecay(5132, 11, 1);  // Xsi_b
0372       ForceParticleDecay(5232, 11, 1);  // Sigma_b
0373       ForceParticleDecay(5332, 11, 1);  // Omega_b
0374       break;
0375 
0376     case kDiElectron:
0377       ForceParticleDecay(113, 11, 2);  // rho
0378       ForceParticleDecay(333, 11, 2);  // phi
0379       ForceParticleDecay(221, 11, 2);  // eta
0380       ForceParticleDecay(223, 11, 2);  // omega
0381       ForceParticleDecay(443, 11, 2);  // J/Psi
0382       ForceParticleDecay(100443, 11, 2);  // Psi'
0383       ForceParticleDecay(553, 11, 2);  // Upsilon
0384       ForceParticleDecay(100553, 11, 2);  // Upsilon'
0385       ForceParticleDecay(200553, 11, 2);  // Upsilon''
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);  // B0   -> J/Psi (Psi') X
0396       ForceParticleDecay(521, products, mult, 2);  // B+/- -> J/Psi (Psi') X
0397       ForceParticleDecay(531, products, mult, 2);  // B_s  -> J/Psi (Psi') X
0398       ForceParticleDecay(5122, products, mult, 2);  // Lambda_b -> J/Psi (Psi')X
0399       ForceParticleDecay(100443, 443, 1);  // Psi'  -> J/Psi X
0400       ForceParticleDecay(443, 13, 2);  // J/Psi -> mu+ mu-
0401       break;
0402 
0403     case kBPsiPrimeDiMuon:
0404       ForceParticleDecay(511, 100443, 1);  // B0
0405       ForceParticleDecay(521, 100443, 1);  // B+/-
0406       ForceParticleDecay(531, 100443, 1);  // B_s
0407       ForceParticleDecay(5122, 100443, 1);  // Lambda_b
0408       ForceParticleDecay(100443, 13, 2);  // Psi'
0409       break;
0410 
0411     case kBJpsiDiElectron:
0412       ForceParticleDecay(511, 443, 1);  // B0
0413       ForceParticleDecay(521, 443, 1);  // B+/-
0414       ForceParticleDecay(531, 443, 1);  // B_s
0415       ForceParticleDecay(5122, 443, 1);  // Lambda_b
0416       ForceParticleDecay(443, 11, 2);  // J/Psi
0417       break;
0418 
0419     case kBJpsi:
0420       ForceParticleDecay(511, 443, 1);  // B0
0421       ForceParticleDecay(521, 443, 1);  // B+/-
0422       ForceParticleDecay(531, 443, 1);  // B_s
0423       ForceParticleDecay(5122, 443, 1);  // Lambda_b
0424       break;
0425 
0426     case kBPsiPrimeDiElectron:
0427       ForceParticleDecay(511, 100443, 1);  // B0
0428       ForceParticleDecay(521, 100443, 1);  // B+/-
0429       ForceParticleDecay(531, 100443, 1);  // B_s
0430       ForceParticleDecay(5122, 100443, 1);  // Lambda_b
0431       ForceParticleDecay(100443, 11, 2);  // Psi'
0432       break;
0433 
0434     case kPiToMu:
0435       ForceParticleDecay(211, 13, 1);  // pi->mu
0436       break;
0437 
0438     case kKaToMu:
0439       ForceParticleDecay(321, 13, 1);  // K->mu
0440       break;
0441 
0442     case kWToMuon:
0443       ForceParticleDecay(24, 13, 1);  // W -> mu
0444       break;
0445 
0446     case kWToCharm:
0447       ForceParticleDecay(24, 4, 1);  // W -> c
0448       break;
0449 
0450     case kWToCharmToMuon:
0451       ForceParticleDecay(24, 4, 1);  // W -> c
0452       ForceParticleDecay(411, 13, 1);  // D+/- -> mu
0453       ForceParticleDecay(421, 13, 1);  // D0  -> mu
0454       ForceParticleDecay(431, 13, 1);  // D_s  -> mu
0455       ForceParticleDecay(4122, 13, 1);  // Lambda_c
0456       ForceParticleDecay(4132, 13, 1);  // Xsi_c
0457       ForceParticleDecay(4232, 13, 1);  // Sigma_c
0458       ForceParticleDecay(4332, 13, 1);  // Omega_c
0459       break;
0460 
0461     case kZDiMuon:
0462       ForceParticleDecay(23, 13, 2);  // Z -> mu+ mu-
0463       break;
0464 
0465     case kHadronicD:
0466       ForceHadronicD();
0467       break;
0468 
0469     case kPhiKK:
0470       ForceParticleDecay(333, 321, 2);  // Phi->K+K-
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0492 
0493 void G4Pythia6Decayer::Decay(G4int pdg, const CLHEP::HepLorentzVector& p)
0494 {
0495   /// Decay a particle of type IDPART (PDG code) and momentum P.
0496 
0497   Pythia6::Instance()->Py1ent(0, pdg, p.e(), p.theta(), p.phi());
0498 }
0499 
0500 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0501 
0502 G4int G4Pythia6Decayer::ImportParticles(ParticleVector* particles)
0503 {
0504   /// Get the decay products into the passed PARTICLES vector
0505 
0506   return Pythia6::Instance()->ImportParticles(particles, "All");
0507 }
0508 
0509 //
0510 // public methods
0511 //
0512 
0513 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0514 
0515 G4DecayProducts* G4Pythia6Decayer::ImportDecayProducts(const G4Track& track)
0516 {
0517   /// Import decay products
0518 
0519   // get particle momentum
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   // get particle PDG
0530   // ask G4Pythia6Decayer to get PDG encoding
0531   // (in order to get PDG from extended TDatabasePDG
0532   // in case the standard PDG code is not defined)
0533   G4ParticleDefinition* particleDef = track.GetDefinition();
0534   G4int pdgEncoding = particleDef->GetPDGEncoding();
0535 
0536   // let Pythia6Decayer decay the particle
0537   // and import the decay products
0538   Decay(pdgEncoding, p);
0539   G4int nofParticles = ImportParticles(fDecayProductsArray);
0540 
0541   if (fVerboseLevel > 0) {
0542     G4cout << "nofParticles: " << nofParticles << G4endl;
0543   }
0544 
0545   // convert decay products Pythia6Particle type
0546   // to G4DecayProducts
0547   G4DecayProducts* decayProducts = new G4DecayProducts(*(track.GetDynamicParticle()));
0548 
0549   G4int counter = 0;
0550   for (G4int i = 0; i < nofParticles; i++) {
0551     // get particle from ParticleVector
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       // pass to tracking final particles only;
0560       // skip neutrinos
0561 
0562       if (fVerboseLevel > 0) {
0563         G4cout << "  " << i << "th particle PDG: " << pdg << "   ";
0564       }
0565 
0566       // create G4DynamicParticle
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         // add dynamicParticle to decayProducts
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0590 
0591 void G4Pythia6Decayer::ForceDecayType(EDecayType decayType)
0592 {
0593   /// Force a given decay type
0594 
0595   // Do nothing if the decay type is not different from current one
0596   if (decayType == fDecayType) return;
0597 
0598   fDecayType = decayType;
0599   ForceDecay(fDecayType);
0600 }