Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:18:52 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!