Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:21:30 

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(new ParticleVector())
0067 {
0068   /// Standard constructor
0069 
0070   
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                                           const 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};
0262   G4int mult[2] = {1, 1};
0263   ForceParticleDecay(iKstar0, products, mult, 2);
0264 
0265   // for Ds -> Phi pi+
0266   G4int iPhi = 333;
0267   ForceParticleDecay(iPhi, iKPlus, 2);  // Phi->K+K-
0268 
0269   G4int decayP1[kNHadrons][3] = {
0270       {iKMinus, iPiPlus, iPiPlus},
0271       {iKMinus, iPiPlus, 0},
0272       {iKPlus, iKstarbar0, 0},
0273       {-1, -1, -1}};
0274   G4int decayP2[kNHadrons][3] = {
0275       {iKstarbar0, iPiPlus, 0},
0276       {-1, -1, -1},
0277       {iPhi, iPiPlus, 0},
0278       {-1, -1, -1}};
0279 
0280   Pythia6* pythia6 = Pythia6::Instance();
0281   for (G4int ihadron = 0; ihadron < kNHadrons; ihadron++)
0282   {
0283     G4int kc = pythia6->Pycomp(hadron[ihadron]);
0284     pythia6->SetMDCY(kc, 1, 1);
0285     G4int ifirst = pythia6->GetMDCY(kc, 2);
0286     G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
0287 
0288     for (channel = ifirst; channel <= ilast; channel++)
0289     {
0290       if ((pythia6->GetKFDP(channel, 1) == decayP1[ihadron][0] &&
0291            pythia6->GetKFDP(channel, 2) == decayP1[ihadron][1] &&
0292            pythia6->GetKFDP(channel, 3) == decayP1[ihadron][2] &&
0293            pythia6->GetKFDP(channel, 4) == 0) ||
0294           (pythia6->GetKFDP(channel, 1) == decayP2[ihadron][0] &&
0295            pythia6->GetKFDP(channel, 2) == decayP2[ihadron][1] &&
0296            pythia6->GetKFDP(channel, 3) == decayP2[ihadron][2] &&
0297            pythia6->GetKFDP(channel, 4) == 0))
0298       {
0299         pythia6->SetMDME(channel, 1, 1);
0300       }
0301       else
0302       {
0303         pythia6->SetMDME(channel, 1, 0);
0304       }  // selected channel ?
0305     }    // decay channels
0306   }      // hadrons
0307 }
0308 
0309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0310 
0311 void G4Pythia6Decayer::ForceOmega()
0312 {
0313   /// Force Omega -> Lambda K- Decay
0314 
0315   Pythia6* pythia6 = Pythia6::Instance();
0316 
0317   G4int iLambda0 = 3122;
0318   G4int iKMinus = -321;
0319 
0320   G4int kc = pythia6->Pycomp(3334);
0321   pythia6->SetMDCY(kc, 1, 1);
0322   G4int ifirst = pythia6->GetMDCY(kc, 2);
0323   G4int ilast = ifirst + pythia6->GetMDCY(kc, 3) - 1;
0324 
0325   for (G4int channel = ifirst; channel <= ilast; channel++)
0326   {
0327     if (pythia6->GetKFDP(channel, 1) == iLambda0 &&
0328         pythia6->GetKFDP(channel, 2) == iKMinus &&
0329         pythia6->GetKFDP(channel, 3) == 0)
0330     {
0331       pythia6->SetMDME(channel, 1, 1);
0332     }
0333     else
0334     {
0335       pythia6->SetMDME(channel, 1, 0);
0336     }
0337     // selected channel ?
0338   }  // decay channels
0339 }
0340 
0341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0342 
0343 void G4Pythia6Decayer::ForceDecay(EDecayType decayType)
0344 {
0345   /// Force a particle decay mode
0346 
0347   Pythia6::Instance()->SetMSTJ(21, 2);
0348 
0349   if (fDecayType == kNoDecayHeavy)
0350   {
0351     return;
0352   }
0353 
0354   //
0355   // select mode
0356   G4int products[3];
0357   G4int mult[3];
0358 
0359   switch (decayType)
0360   {
0361   case kHardMuons:
0362     products[0] = 13;
0363     products[1] = 443;
0364     products[2] = 100443;
0365     mult[0] = 1;
0366     mult[1] = 1;
0367     mult[2] = 1;
0368     ForceParticleDecay(511, products, mult, 3);
0369     ForceParticleDecay(521, products, mult, 3);
0370     ForceParticleDecay(531, products, mult, 3);
0371     ForceParticleDecay(5122, products, mult, 3);
0372     ForceParticleDecay(5132, products, mult, 3);
0373     ForceParticleDecay(5232, products, mult, 3);
0374     ForceParticleDecay(5332, products, mult, 3);
0375     ForceParticleDecay(100443, 443, 1);  // Psi'  -> J/Psi X
0376     ForceParticleDecay(443, 13, 2);      // J/Psi -> mu+ mu-
0377 
0378     ForceParticleDecay(411, 13, 1);   // D+/-
0379     ForceParticleDecay(421, 13, 1);   // D0
0380     ForceParticleDecay(431, 13, 1);   // D_s
0381     ForceParticleDecay(4122, 13, 1);  // Lambda_c
0382     ForceParticleDecay(4132, 13, 1);  // Xsi_c
0383     ForceParticleDecay(4232, 13, 1);  // Sigma_c
0384     ForceParticleDecay(4332, 13, 1);  // Omega_c
0385     break;
0386 
0387   case kSemiMuonic:
0388     ForceParticleDecay(411, 13, 1);   // D+/-
0389     ForceParticleDecay(421, 13, 1);   // D0
0390     ForceParticleDecay(431, 13, 1);   // D_s
0391     ForceParticleDecay(4122, 13, 1);  // Lambda_c
0392     ForceParticleDecay(4132, 13, 1);  // Xsi_c
0393     ForceParticleDecay(4232, 13, 1);  // Sigma_c
0394     ForceParticleDecay(4332, 13, 1);  // Omega_c
0395     ForceParticleDecay(511, 13, 1);   // B0
0396     ForceParticleDecay(521, 13, 1);   // B+/-
0397     ForceParticleDecay(531, 13, 1);   // B_s
0398     ForceParticleDecay(5122, 13, 1);  // Lambda_b
0399     ForceParticleDecay(5132, 13, 1);  // Xsi_b
0400     ForceParticleDecay(5232, 13, 1);  // Sigma_b
0401     ForceParticleDecay(5332, 13, 1);  // Omega_b
0402     break;
0403 
0404   case kDiMuon:
0405     ForceParticleDecay(113, 13, 2);     // rho
0406     ForceParticleDecay(221, 13, 2);     // eta
0407     ForceParticleDecay(223, 13, 2);     // omega
0408     ForceParticleDecay(333, 13, 2);     // phi
0409     ForceParticleDecay(443, 13, 2);     // J/Psi
0410     ForceParticleDecay(100443, 13, 2);  // Psi'
0411     ForceParticleDecay(553, 13, 2);     // Upsilon
0412     ForceParticleDecay(100553, 13, 2);  // Upsilon'
0413     ForceParticleDecay(200553, 13, 2);  // Upsilon''
0414     break;
0415 
0416   case kSemiElectronic:
0417     ForceParticleDecay(411, 11, 1);   // D+/-
0418     ForceParticleDecay(421, 11, 1);   // D0
0419     ForceParticleDecay(431, 11, 1);   // D_s
0420     ForceParticleDecay(4122, 11, 1);  // Lambda_c
0421     ForceParticleDecay(4132, 11, 1);  // Xsi_c
0422     ForceParticleDecay(4232, 11, 1);  // Sigma_c
0423     ForceParticleDecay(4332, 11, 1);  // Omega_c
0424     ForceParticleDecay(511, 11, 1);   // B0
0425     ForceParticleDecay(521, 11, 1);   // B+/-
0426     ForceParticleDecay(531, 11, 1);   // B_s
0427     ForceParticleDecay(5122, 11, 1);  // Lambda_b
0428     ForceParticleDecay(5132, 11, 1);  // Xsi_b
0429     ForceParticleDecay(5232, 11, 1);  // Sigma_b
0430     ForceParticleDecay(5332, 11, 1);  // Omega_b
0431     break;
0432 
0433   case kDiElectron:
0434     ForceParticleDecay(113, 11, 2);     // rho
0435     ForceParticleDecay(333, 11, 2);     // phi
0436     ForceParticleDecay(221, 11, 2);     // eta
0437     ForceParticleDecay(223, 11, 2);     // omega
0438     ForceParticleDecay(443, 11, 2);     // J/Psi
0439     ForceParticleDecay(100443, 11, 2);  // Psi'
0440     ForceParticleDecay(553, 11, 2);     // Upsilon
0441     ForceParticleDecay(100553, 11, 2);  // Upsilon'
0442     ForceParticleDecay(200553, 11, 2);  // Upsilon''
0443     break;
0444 
0445   case kBJpsiDiMuon:
0446 
0447     products[0] = 443;
0448     products[1] = 100443;
0449     mult[0] = 1;
0450     mult[1] = 1;
0451 
0452     ForceParticleDecay(511, products, mult, 2);   // B0   -> J/Psi (Psi') X
0453     ForceParticleDecay(521, products, mult, 2);   // B+/- -> J/Psi (Psi') X
0454     ForceParticleDecay(531, products, mult, 2);   // B_s  -> J/Psi (Psi') X
0455     ForceParticleDecay(5122, products, mult, 2);  // Lambda_b -> J/Psi (Psi')X
0456     ForceParticleDecay(100443, 443, 1);           // Psi'  -> J/Psi X
0457     ForceParticleDecay(443, 13, 2);               // J/Psi -> mu+ mu-
0458     break;
0459 
0460   case kBPsiPrimeDiMuon:
0461     ForceParticleDecay(511, 100443, 1);   // B0
0462     ForceParticleDecay(521, 100443, 1);   // B+/-
0463     ForceParticleDecay(531, 100443, 1);   // B_s
0464     ForceParticleDecay(5122, 100443, 1);  // Lambda_b
0465     ForceParticleDecay(100443, 13, 2);    // Psi'
0466     break;
0467 
0468   case kBJpsiDiElectron:
0469     ForceParticleDecay(511, 443, 1);   // B0
0470     ForceParticleDecay(521, 443, 1);   // B+/-
0471     ForceParticleDecay(531, 443, 1);   // B_s
0472     ForceParticleDecay(5122, 443, 1);  // Lambda_b
0473     ForceParticleDecay(443, 11, 2);    // J/Psi
0474     break;
0475 
0476   case kBJpsi:
0477     ForceParticleDecay(511, 443, 1);   // B0
0478     ForceParticleDecay(521, 443, 1);   // B+/-
0479     ForceParticleDecay(531, 443, 1);   // B_s
0480     ForceParticleDecay(5122, 443, 1);  // Lambda_b
0481     break;
0482 
0483   case kBPsiPrimeDiElectron:
0484     ForceParticleDecay(511, 100443, 1);   // B0
0485     ForceParticleDecay(521, 100443, 1);   // B+/-
0486     ForceParticleDecay(531, 100443, 1);   // B_s
0487     ForceParticleDecay(5122, 100443, 1);  // Lambda_b
0488     ForceParticleDecay(100443, 11, 2);    // Psi'
0489     break;
0490 
0491   case kPiToMu:
0492     ForceParticleDecay(211, 13, 1);  // pi->mu
0493     break;
0494 
0495   case kKaToMu:
0496     ForceParticleDecay(321, 13, 1);  // K->mu
0497     break;
0498 
0499   case kWToMuon:
0500     ForceParticleDecay(24, 13, 1);  // W -> mu
0501     break;
0502 
0503   case kWToCharm:
0504     ForceParticleDecay(24, 4, 1);  // W -> c
0505     break;
0506 
0507   case kWToCharmToMuon:
0508     ForceParticleDecay(24, 4, 1);     // W -> c
0509     ForceParticleDecay(411, 13, 1);   // D+/- -> mu
0510     ForceParticleDecay(421, 13, 1);   // D0  -> mu
0511     ForceParticleDecay(431, 13, 1);   // D_s  -> mu
0512     ForceParticleDecay(4122, 13, 1);  // Lambda_c
0513     ForceParticleDecay(4132, 13, 1);  // Xsi_c
0514     ForceParticleDecay(4232, 13, 1);  // Sigma_c
0515     ForceParticleDecay(4332, 13, 1);  // Omega_c
0516     break;
0517 
0518   case kZDiMuon:
0519     ForceParticleDecay(23, 13, 2);  // Z -> mu+ mu-
0520     break;
0521 
0522   case kHadronicD:
0523     ForceHadronicD();
0524     break;
0525 
0526   case kPhiKK:
0527     ForceParticleDecay(333, 321, 2);  // Phi->K+K-
0528     break;
0529 
0530   case kOmega:
0531     ForceOmega();
0532 
0533   case kAll:
0534     break;
0535 
0536   case kNoDecay:
0537     Pythia6::Instance()->SetMSTJ(21, 0);
0538     break;
0539 
0540   case kNoDecayHeavy:
0541   case kMaxDecay:
0542     break;
0543   }
0544 }
0545 
0546 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0547 
0548 void G4Pythia6Decayer::Decay(G4int pdg, const CLHEP::HepLorentzVector& p)
0549 {
0550   /// Decay a particle of type IDPART (PDG code) and momentum P.
0551 
0552   Pythia6::Instance()->Py1ent(0, pdg, p.e(), p.theta(), p.phi());
0553 }
0554 
0555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0556 
0557 G4int G4Pythia6Decayer::ImportParticles(ParticleVector* particles)
0558 {
0559   /// Get the decay products into the passed PARTICLES vector
0560 
0561   return Pythia6::Instance()->ImportParticles(particles, "All");
0562 }
0563 
0564 //
0565 // public methods
0566 //
0567 
0568 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0569 
0570 G4DecayProducts* G4Pythia6Decayer::ImportDecayProducts(const G4Track& track)
0571 {
0572   /// Import decay products
0573 
0574   // get particle momentum
0575   G4ThreeVector momentum = track.GetMomentum();
0576   G4double etot = track.GetDynamicParticle()->GetTotalEnergy();
0577   ;
0578   CLHEP::HepLorentzVector p;
0579   p[0] = momentum.x() / GeV;
0580   p[1] = momentum.y() / GeV;
0581   p[2] = momentum.z() / GeV;
0582   p[3] = etot / GeV;
0583 
0584   // get particle PDG
0585   // ask G4Pythia6Decayer to get PDG encoding
0586   // (in order to get PDG from extended TDatabasePDG
0587   // in case the standard PDG code is not defined)
0588   G4ParticleDefinition* particleDef = track.GetDefinition();
0589   G4int pdgEncoding = particleDef->GetPDGEncoding();
0590 
0591   // let Pythia6Decayer decay the particle
0592   // and import the decay products
0593   Decay(pdgEncoding, p);
0594   G4int nofParticles = ImportParticles(fDecayProductsArray);
0595 
0596   if (fVerboseLevel > 0)
0597   {
0598     std::cout << "nofParticles: " << nofParticles << std::endl;
0599   }
0600 
0601   // convert decay products Pythia6Particle type
0602   // to G4DecayProducts
0603   G4DecayProducts* decayProducts = new G4DecayProducts(*(track.GetDynamicParticle()));
0604 
0605   G4int counter = 0;
0606   for (G4int i = 0; i < nofParticles; i++)
0607   {
0608     // get particle from ParticleVector
0609     Pythia6Particle* particle = (*fDecayProductsArray)[i];
0610 
0611     G4int status = particle->fKS;
0612     G4int pdg = particle->fKF;
0613     if (status > 0 && status < 11 &&
0614         std::abs(pdg) != 12 && std::abs(pdg) != 14 && std::abs(pdg) != 16)
0615     {
0616       // pass to tracking final particles only;
0617       // skip neutrinos
0618 
0619       if (fVerboseLevel > 0)
0620       {
0621         std::cout << "  " << i << "th particle PDG: " << pdg << "   ";
0622       }
0623 
0624       // create G4DynamicParticle
0625       G4DynamicParticle* dynamicParticle = CreateDynamicParticle(particle);
0626 
0627       if (dynamicParticle)
0628       {
0629         if (fVerboseLevel > 0)
0630         {
0631           std::cout << "  G4 particle name: "
0632                     << dynamicParticle->GetDefinition()->GetParticleName()
0633                     << std::endl;
0634         }
0635 
0636         // add dynamicParticle to decayProducts
0637         decayProducts->PushProducts(dynamicParticle);
0638 
0639         counter++;
0640       }
0641     }
0642   }
0643   if (fVerboseLevel > 0)
0644   {
0645     std::cout << "nofParticles for tracking: " << counter << std::endl;
0646   }
0647 
0648   return decayProducts;
0649 }
0650 
0651 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0652 
0653 void G4Pythia6Decayer::ForceDecayType(EDecayType decayType)
0654 {
0655   /// Force a given decay type
0656 
0657   // Do nothing if the decay type is not different from current one
0658   if (decayType == fDecayType)
0659   {
0660     return;
0661   }
0662 
0663   fDecayType = decayType;
0664   ForceDecay(fDecayType);
0665 }
[ 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!