Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:24 

0001 #include "SimG4Core/CustomPhysics/interface/CustomParticleFactory.h"
0002 #include "SimG4Core/CustomPhysics/interface/CustomPDGParser.h"
0003 #include "SimG4Core/CustomPhysics/interface/CustomParticle.h"
0004 
0005 #include "FWCore/ParameterSet/interface/FileInPath.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 
0008 #include "G4ParticleTable.hh"
0009 #include "G4DecayTable.hh"
0010 #include "G4PhaseSpaceDecayChannel.hh"
0011 #include "G4ProcessManager.hh"
0012 #include "G4ParticleDefinition.hh"
0013 #include <CLHEP/Units/SystemOfUnits.h>
0014 
0015 #include "SimG4Core/CustomPhysics/interface/CMSSIMP.h"
0016 #include "SimG4Core/CustomPhysics/interface/CMSAntiSIMP.h"
0017 
0018 #include <iomanip>
0019 #include <iostream>
0020 #include <sstream>
0021 
0022 bool CustomParticleFactory::loaded = false;
0023 std::vector<G4ParticleDefinition *> CustomParticleFactory::m_particles;
0024 
0025 #ifdef G4MULTITHREADED
0026 G4Mutex CustomParticleFactory::customParticleFactoryMutex = G4MUTEX_INITIALIZER;
0027 #endif
0028 
0029 CustomParticleFactory::CustomParticleFactory() {}
0030 
0031 void CustomParticleFactory::loadCustomParticles(const std::string &filePath) {
0032   if (loaded) {
0033     return;
0034   }
0035 #ifdef G4MULTITHREADED
0036   G4MUTEXLOCK(&customParticleFactoryMutex);
0037   if (loaded) {
0038     return;
0039   }
0040 #endif
0041 
0042   // loading once
0043   loaded = true;
0044   std::ifstream configFile(filePath.c_str());
0045 
0046   std::string line;
0047   edm::LogVerbatim("SimG4CoreCustomPhysics")
0048       << "CustomParticleFactory: Reading Custom Particle and G4DecayTable from \n"
0049       << filePath;
0050   // This should be compatible IMO to SLHA
0051   G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0052   while (getline(configFile, line)) {
0053     line.erase(0, line.find_first_not_of(" \t"));  // Remove leading whitespace.
0054     if (line.length() == 0 || line.at(0) == '#') {
0055       continue;
0056     }  // Skip blank lines and comments.
0057     // The mass table begins with a line containing "BLOCK MASS"
0058     if (ToLower(line).find("block") < line.npos && ToLower(line).find("mass") < line.npos) {
0059       edm::LogVerbatim("SimG4CoreCustomPhysics") << "CustomParticleFactory: Retrieving mass table.";
0060       getMassTable(&configFile);
0061     }
0062     if (line.find("DECAY") < line.npos) {
0063       int pdgId;
0064       double width;
0065       std::string tmpString;
0066       std::stringstream lineStream(line);
0067       lineStream >> tmpString >> pdgId >> width;
0068       // assume SLHA format, e.g.: DECAY  1000021  5.50675438E+00   # gluino decays
0069       edm::LogVerbatim("SimG4CoreCustomPhysics")
0070           << "CustomParticleFactory: entry to G4DecayTable: pdgID, width " << pdgId << ",  " << width;
0071       G4ParticleDefinition *aParticle = theParticleTable->FindParticle(pdgId);
0072       if (nullptr == aParticle || width == 0.0) {
0073         continue;
0074       }
0075       G4DecayTable *aDecayTable = getDecayTable(&configFile, pdgId);
0076       aParticle->SetDecayTable(aDecayTable);
0077       aParticle->SetPDGStable(false);
0078       aParticle->SetPDGLifeTime(1.0 / (width * CLHEP::GeV) * 6.582122e-22 * CLHEP::MeV * CLHEP::s);
0079       G4ParticleDefinition *aAntiParticle = theParticleTable->FindAntiParticle(pdgId);
0080       if (nullptr != aAntiParticle && aAntiParticle->GetPDGEncoding() != pdgId) {
0081         aAntiParticle->SetDecayTable(getAntiDecayTable(pdgId, aDecayTable));
0082         aAntiParticle->SetPDGStable(false);
0083         aAntiParticle->SetPDGLifeTime(1.0 / (width * CLHEP::GeV) * 6.582122e-22 * CLHEP::MeV * CLHEP::s);
0084       }
0085     }
0086   }
0087 #ifdef G4MULTITHREADED
0088   G4MUTEXUNLOCK(&customParticleFactoryMutex);
0089 #endif
0090 }
0091 
0092 void CustomParticleFactory::addCustomParticle(int pdgCode, double mass, const std::string &name) {
0093   if (std::abs(pdgCode) % 100 < 14 && std::abs(pdgCode) / 1000000 == 0) {
0094     edm::LogError("CustomParticleFactory::addCustomParticle")
0095         << "Pdg code too low " << pdgCode << " " << std::abs(pdgCode) / 1000000;
0096     return;
0097   }
0098 
0099   if (CustomPDGParser::s_isSIMP(pdgCode)) {
0100     CMSSIMP *simp = CMSSIMP::Definition(mass * CLHEP::GeV);
0101     CMSAntiSIMP *antisimp = CMSAntiSIMP::Definition(mass * CLHEP::GeV);
0102     m_particles.push_back(simp);
0103     m_particles.push_back(antisimp);
0104     return;
0105   }
0106 
0107   /////////////////////// Check!!!!!!!!!!!!!
0108   G4String pType = "custom";
0109   G4String pSubType = "";
0110   G4double spectatormass = 0.0;
0111   G4ParticleDefinition *spectator = nullptr;
0112   //////////////////////
0113   if (CustomPDGParser::s_isgluinoHadron(pdgCode)) {
0114     pType = "rhadron";
0115   }
0116   if (CustomPDGParser::s_isSLepton(pdgCode)) {
0117     pType = "sLepton";
0118   }
0119   if (CustomPDGParser::s_isMesonino(pdgCode)) {
0120     pType = "mesonino";
0121   }
0122   if (CustomPDGParser::s_isSbaryon(pdgCode)) {
0123     pType = "sbaryon";
0124   }
0125 
0126   double massGeV = mass * CLHEP::GeV;
0127   double width = 0.0;
0128   double charge = CLHEP::eplus * CustomPDGParser::s_charge(pdgCode);
0129   if (name.compare(0, 4, "~HIP") == 0) {
0130     if ((name.compare(0, 7, "~HIPbar") == 0)) {
0131       std::string str = name.substr(7);
0132       charge = CLHEP::eplus * atoi(str.c_str()) / 3.;
0133     } else {
0134       std::string str = name.substr(4);
0135       charge = -CLHEP::eplus * atoi(str.c_str()) / 3.;
0136     }
0137   }
0138   if (name.compare(0, 9, "anti_~HIP") == 0) {
0139     if ((name.compare(0, 12, "anti_~HIPbar") == 0)) {
0140       std::string str = name.substr(12);
0141       charge = -CLHEP::eplus * atoi(str.c_str()) / 3.;
0142     } else {
0143       std::string str = name.substr(9);
0144       charge = CLHEP::eplus * atoi(str.c_str()) / 3.;
0145     }
0146   }
0147   int spin = (int)CustomPDGParser::s_spin(pdgCode) - 1;
0148   int parity = +1;
0149   int conjugation = 0;
0150   int isospin = 0;
0151   int isospinZ = 0;
0152   int gParity = 0;
0153   int lepton = 0;  //FIXME:
0154   int baryon = 1;  //FIXME:
0155   bool stable = true;
0156   double lifetime = -1;
0157 
0158   if (CustomPDGParser::s_isDphoton(pdgCode)) {
0159     pType = "darkpho";
0160     spin = 2;
0161     parity = -1;
0162     conjugation = -1;
0163     isospin = 0;
0164     isospinZ = 0;
0165     gParity = 0;
0166     lepton = 0;
0167     baryon = 0;
0168     stable = true;
0169     lifetime = -1;
0170   }
0171 
0172   G4DecayTable *decaytable = nullptr;
0173   G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0174 
0175   CustomParticle *particle = new CustomParticle(name,
0176                                                 massGeV,
0177                                                 width,
0178                                                 charge,
0179                                                 spin,
0180                                                 parity,
0181                                                 conjugation,
0182                                                 isospin,
0183                                                 isospinZ,
0184                                                 gParity,
0185                                                 pType,
0186                                                 lepton,
0187                                                 baryon,
0188                                                 pdgCode,
0189                                                 stable,
0190                                                 lifetime,
0191                                                 decaytable);
0192 
0193   if (pType == "rhadron" && name != "~g") {
0194     G4String cloudname = name + "cloud";
0195     G4String cloudtype = pType + "cloud";
0196     spectator = theParticleTable->FindParticle(1000021);
0197     spectatormass = spectator->GetPDGMass();
0198     G4double cloudmass = massGeV - spectatormass;
0199     CustomParticle *tmpParticle =
0200         new CustomParticle(cloudname, cloudmass, 0.0, 0, 0, +1, 0, 0, 0, 0, cloudtype, 0, +1, 0, true, -1.0, nullptr);
0201     particle->SetCloud(tmpParticle);
0202     particle->SetSpectator(spectator);
0203 
0204     edm::LogVerbatim("SimG4CoreCustomPhysics")
0205         << "CustomParticleFactory: " << name << " being assigned spectator" << spectator->GetParticleName()
0206         << " and cloud " << cloudname << "\n"
0207         << "                        Masses: " << mass << " Gev, " << spectatormass / CLHEP::GeV << " GeV and "
0208         << cloudmass / CLHEP::GeV << " GeV.";
0209   } else if (pType == "mesonino" || pType == "sbaryon") {
0210     int sign = 1;
0211     if (pdgCode < 0)
0212       sign = -1;
0213 
0214     G4String cloudname = name + "cloud";
0215     G4String cloudtype = pType + "cloud";
0216     if (CustomPDGParser::s_isstopHadron(pdgCode)) {
0217       spectator = theParticleTable->FindParticle(1000006 * sign);
0218     } else {
0219       if (CustomPDGParser::s_issbottomHadron(pdgCode)) {
0220         spectator = theParticleTable->FindParticle(1000005 * sign);
0221       } else {
0222         spectator = nullptr;
0223         edm::LogError("SimG4CoreCustomPhysics") << "CustomParticleFactory: Cannot find spectator parton";
0224       }
0225     }
0226     if (nullptr != spectator) {
0227       spectatormass = spectator->GetPDGMass();
0228     }
0229     G4double cloudmass = massGeV - spectatormass;
0230     CustomParticle *tmpParticle =
0231         new CustomParticle(cloudname, cloudmass, 0.0, 0, 0, +1, 0, 0, 0, 0, cloudtype, 0, +1, 0, true, -1.0, nullptr);
0232     particle->SetCloud(tmpParticle);
0233     particle->SetSpectator(spectator);
0234 
0235     if (nullptr != spectator)
0236       edm::LogVerbatim("SimG4CoreCustomPhysics")
0237           << "CustomParticleFactory: " << name << " being assigned spectator" << spectator->GetParticleName()
0238           << " and cloud " << cloudname << "\n"
0239           << "                        Masses: " << mass << " Gev, " << spectatormass / CLHEP::GeV << " GeV and "
0240           << cloudmass / CLHEP::GeV << " GeV.";
0241   } else {
0242     particle->SetCloud(nullptr);
0243     particle->SetSpectator(nullptr);
0244   }
0245   m_particles.push_back(particle);
0246 }
0247 
0248 void CustomParticleFactory::getMassTable(std::ifstream *configFile) {
0249   int pdgId;
0250   double mass;
0251   std::string name, tmp;
0252   std::string line;
0253   G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0254 
0255   // This should be compatible IMO to SLHA
0256   while (getline(*configFile, line)) {
0257     line.erase(0, line.find_first_not_of(" \t"));  // remove leading whitespace
0258     if (line.length() == 0 || line.at(0) == '#')
0259       continue;  // skip blank lines and comments
0260     if (ToLower(line).find("block") < line.npos) {
0261       edm::LogInfo("SimG4CoreCustomPhysics") << "CustomParticleFactory: Finished the Mass Table ";
0262       break;
0263     }
0264     std::stringstream sstr(line);
0265     sstr >> pdgId >> mass >> tmp >> name;  // Assume SLHA format, e.g.: 1000001 5.68441109E+02 # ~d_L
0266 
0267     mass = std::max(mass, 0.0);
0268     if (theParticleTable->FindParticle(pdgId)) {
0269       continue;
0270     }
0271 
0272     edm::LogVerbatim("SimG4CoreCustomPhysics")
0273         << "CustomParticleFactory: Calling addCustomParticle for pdgId: " << pdgId << ", mass " << mass << " GeV  "
0274         << name << ", isgluinoHadron: " << CustomPDGParser::s_isgluinoHadron(pdgId)
0275         << ", isstopHadron: " << CustomPDGParser::s_isstopHadron(pdgId)
0276         << ", isSIMP: " << CustomPDGParser::s_isSIMP(pdgId);
0277     addCustomParticle(pdgId, mass, name);
0278 
0279     ////Find SM particle partner and check for the antiparticle.
0280     int pdgIdPartner = pdgId % 100;
0281     G4ParticleDefinition *aParticle = theParticleTable->FindParticle(pdgIdPartner);
0282     //Add antiparticles for SUSY particles only, not for rHadrons.
0283     if (nullptr != aParticle) {
0284       edm::LogVerbatim("SimG4CoreCustomPhysics")
0285           << "CustomParticleFactory: Found partner particle for "
0286           << " pdgId= " << pdgId << ", pdgIdPartner= " << pdgIdPartner << " " << aParticle->GetParticleName();
0287     }
0288 
0289     if (nullptr != aParticle && !CustomPDGParser::s_isgluinoHadron(pdgId) && !CustomPDGParser::s_isstopHadron(pdgId) &&
0290         !CustomPDGParser::s_isSIMP(pdgId) && pdgId != 1000006 && pdgId != -1000006 && pdgId != 25 && pdgId != 35 &&
0291         pdgId != 36 && pdgId != 37) {
0292       int sign = aParticle->GetAntiPDGEncoding() / pdgIdPartner;
0293       edm::LogVerbatim("SimG4CoreCustomPhysics")
0294           << "CustomParticleFactory: For " << aParticle->GetParticleName() << " pdg= " << pdgIdPartner
0295           << ", GetAntiPDGEncoding() " << aParticle->GetAntiPDGEncoding() << " sign= " << sign;
0296 
0297       if (std::abs(sign) != 1) {
0298         edm::LogVerbatim("SimG4CoreCustomPhysics")
0299             << "CustomParticleFactory: sign= " << sign << " a: " << aParticle->GetAntiPDGEncoding()
0300             << " b: " << pdgIdPartner;
0301         aParticle->DumpTable();
0302       }
0303       if (sign == -1 && pdgId != 25 && pdgId != 35 && pdgId != 36 && pdgId != 37 && pdgId != 1000039) {
0304         tmp = "anti_" + name;
0305         if (theParticleTable->FindParticle(-pdgId) == nullptr) {
0306           edm::LogVerbatim("SimG4CoreCustomPhysics")
0307               << "CustomParticleFactory: Calling addCustomParticle for antiparticle with pdgId: " << -pdgId << ", mass "
0308               << mass << " GeV, name " << tmp;
0309           addCustomParticle(-pdgId, mass, tmp);
0310           theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(-pdgId);
0311         }
0312       } else
0313         theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(pdgId);
0314     }
0315 
0316     if (pdgId == 1000039)
0317       theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(pdgId);  // gravitino
0318     if (pdgId == 1000024 || pdgId == 1000037 || pdgId == 37) {
0319       tmp = "anti_" + name;
0320       if (!theParticleTable->FindParticle(-pdgId)) {
0321         edm::LogVerbatim("SimG4CoreCustomPhysics")
0322             << "CustomParticleFactory: Calling addCustomParticle for antiparticle (2) with pdgId: " << -pdgId
0323             << ", mass " << mass << " GeV, name " << tmp;
0324         addCustomParticle(-pdgId, mass, tmp);
0325         theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(-pdgId);
0326       }
0327     }
0328     if (pdgId == 9000006) {
0329       tmp = name + "bar";
0330       edm::LogVerbatim("CustomPhysics") << "Calling addCustomParticle for antiparticle with pdgId: " << -pdgId
0331                                         << ", mass " << mass << ", name " << tmp;
0332       addCustomParticle(-pdgId, mass, tmp);
0333       theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(-pdgId);
0334     }
0335   }
0336 }
0337 
0338 G4DecayTable *CustomParticleFactory::getDecayTable(std::ifstream *configFile, int pdgId) {
0339   double br;
0340   int nDaughters;
0341   int pdg[4] = {0};
0342   std::string line;
0343 
0344   G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0345 
0346   std::string parentName = theParticleTable->FindParticle(pdgId)->GetParticleName();
0347   G4DecayTable *decaytable = new G4DecayTable();
0348 
0349   while (getline(*configFile, line)) {
0350     line.erase(0, line.find_first_not_of(" \t"));  // remove leading whitespace
0351     if (line.length() == 0)
0352       continue;  // skip blank lines
0353     if (line.at(0) == '#' && ToLower(line).find("br") < line.npos && ToLower(line).find("nda") < line.npos)
0354       continue;  // skip a comment of the form:  # BR  NDA  ID1  ID2
0355     if (line.at(0) == '#' || ToLower(line).find("block") < line.npos) {
0356       edm::LogInfo("SimG4CoreCustomPhysics") << "CustomParticleFactory: Finished the Decay Table ";
0357       break;
0358     }
0359 
0360     std::stringstream sstr(line);
0361     sstr >> br >> nDaughters;  // assume SLHA format, e.g.:  1.49435135E-01  2  -15  16  # BR(H+ -> tau+ nu_tau)
0362     LogDebug("SimG4CoreCustomPhysics") << "CustomParticleFactory: Branching Ratio: " << br
0363                                        << ", Number of Daughters: " << nDaughters;
0364     if (nDaughters > 4) {
0365       edm::LogError("SimG4CoreCustomPhysics")
0366           << "CustomParticleFactory: Number of daughters is too large (max = 4): " << nDaughters
0367           << " for pdgId: " << pdgId;
0368       break;
0369     }
0370     if (br <= 0.0) {
0371       edm::LogError("SimG4CoreCustomPhysics")
0372           << "CustomParticleFactory: Branching ratio is " << br << " for pdgId: " << pdgId;
0373       break;
0374     }
0375     std::string name[4] = {""};
0376     for (int i = 0; i < nDaughters; ++i) {
0377       sstr >> pdg[i];
0378       LogDebug("SimG4CoreCustomPhysics") << "CustomParticleFactory: Daughter ID " << pdg[i];
0379       const G4ParticleDefinition *part = theParticleTable->FindParticle(pdg[i]);
0380       if (!part) {
0381         edm::LogWarning("SimG4CoreCustomPhysics")
0382             << "CustomParticleFactory: particle with PDG code" << pdg[i] << " not found!";
0383         continue;
0384       }
0385       name[i] = part->GetParticleName();
0386     }
0387     ////Set the G4 decay
0388     G4PhaseSpaceDecayChannel *aDecayChannel =
0389         new G4PhaseSpaceDecayChannel(parentName, br, nDaughters, name[0], name[1], name[2], name[3]);
0390     decaytable->Insert(aDecayChannel);
0391     edm::LogVerbatim("SimG4CoreCustomPhysics")
0392         << "CustomParticleFactory: inserted decay channel "
0393         << " for pdgID= " << pdgId << "  " << parentName << " BR= " << br << " Daughters: " << name[0] << "  "
0394         << name[1] << "  " << name[2] << "  " << name[3];
0395   }
0396   return decaytable;
0397 }
0398 
0399 G4DecayTable *CustomParticleFactory::getAntiDecayTable(int pdgId, G4DecayTable *theDecayTable) {
0400   std::string name[4] = {""};
0401   G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0402 
0403   std::string parentName = theParticleTable->FindParticle(-pdgId)->GetParticleName();
0404   G4DecayTable *decaytable = new G4DecayTable();
0405 
0406   for (int i = 0; i < theDecayTable->entries(); ++i) {
0407     G4VDecayChannel *theDecayChannel = theDecayTable->GetDecayChannel(i);
0408     int nd = std::min(4, theDecayChannel->GetNumberOfDaughters());
0409     for (int j = 0; j < nd; ++j) {
0410       int id = theDecayChannel->GetDaughter(j)->GetAntiPDGEncoding();
0411       G4ParticleDefinition *part = theParticleTable->FindParticle(id);
0412       if (nullptr == part) {
0413         edm::LogWarning("SimG4CoreCustomPhysics")
0414             << "CustomParticleFactory: antiparticle with PDG code" << id << " not found!";
0415         continue;
0416       }
0417       name[j] = part->GetParticleName();
0418     }
0419     G4PhaseSpaceDecayChannel *aDecayChannel =
0420         new G4PhaseSpaceDecayChannel(parentName, theDecayChannel->GetBR(), nd, name[0], name[1], name[2], name[3]);
0421     decaytable->Insert(aDecayChannel);
0422     edm::LogVerbatim("SimG4CoreCustomPhysics")
0423         << "CustomParticleFactory: inserted decay channel "
0424         << " for pdgID= " << -pdgId << "  " << parentName << " BR= " << theDecayChannel->GetBR()
0425         << " Daughters: " << name[0] << "  " << name[1] << "  " << name[2] << "  " << name[3];
0426   }
0427   return decaytable;
0428 }
0429 
0430 std::string CustomParticleFactory::ToLower(std::string str) {
0431   std::locale loc;
0432   for (std::string::size_type i = 0; i < str.length(); ++i)
0433     str.at(i) = std::tolower(str.at(i), loc);
0434   return str;
0435 }
0436 
0437 const std::vector<G4ParticleDefinition *> &CustomParticleFactory::getCustomParticles() { return m_particles; }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!