Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Read generic slha file format. #7785

Merged
merged 2 commits into from Feb 19, 2015
Merged

Read generic slha file format. #7785

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7de89a5
allow to run on generic slha format
Feb 11, 2015
f3c710e
Change the end condition of getDecayTable()
Feb 11, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 2 additions & 0 deletions SimG4Core/CustomPhysics/interface/CustomParticleFactory.h
View file
Open in desktop
Expand Up @@ -27,6 +27,8 @@ class CustomParticleFactory {

static bool loaded;
static std::set<G4ParticleDefinition *> m_particles;

static std::string ToLower(std::string str);

};

Expand Down
231 changes: 140 additions & 91 deletions SimG4Core/CustomPhysics/src/CustomParticleFactory.cc
View file
Open in desktop
Expand Up @@ -5,6 +5,7 @@
#include <vector>
#include <sstream>
#include <set>
#include <locale>

#include "SimG4Core/CustomPhysics/interface/CustomPDGParser.h"
#include "SimG4Core/CustomPhysics/interface/CustomParticle.h"
Expand Down Expand Up @@ -36,27 +37,35 @@ void CustomParticleFactory::loadCustomParticles(const std::string & filePath){
std::ifstream configFile(filePath.c_str());

std::string line;
edm::LogInfo("CustomPhysics") << "Reading Custom Particle and G4DecayTable from " << filePath << std::endl;
// This should be compatible IMO to SLHA
G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
while(getline(configFile,line)){
if(line.find("PDG code")<line.npos) getMassTable(&configFile);
line.erase(0, line.find_first_not_of(" \t")); // Remove leading whitespace.
if (line.length()==0 || line.at(0) == '#') continue; // Skip blank lines and comments.
if (ToLower(line).find("block") < line.npos && // The mass table begins with a line containing "BLOCK MASS".
ToLower(line).find("mass") < line.npos) {
edm::LogInfo("CustomPhysics") << " Retrieving mass table." << std::endl;
getMassTable(&configFile);
}
if(line.find("DECAY")<line.npos){
int pdgId;
double width;
std::string tmpString;
std::stringstream lineStream(line);
lineStream>>tmpString>>pdgId>>width;
G4DecayTable* aDecayTable = getDecayTable(&configFile, pdgId);
G4ParticleDefinition *aParticle = theParticleTable->FindParticle(pdgId);
int pdgId;
double width;
std::string tmpString;
std::stringstream lineStream(line);
lineStream >> tmpString >> pdgId >> width; // assume SLHA format, e.g.: DECAY 1000021 5.50675438E+00 # gluino decays
edm::LogInfo("CustomPhysics") << "G4DecayTable: pdgID, width " << pdgId << ", " << width << std::endl;
G4DecayTable* aDecayTable = getDecayTable(&configFile, pdgId);
G4ParticleDefinition *aParticle = theParticleTable->FindParticle(pdgId);
G4ParticleDefinition *aAntiParticle = theParticleTable->FindAntiParticle(pdgId);
if(!aParticle) continue;
if (!aParticle) continue;
aParticle->SetDecayTable(aDecayTable);
aParticle->SetPDGStable(false);
aParticle->SetPDGLifeTime(1.0/(width*GeV)*6.582122e-22*MeV*s);
if(aAntiParticle && aAntiParticle->GetPDGEncoding()!=pdgId){
//aAntiParticle->SetDecayTable(getAntiDecayTable(pdgId,aDecayTable));
//aAntiParticle->SetDecayTable(getAntiDecayTable(pdgId,aDecayTable)); // Why is this commented out??
aAntiParticle->SetPDGStable(false);
aParticle->SetPDGLifeTime(1.0/(width*GeV)*6.582122e-22*MeV*s);
aAntiParticle->SetPDGLifeTime(1.0/(width*GeV)*6.582122e-22*MeV*s);
}
}
}
Expand All @@ -65,8 +74,7 @@ void CustomParticleFactory::loadCustomParticles(const std::string & filePath){

void CustomParticleFactory::addCustomParticle(int pdgCode, double mass, const std::string & name){


if(abs(pdgCode)%100 <14 && abs(pdgCode) / 1000000 == 0){
if(abs(pdgCode)%100 <14 && abs(pdgCode) / 1000000 == 0){
edm::LogError("") << "Pdg code too low " << pdgCode << " "<<abs(pdgCode) / 1000000 << std::endl;
return;
}
Expand All @@ -86,17 +94,17 @@ void CustomParticleFactory::addCustomParticle(int pdgCode, double mass, const st
double massGeV =mass*GeV;
double width = 0.0*MeV;
double charge = eplus* CustomPDGParser::s_charge(pdgCode);
if (name.compare(0,4,"~HIP") == 0)
if (name.compare(0,4,"~HIP") == 0)
{

if ((name.compare(0,7,"~HIPbar") == 0)) {std::string str = name.substr (7); charge=eplus*atoi(str.c_str())/3.;}
else {std::string str = name.substr (4); charge=eplus*atoi(str.c_str())*-1./3.; }
if ((name.compare(0,7,"~HIPbar") == 0)) {std::string str = name.substr (7); charge=eplus*atoi(str.c_str())/3.;}
else {std::string str = name.substr (4); charge=eplus*atoi(str.c_str())*-1./3.; }
}
if (name.compare(0,9,"anti_~HIP") == 0)
if (name.compare(0,9,"anti_~HIP") == 0)
{

if ((name.compare(0,12,"anti_~HIPbar") == 0)) {std::string str = name.substr (12); charge=eplus*atoi(str.c_str())*-1./3.;}
else {std::string str = name.substr (9); charge=eplus*atoi(str.c_str())*1./3.; }
if ((name.compare(0,12,"anti_~HIPbar") == 0)) {std::string str = name.substr (12); charge=eplus*atoi(str.c_str())*-1./3.;}
else {std::string str = name.substr (9); charge=eplus*atoi(str.c_str())*1./3.; }
}
int spin = (int)CustomPDGParser::s_spin(pdgCode)-1;
int parity = +1;
Expand Down Expand Up @@ -133,35 +141,30 @@ void CustomParticleFactory::addCustomParticle(int pdgCode, double mass, const st
particle->SetSpectator(spectator);

edm::LogInfo("CustomPhysics")<<name<<" being assigned "
<<particle->GetCloud()->GetParticleName()
<<" and "<<particle->GetSpectator()->GetParticleName()<<std::endl;
<<particle->GetCloud()->GetParticleName()
<<" and "<<particle->GetSpectator()->GetParticleName()<<std::endl;
edm::LogInfo("CustomPhysics")<<"Masses: "
<<particle->GetPDGMass()/GeV<<" Gev, "
<<particle->GetCloud()->GetPDGMass()/GeV<<" GeV and "
<<particle->GetSpectator()->GetPDGMass()/GeV<<" GeV."
<<std::endl;
}else if(pType == "mesonino" || pType == "sbaryon")
{
int sign=1;
if(pdgCode < 0 ) sign=-1;
<<particle->GetPDGMass()/GeV<<" Gev, "
<<particle->GetCloud()->GetPDGMass()/GeV<<" GeV and "
<<particle->GetSpectator()->GetPDGMass()/GeV<<" GeV."
<<std::endl;
} else if(pType == "mesonino" || pType == "sbaryon") {
int sign=1;
if(pdgCode < 0 ) sign=-1;

G4String cloudname = name+"cloud";
G4String cloudtype = pType+"cloud";
if(CustomPDGParser::s_isstopHadron(pdgCode))
{
spectator = theParticleTable->FindParticle(1000006*sign);
if(CustomPDGParser::s_isstopHadron(pdgCode)) {
spectator = theParticleTable->FindParticle(1000006*sign);
}
else
{
if(CustomPDGParser::s_issbottomHadron(pdgCode)) {
spectator = theParticleTable->FindParticle(1000005*sign);
}
else
{
else {
if (CustomPDGParser::s_issbottomHadron(pdgCode)) {
spectator = theParticleTable->FindParticle(1000005*sign);
} else {
spectator = 0;
edm::LogError("CustomPhysics")<< " Cannot find spectator parton";
}
}
}
spectatormass = spectator->GetPDGMass();
G4double cloudmass = mass-spectatormass/GeV;
CustomParticle *tmpParticle = new CustomParticle(
Expand All @@ -174,13 +177,13 @@ void CustomParticleFactory::addCustomParticle(int pdgCode, double mass, const st
particle->SetSpectator(spectator);

edm::LogInfo("CustomPhysics")<<name<<" being assigned "
<<particle->GetCloud()->GetParticleName()
<<" and "<<particle->GetSpectator()->GetParticleName()<<std::endl;
<<particle->GetCloud()->GetParticleName()
<<" and "<<particle->GetSpectator()->GetParticleName()<<std::endl;
edm::LogInfo("CustomPhysics")<<"Masses: "
<<particle->GetPDGMass()/GeV<<" Gev, "
<<particle->GetCloud()->GetPDGMass()/GeV<<" GeV and "
<<particle->GetSpectator()->GetPDGMass()/GeV<<" GeV."
<<std::endl;
<<particle->GetPDGMass()/GeV<<" Gev, "
<<particle->GetCloud()->GetPDGMass()/GeV<<" GeV and "
<<particle->GetSpectator()->GetPDGMass()/GeV<<" GeV."
<<std::endl;
}
else{
particle->SetCloud(0);
Expand All @@ -196,48 +199,74 @@ void CustomParticleFactory::getMassTable(std::ifstream *configFile) {
std::string name, tmp;
std::string line;
// This should be compatible IMO to SLHA
while(getline(*configFile,line))
{
if(tmp.find("Blo")<tmp.npos) break;
while (getline(*configFile,line)) {
line.erase(0, line.find_first_not_of(" \t")); // remove leading whitespace
if (line.length()==0 || line.at(0) == '#') continue; // skip blank lines and comments
if (ToLower(line).find("block") < line.npos) {
edm::LogInfo("CustomPhysics") << " Finished the Mass Table " << std::endl;
break;
}
std::stringstream sstr(line);
sstr >>pdgId>>mass>>tmp>>name;
sstr >> pdgId >> mass >> tmp >> name; // Assume SLHA format, e.g.: 1000001 5.68441109E+02 # ~d_L

addCustomParticle(pdgId, fabs(mass), name);
edm::LogInfo("CustomPhysics") << "Calling addCustomParticle for pdgId: " << pdgId
<< ", mass " << mass << ", name " << name
<< std::endl;
addCustomParticle(pdgId, fabs(mass), name);
////Find SM particle partner and check for the antiparticle.
int pdgIdPartner = pdgId%100;
G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition *aParticle = theParticleTable->FindParticle(pdgIdPartner);
//Add antiparticles for SUSY particles only, not for rHadrons.
if(aParticle && !CustomPDGParser::s_isRHadron(pdgId) && !CustomPDGParser::s_isstopHadron(pdgId)&& pdgId!=1000006 && pdgId!=-1000006 && pdgId!=25 && pdgId!=35 && pdgId!=36 && pdgId!=37){
int sign = aParticle->GetAntiPDGEncoding()/pdgIdPartner;
edm::LogInfo("CustomPhysics") << "Found aParticle = " << aParticle
<< ", pdgId = " << pdgId
<< ", pdgIdPartner = " << pdgIdPartner
<< ", CustomPDGParser::s_isRHadron(pdgId) = " << CustomPDGParser::s_isRHadron(pdgId)
<< ", CustomPDGParser::s_isstopHadron(pdgId) = " << CustomPDGParser::s_isstopHadron(pdgId)
<< std::endl;

if (aParticle &&
!CustomPDGParser::s_isRHadron(pdgId) &&
!CustomPDGParser::s_isstopHadron(pdgId) &&
pdgId!=1000006 &&
pdgId!=-1000006 &&
pdgId!=25 &&
pdgId!=35 &&
pdgId!=36 &&
pdgId!=37){
int sign = aParticle->GetAntiPDGEncoding()/pdgIdPartner;
edm::LogInfo("CustomPhysics") << "Found sign = " << sign
<< ", aParticle->GetAntiPDGEncoding() " << aParticle->GetAntiPDGEncoding()
<< ", pdgIdPartner = " << pdgIdPartner
<< std::endl;
if(abs(sign)!=1) {
std::cout<<"sgn: "<<sign<<" a "
<<aParticle->GetAntiPDGEncoding()
<<" b "<<pdgIdPartner
<<std::endl;
edm::LogInfo("CustomPhysics")<<"sgn: "<<sign<<" a "
<<aParticle->GetAntiPDGEncoding()
<<" b "<<pdgIdPartner
<<std::endl;
aParticle->DumpTable();
}
if(sign==-1 && pdgId!=25 && pdgId!=35 && pdgId!=36 && pdgId!=37 && pdgId!=1000039){
tmp = "anti_"+name;
addCustomParticle(-pdgId, mass, tmp);
theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(-pdgId);
}
else theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(pdgId);
tmp = "anti_"+name;
edm::LogInfo("CustomPhysics") << "Calling addCustomParticle for antiparticle with pdgId: " << -pdgId
<< ", mass " << mass << ", name " << tmp
<< std::endl;
addCustomParticle(-pdgId, mass, tmp);
theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(-pdgId);
}
else theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(pdgId);
}

if(pdgId==1000039) theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(pdgId);
if(pdgId==1000024 || pdgId==1000037 || pdgId==37) {
if(pdgId==1000039) theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(pdgId); // gravitino
if(pdgId==1000024 || pdgId==1000037 || pdgId==37) {
tmp = "anti_"+name;
edm::LogInfo("CustomPhysics") << "Calling addCustomParticle for antiparticle (2) with pdgId: " << -pdgId
<< ", mass " << mass << ", name " << tmp
<< std::endl;
addCustomParticle(-pdgId, mass, tmp);
theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(-pdgId);
}

/* getline(*configFile,tmp);
char text[100];
configFile->get(text,3);
tmp.clear();
tmp.append(text);
if(tmp.find("Bl")<tmp.npos) break;*/
}
}

Expand All @@ -246,49 +275,60 @@ G4DecayTable* CustomParticleFactory::getDecayTable(std::ifstream *configFile, i
double br;
int nDaughters;
std::vector<int> pdg(4);
std::string tmp;
std::string line;
std::vector<std::string> name(4);

G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();

std::string parentName = theParticleTable->FindParticle(pdgId)->GetParticleName();
G4DecayTable *decaytable= new G4DecayTable();

getline(*configFile,tmp);

while(!configFile->eof()){
while(getline(*configFile,line)) {

line.erase(0, line.find_first_not_of(" \t")); // remove leading whitespace
if (line.length()==0) continue; // skip blank lines
if (line.at(0) == '#' &&
ToLower(line).find("br") < line.npos &&
ToLower(line).find("nda") < line.npos) continue; // skip a comment of the form: # BR NDA ID1 ID2
if (line.at(0) == '#') { // other comments signal the end of the decay block
edm::LogInfo("CustomPhysics") << " Finished the Decay Table " << std::endl;
break;
}

pdg.clear();
name.clear();
(*configFile)>>br>>nDaughters;
for(int i=0;i<nDaughters;i++) (*configFile)>>pdg[i];
getline(*configFile,tmp);
for(int i=0;i<nDaughters;i++){
if(!theParticleTable->FindParticle(pdg[i])){
//std::cout<<pdg[i]<<" CustomParticleFactory::getDecayTable(): not found in the table!"<<std::endl;
continue;

std::stringstream sstr(line);
sstr >> br >> nDaughters; // assume SLHA format, e.g.: 1.49435135E-01 2 -15 16 # BR(H+ -> tau+ nu_tau)
edm::LogInfo("CustomPhysics") << " Branching Ratio: " << br << ", Number of Daughters: " << nDaughters << std::endl;
if (nDaughters > 4) {
edm::LogError("CustomPhysics") << "Number of daughters is too large (max = 4): " << nDaughters << " for pdgId: " << pdgId << std::endl;
break;
}
for(int i=0; i<nDaughters; i++) {
sstr >> pdg[i];
edm::LogInfo("CustomPhysics") << " Daughter ID " << pdg[i];
}
for (int i=0;i<nDaughters;i++) {
if (!theParticleTable->FindParticle(pdg[i])) {
edm::LogWarning("CustomPhysics")<<pdg[i]<<" CustomParticleFactory::getDecayTable(): not found in the table!"<<std::endl;
continue;
}
name[i] = theParticleTable->FindParticle(pdg[i])->GetParticleName();
}
////Set the G4 decay
////Set the G4 decay
G4PhaseSpaceDecayChannel *aDecayChannel = new G4PhaseSpaceDecayChannel(parentName, br, nDaughters,
name[0],name[1],name[2],name[3]);
decaytable->Insert(aDecayChannel);

/////////////////////////

char text[200];
configFile->get(text,2);
tmp.clear();
tmp.append(text);
if(tmp.find("#")<tmp.npos) break;
}

return decaytable;
}

G4DecayTable* CustomParticleFactory::getAntiDecayTable(int pdgId, G4DecayTable *theDecayTable) {

std::vector<std::string> name(4);
std::vector<std::string> name(4);
G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();

std::string parentName = theParticleTable->FindParticle(-pdgId)->GetParticleName();
Expand All @@ -311,3 +351,12 @@ G4DecayTable* CustomParticleFactory::getAntiDecayTable(int pdgId, G4DecayTable
}
return decaytable;
}


std::string CustomParticleFactory::ToLower(std::string str) {
std::locale loc;
for (std::string::size_type i=0; i<str.length(); ++i)
str.at(i) = std::tolower(str.at(i),loc);
return str;
}