DDG4/plugins/Geant4ExtraParticles.cpp

0001 //
0002 // Authors: Tomohiko Tanabe <tomohiko@icepp.s.u-tokyo.ac.jp>
0003 //          Taikan Suehara <suehara@icepp.s.u-tokyo.ac.jp>
0004 // Ported from Mokka by A.Sailer (CERN )
0005 //
0006
0007 /** \addtogroup Geant4PhysicsConstructor
0008  *
0009  * @{
0010  * \package Geant4ExtraParticles
0011  *
0012  * \brief PhysicsConstructor to add additional particles to geant
0013
0014  When enabled this constructor will read a particle.tbl file and add all particles in it to the geant4 particle table.
0015  This can be used to let geant know about some B-baryons with non-zero lifetime for example.
0016  *
0017  * @}
0018  */
0019
0020 #include "Geant4ExtraParticles.h"
0021 #include <DDG4/Geant4PhysicsConstructor.h>
0022 #include <DDG4/Geant4Kernel.h>
0023 #include <DDG4/Factories.h>
0024
0025 #include <G4ParticleTable.hh>
0026 #include <G4ParticleDefinition.hh>
0027 #include <G4PhysicalConstants.hh>
0028 #include <G4Version.hh>
0029
0030 #include <CLHEP/Units/SystemOfUnits.h>
0031 #include <CLHEP/Units/PhysicalConstants.h>
0032
0033 #include <fstream>
0034 #include <sstream>
0035 #include <string>
0036
0037 using namespace dd4hep::sim;
0038
0039 Geant4ExtraParticles::Geant4ExtraParticles(Geant4Context* ctxt, const std::string& nam)
0040   : Geant4PhysicsConstructor(ctxt, nam)
0041 {
0042   declareProperty("pdgfile", m_pdgfile);
0043 }
0044
0045 Geant4ExtraParticles::~Geant4ExtraParticles() {}
0046
0047
0048 void Geant4ExtraParticles::constructParticle(Constructor& ) {
0049   if (m_pdgfile.empty()) return;
0050
0051   info("pdgfile: %s",m_pdgfile.c_str());
0052
0053   G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0054   std::ifstream pdgFile( m_pdgfile.c_str(), std::ifstream::in );
0055
0056   if (!pdgFile.is_open()) {
0057     except("Could not open pdgfile: %s",m_pdgfile.c_str());
0058     return;
0059   }
0060
0061   info("opened pdgfile: %s",m_pdgfile.c_str());
0062
0063   while ( !pdgFile.eof() ) {
0064     // read line
0065     std::string linebuf;
0066     getline( pdgFile, linebuf );
0067
0068     // ignore comments
0069     if (linebuf.substr(0,1) == "#") continue;
0070     if (linebuf.substr(0,2) == "//") continue;
0071
0072     // ignore empty lines
0073     if (linebuf.empty()) continue;
0074
0075     // parse line
0076     int pdg;
0077     std::string nam;
0078     double charge;
0079     double mass;
0080     double width;
0081     double lifetime;
0082
0083     std::istringstream istr(linebuf);
0084
0085     istr >> pdg >> nam >> charge >> mass >> width >> lifetime;
0086
0087     // do add particles that don't fly
0088     // if (lifetime == 0) continue;
0089
0090     //unstable particles with too small lifetime have width -1
0091     bool stable = (width > 0 or width < -0.5) ? false : true;
0092
0093     if(width<0) width = 0;
0094
0095     // normalize to G4 units
0096     mass *= CLHEP::GeV;
0097
0098     if (charge != 0) {
0099       charge /= 3.;
0100     }
0101
0102     if (lifetime > 0) {
0103       lifetime = lifetime*CLHEP::mm/CLHEP::c_light;
0104     }
0105
0106     if (width == 0 && lifetime > 0) {
0107       width = CLHEP::hbar_Planck/lifetime;
0108     }
0109
0110     // don't add if the particle already exists
0111     G4ParticleDefinition* p = theParticleTable->FindParticle(pdg);
0112     if ( !p ) {
0113
0114       if (abs(pdg)>80 && abs(pdg)<=100) {
0115         // don't add generator-specific particles
0116       } else {
0117         /*
0118           if (pdg==5122) {
0119           G4cout << "Lambda_b0: " << "PDG=" << pdg << ", name=" << name << ", chrg=" << charge
0120           << ", mass=" << mass << ", width=" << width << ", lifetime=" << lifetime << "\n";
0121           G4cout << "debug: mass=" << 5.62 << ", width =" << 1.39e-12/6.582e-16 << ", lifetime=" << 1.39e-12 << "\n";
0122           }
0123         //*/
0124         p = new G4ParticleDefinition(
0125                                      nam,        // name
0126                                      mass,       // mass
0127                                      width,      // width
0128                                      charge,     // charge
0129                                      0,                                      // 2*spin
0130                                      0,          // parity
0131                                      0,          // C-conjugation
0132                                      0,          // 2*isospin
0133                                      0,          // 2*isospin3
0134                                      0,          // G-parity
0135                                      "extra",    // type
0136                                      0,          // lepton number
0137                                      0,          // baryon number
0138                                      pdg,        // PDG encoding
0139                                      stable,     // stable
0140                                      lifetime,   // lifetime
0141                                      NULL,       // decay table
0142                                      false);      // short lived
0143       }
0144     }
0145   }
0146
0147   G4cout << "Loaded extra particles using file: " << m_pdgfile << G4endl;
0148 }
0149
0150 void Geant4ExtraParticles::constructProcess(Constructor& ctor) {
0151   G4ParticleTable::G4PTblDicIterator* ParticleIterator = ctor.particleIterator();
0152   while((*ParticleIterator)()) {
0153     G4ParticleDefinition* pdef = ParticleIterator->value();
0154     G4ProcessManager* pmgr = pdef->GetProcessManager();
0155     if (pdef->GetParticleType() == "extra") {
0156       if (pdef->GetPDGCharge() != 0) {
0157         pmgr->AddProcess(new G4hMultipleScattering(), -1,  1, 1); //multiple scattering
0158         pmgr->AddProcess(new G4hIonisation(),  -1,  2, 2); // ionisation
0159       } else {
0160         //nothing to do
0161       }
0162     }
0163   }
0164 }
0165
0166 DECLARE_GEANT4ACTION(Geant4ExtraParticles)