decayer6/src/G4Pythia6Decayer.cc

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