CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Static Private Attributes
CustomParticleFactory Class Reference

#include <CustomParticleFactory.h>

Public Member Functions

 CustomParticleFactory ()
 
const std::vector< G4ParticleDefinition * > & GetCustomParticles ()
 
void loadCustomParticles (const std::string &filePath)
 
 ~CustomParticleFactory ()
 

Private Member Functions

void addCustomParticle (int pdgCode, double mass, const std::string &name)
 
G4DecayTable * getAntiDecayTable (int pdgId, G4DecayTable *theDecayTable)
 
G4DecayTable * getDecayTable (std::ifstream *configFile, int pdgId)
 
void getMassTable (std::ifstream *configFile)
 
std::string ToLower (std::string str)
 

Static Private Attributes

static bool loaded = false
 
static std::vector< G4ParticleDefinition * > m_particles
 

Detailed Description

Definition at line 12 of file CustomParticleFactory.h.

Constructor & Destructor Documentation

CustomParticleFactory::CustomParticleFactory ( )
explicit

Definition at line 26 of file CustomParticleFactory.cc.

27 {}
CustomParticleFactory::~CustomParticleFactory ( )

Definition at line 29 of file CustomParticleFactory.cc.

30 {}

Member Function Documentation

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

Definition at line 88 of file CustomParticleFactory.cc.

References funct::abs(), ALCARECOTkAlJpsiMuMu_cff::charge, GeV, createfilelist::int, m_particles, CustomPDGParser::s_charge(), CustomPDGParser::s_isDphoton(), CustomPDGParser::s_isMesonino(), CustomPDGParser::s_isRHadron(), CustomPDGParser::s_isSbaryon(), CustomPDGParser::s_issbottomHadron(), CustomPDGParser::s_isSLepton(), CustomPDGParser::s_isstopHadron(), CustomPDGParser::s_spin(), CustomParticle::SetCloud(), CustomParticle::SetSpectator(), Validation_hcalonly_cfi::sign, spin(), TopDecayID::stable, str, AlCaHLTBitMon_QueryRunRegistry::string, and ApeEstimator_cff::width.

Referenced by getMassTable().

89  {
90 
91  if(std::abs(pdgCode)%100 <14 && std::abs(pdgCode)/1000000 == 0){
92  edm::LogError("CustomParticleFactory::addCustomParticle")
93  << "Pdg code too low " << pdgCode << " "<<std::abs(pdgCode)/1000000;
94  return;
95  }
96 
98  G4String pType="custom";
99  G4String pSubType="";
100  G4double spectatormass = 0.0;
101  G4ParticleDefinition* spectator = nullptr;
103  if(CustomPDGParser::s_isRHadron(pdgCode)) { pType = "rhadron"; }
104  if(CustomPDGParser::s_isSLepton(pdgCode)) { pType = "sLepton"; }
105  if(CustomPDGParser::s_isMesonino(pdgCode)){ pType = "mesonino"; }
106  if(CustomPDGParser::s_isSbaryon(pdgCode)) { pType = "sbaryon"; }
107 
108  double massGeV = mass*CLHEP::GeV;
109  double width = 0.0;
110  double charge = CLHEP::eplus* CustomPDGParser::s_charge(pdgCode);
111  if (name.compare(0,4,"~HIP") == 0)
112  {
113  if ((name.compare(0,7,"~HIPbar") == 0)) {
114  std::string str = name.substr(7);
115  charge=CLHEP::eplus*atoi(str.c_str())/3.;
116  } else {
117  std::string str = name.substr(4);
118  charge=-CLHEP::eplus*atoi(str.c_str())/3.;
119  }
120  }
121  if (name.compare(0,9,"anti_~HIP") == 0)
122  {
123  if ((name.compare(0,12,"anti_~HIPbar") == 0)) {
124  std::string str = name.substr (12);
125  charge=-CLHEP::eplus*atoi(str.c_str())/3.;
126  } else {
127  std::string str = name.substr (9);
128  charge=CLHEP::eplus*atoi(str.c_str())/3.;
129  }
130  }
131  int spin = (int)CustomPDGParser::s_spin(pdgCode)-1;
132  int parity = +1;
133  int conjugation = 0;
134  int isospin = 0;
135  int isospinZ = 0;
136  int gParity = 0;
137  int lepton = 0; //FIXME:
138  int baryon = 1; //FIXME:
139  bool stable = true;
140  double lifetime = -1;
141 
142  if(CustomPDGParser::s_isDphoton(pdgCode)){
143  pType = "darkpho";
144  spin = 2;
145  parity = -1;
146  conjugation = -1;
147  isospin = 0;
148  isospinZ = 0;
149  gParity = 0;
150  lepton = 0;
151  baryon =0;
152  stable = true;
153  lifetime = -1;
154  }
155 
156  G4DecayTable *decaytable = nullptr;
157  G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
158 
159  CustomParticle *particle = new CustomParticle(name, massGeV, width, charge, spin,
160  parity, conjugation, isospin, isospinZ,
161  gParity, pType, lepton, baryon, pdgCode,
162  stable, lifetime, decaytable);
163 
164  if(pType == "rhadron" && name!="~g"){
165  G4String cloudname = name+"cloud";
166  G4String cloudtype = pType+"cloud";
167  spectator = theParticleTable->FindParticle(1000021);
168  spectatormass = spectator->GetPDGMass();
169  G4double cloudmass = massGeV - spectatormass;
170  CustomParticle *tmpParticle = new CustomParticle(
171  cloudname, cloudmass, 0.0, 0,
172  0, +1, 0,
173  0, 0, 0,
174  cloudtype, 0, +1, 0,
175  true, -1.0, nullptr );
176  particle->SetCloud(tmpParticle);
177  particle->SetSpectator(spectator);
178 
179  edm::LogInfo("SimG4CoreCustomPhysics")
180  <<"CustomParticleFactory: " <<name<<" being assigned spectator"
181  << spectator->GetParticleName() << " and cloud " <<cloudname << "\n"
182  <<" Masses: "<<mass<<" Gev, "
183  <<spectatormass/CLHEP::GeV<<" GeV and "<<cloudmass/CLHEP::GeV<<" GeV.";
184  } else if(pType == "mesonino" || pType == "sbaryon") {
185  int sign=1;
186  if(pdgCode < 0 ) sign=-1;
187 
188  G4String cloudname = name+"cloud";
189  G4String cloudtype = pType+"cloud";
190  if(CustomPDGParser::s_isstopHadron(pdgCode)) {
191  spectator = theParticleTable->FindParticle(1000006*sign);
192  }
193  else {
194  if (CustomPDGParser::s_issbottomHadron(pdgCode)) {
195  spectator = theParticleTable->FindParticle(1000005*sign);
196  } else {
197  spectator = nullptr;
198  edm::LogError("SimG4CoreCustomPhysics")<< "CustomParticleFactory: Cannot find spectator parton";
199  }
200  }
201  if(spectator) { spectatormass = spectator->GetPDGMass(); }
202  G4double cloudmass = massGeV - spectatormass;
203  CustomParticle *tmpParticle = new CustomParticle(
204  cloudname, cloudmass, 0.0, 0 ,
205  0, +1, 0,
206  0, 0, 0,
207  cloudtype, 0, +1, 0,
208  true, -1.0, nullptr);
209  particle->SetCloud(tmpParticle);
210  particle->SetSpectator(spectator);
211 
212  edm::LogInfo("SimG4CoreCustomPhysics")
213  <<"CustomParticleFactory: " <<name<<" being assigned spectator"
214  << spectator->GetParticleName() << " and cloud " <<cloudname << "\n"
215  <<" Masses: "<<mass<<" Gev, "
216  <<spectatormass/CLHEP::GeV<<" GeV and "<<cloudmass/CLHEP::GeV<<" GeV.";
217  }
218  else{
219  particle->SetCloud(nullptr);
220  particle->SetSpectator(nullptr);
221  }
222  m_particles.push_back(particle);
223 }
const double GeV
Definition: MathUtil.h:16
static bool s_isDphoton(int pdg)
static double s_charge(int pdg)
static const int stable
Definition: TopGenEvent.h:11
static bool s_isSbaryon(int pdg)
void SetCloud(G4ParticleDefinition *theCloud)
static std::vector< G4ParticleDefinition * > m_particles
static bool s_isMesonino(int pdg)
static bool s_isSLepton(int pdg)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void SetSpectator(G4ParticleDefinition *theSpectator)
static double s_spin(int pdg)
static bool s_isRHadron(int pdg)
static bool s_issbottomHadron(int pdg)
static bool s_isstopHadron(int pdg)
#define str(s)
float spin(float ph)
G4DecayTable * CustomParticleFactory::getAntiDecayTable ( int  pdgId,
G4DecayTable *  theDecayTable 
)
private

Definition at line 378 of file CustomParticleFactory.cc.

References mps_fire::i, min(), dataset::name, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by loadCustomParticles().

378  {
379 
380  std::string name[4] = {""};
381  G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
382 
383  std::string parentName = theParticleTable->FindParticle(-pdgId)->GetParticleName();
384  G4DecayTable *decaytable= new G4DecayTable();
385 
386  for(int i=0;i<theDecayTable->entries(); ++i){
387  G4VDecayChannel *theDecayChannel = theDecayTable->GetDecayChannel(i);
388  int nd = std::min(4, theDecayChannel->GetNumberOfDaughters());
389  for(int j=0; j<nd; ++j){
390  int id = theDecayChannel->GetDaughter(j)->GetAntiPDGEncoding();
391  const G4ParticleDefinition* part = theParticleTable->FindParticle(id);
392  if (!part) {
393  edm::LogWarning("SimG4CoreCustomPhysics")
394  <<"CustomParticleFactory: antiparticle with PDG code"<<id <<" not found!";
395  continue;
396  }
397  name[j] = part->GetParticleName();
398  }
399  G4PhaseSpaceDecayChannel *aDecayChannel =
400  new G4PhaseSpaceDecayChannel(parentName,
401  theDecayChannel->GetBR(), nd,
402  name[0],name[1],name[2],name[3]);
403  decaytable->Insert(aDecayChannel);
404  edm::LogInfo("SimG4CoreCustomPhysics") <<"CustomParticleFactory: inserted decay channel "
405  <<" for pdgID= " << -pdgId << " " << parentName
406  <<" BR= " << theDecayChannel->GetBR()
407  << " Daughters: " << name[0]
408  <<" " << name[1]<< " " << name[2]<< " " << name[3];
409  }
410  return decaytable;
411 }
T min(T a, T b)
Definition: MathUtil.h:58
part
Definition: HCALResponse.h:20
const std::vector< G4ParticleDefinition * > & CustomParticleFactory::GetCustomParticles ( )

Definition at line 420 of file CustomParticleFactory.cc.

References m_particles.

421 {
422  return m_particles;
423 }
static std::vector< G4ParticleDefinition * > m_particles
G4DecayTable * CustomParticleFactory::getDecayTable ( std::ifstream *  configFile,
int  pdgId 
)
private

Definition at line 314 of file CustomParticleFactory.cc.

References beamvalidation::br, mps_fire::i, geometryCSVtoXML::line, LogDebug, dataset::name, common_cff::pdgId, AlCaHLTBitMon_QueryRunRegistry::string, and ToLower().

Referenced by loadCustomParticles().

314  {
315 
316  double br;
317  int nDaughters;
318  int pdg[4] = {0};
320 
321  G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
322 
323  std::string parentName = theParticleTable->FindParticle(pdgId)->GetParticleName();
324  G4DecayTable *decaytable= new G4DecayTable();
325 
326  while(getline(*configFile,line)) {
327 
328  line.erase(0, line.find_first_not_of(" \t")); // remove leading whitespace
329  if (line.length()==0) continue; // skip blank lines
330  if (line.at(0) == '#' &&
331  ToLower(line).find("br") < line.npos &&
332  ToLower(line).find("nda") < line.npos) continue; // skip a comment of the form: # BR NDA ID1 ID2
333  if (line.at(0) == '#' || ToLower(line).find("block") < line.npos) {
334  edm::LogInfo("SimG4CoreCustomPhysics") <<"CustomParticleFactory: Finished the Decay Table ";
335  break;
336  }
337 
338  std::stringstream sstr(line);
339  sstr >> br >> nDaughters; // assume SLHA format, e.g.: 1.49435135E-01 2 -15 16 # BR(H+ -> tau+ nu_tau)
340  LogDebug("SimG4CoreCustomPhysics")
341  <<"CustomParticleFactory: Branching Ratio: " << br << ", Number of Daughters: " << nDaughters;
342  if (nDaughters > 4) {
343  edm::LogError("SimG4CoreCustomPhysics")
344  <<"CustomParticleFactory: Number of daughters is too large (max = 4): " << nDaughters
345  << " for pdgId: " << pdgId;
346  break;
347  }
348  if (br <= 0.0) {
349  edm::LogError("SimG4CoreCustomPhysics")
350  <<"CustomParticleFactory: Branching ratio is " << br
351  << " for pdgId: " << pdgId;
352  break;
353  }
354  std::string name[4] = {""};
355  for(int i=0; i<nDaughters; ++i) {
356  sstr >> pdg[i];
357  LogDebug("SimG4CoreCustomPhysics") <<"CustomParticleFactory: Daughter ID " << pdg[i];
358  const G4ParticleDefinition* part = theParticleTable->FindParticle(pdg[i]);
359  if (!part) {
360  edm::LogWarning("SimG4CoreCustomPhysics")
361  <<"CustomParticleFactory: particle with PDG code"<<pdg[i] <<" not found!";
362  continue;
363  }
364  name[i] = part->GetParticleName();
365  }
367  G4PhaseSpaceDecayChannel *aDecayChannel = new G4PhaseSpaceDecayChannel(parentName, br, nDaughters,
368  name[0],name[1],name[2],name[3]);
369  decaytable->Insert(aDecayChannel);
370  edm::LogInfo("SimG4CoreCustomPhysics") <<"CustomParticleFactory: inserted decay channel "
371  <<" for pdgID= " << pdgId << " " << parentName
372  <<" BR= " << br << " Daughters: " << name[0]
373  <<" " << name[1]<< " " << name[2]<< " " << name[3];
374  }
375  return decaytable;
376 }
#define LogDebug(id)
std::string ToLower(std::string str)
part
Definition: HCALResponse.h:20
void CustomParticleFactory::getMassTable ( std::ifstream *  configFile)
private

Definition at line 225 of file CustomParticleFactory.cc.

References funct::abs(), addCustomParticle(), geometryCSVtoXML::line, ResonanceBuilder::mass, SiStripPI::max, dataset::name, common_cff::pdgId, CustomPDGParser::s_isRHadron(), CustomPDGParser::s_isstopHadron(), Validation_hcalonly_cfi::sign, AlCaHLTBitMon_QueryRunRegistry::string, tmp, and ToLower().

Referenced by loadCustomParticles().

225  {
226 
227  int pdgId;
228  double mass;
231  G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
232 
233  // This should be compatible IMO to SLHA
234  while (getline(*configFile,line)) {
235  line.erase(0, line.find_first_not_of(" \t")); // remove leading whitespace
236  if (line.length()==0 || line.at(0) == '#') continue; // skip blank lines and comments
237  if (ToLower(line).find("block") < line.npos) {
238  edm::LogInfo("SimG4CoreCustomPhysics")
239  <<"CustomParticleFactory: Finished the Mass Table ";
240  break;
241  }
242  std::stringstream sstr(line);
243  sstr >> pdgId >> mass >> tmp >> name; // Assume SLHA format, e.g.: 1000001 5.68441109E+02 # ~d_L
244 
245  mass = std::max(mass, 0.0);
246  if(theParticleTable->FindParticle(pdgId)) { continue; }
247 
248  edm::LogInfo("SimG4CoreCustomPhysics")
249  <<"CustomParticleFactory: Calling addCustomParticle for pdgId: " << pdgId
250  << ", mass " << mass << " GeV " << name
251  << ", isRHadron: " << CustomPDGParser::s_isRHadron(pdgId)
252  << ", isstopHadron: " << CustomPDGParser::s_isstopHadron(pdgId);
253  addCustomParticle(pdgId, mass, name);
254 
256  int pdgIdPartner = pdgId%100;
257  G4ParticleDefinition *aParticle = theParticleTable->FindParticle(pdgIdPartner);
258  //Add antiparticles for SUSY particles only, not for rHadrons.
259  if(aParticle) {
260  edm::LogInfo("SimG4CoreCustomPhysics")
261  <<"CustomParticleFactory: Found partner particle for " << " pdgId= " << pdgId
262  << ", pdgIdPartner= " << pdgIdPartner << " " << aParticle->GetParticleName();
263  }
264 
265  if (aParticle &&
266  !CustomPDGParser::s_isRHadron(pdgId) &&
268  pdgId!=1000006 &&
269  pdgId!=-1000006 &&
270  pdgId!=25 &&
271  pdgId!=35 &&
272  pdgId!=36 &&
273  pdgId!=37){
274  int sign = aParticle->GetAntiPDGEncoding()/pdgIdPartner;
275  edm::LogInfo("SimG4CoreCustomPhysics")
276  <<"CustomParticleFactory: For " << aParticle->GetParticleName()
277  <<" pdg= " << pdgIdPartner
278  << ", GetAntiPDGEncoding() " << aParticle->GetAntiPDGEncoding()
279  << " sign= " << sign;
280 
281  if(std::abs(sign)!=1) {
282  edm::LogInfo("SimG4CoreCustomPhysics")
283  <<"CustomParticleFactory: sign= "<<sign<<" a: "
284  <<aParticle->GetAntiPDGEncoding()<<" b: "<<pdgIdPartner;
285  aParticle->DumpTable();
286  }
287  if(sign==-1 && pdgId!=25 && pdgId!=35 && pdgId!=36 && pdgId!=37 && pdgId!=1000039){
288  tmp = "anti_"+name;
289  if(!theParticleTable->FindParticle(-pdgId)) {
290  edm::LogInfo("SimG4CoreCustomPhysics")
291  <<"CustomParticleFactory: Calling addCustomParticle for antiparticle with pdgId: "
292  << -pdgId << ", mass " << mass << " GeV, name " << tmp;
293  addCustomParticle(-pdgId, mass, tmp);
294  theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(-pdgId);
295  }
296  }
297  else theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(pdgId);
298  }
299 
300  if(pdgId==1000039) theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(pdgId); // gravitino
301  if(pdgId==1000024 || pdgId==1000037 || pdgId==37) {
302  tmp = "anti_"+name;
303  if(!theParticleTable->FindParticle(-pdgId)) {
304  edm::LogInfo("SimG4CoreCustomPhysics")
305  <<"CustomParticleFactory: Calling addCustomParticle for antiparticle (2) with pdgId: "
306  << -pdgId << ", mass " << mass << " GeV, name " << tmp;
307  addCustomParticle(-pdgId, mass, tmp);
308  theParticleTable->FindParticle(pdgId)->SetAntiPDGEncoding(-pdgId);
309  }
310  }
311  }
312 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::string ToLower(std::string str)
void addCustomParticle(int pdgCode, double mass, const std::string &name)
static bool s_isRHadron(int pdg)
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
static bool s_isstopHadron(int pdg)
void CustomParticleFactory::loadCustomParticles ( const std::string &  filePath)

Definition at line 32 of file CustomParticleFactory.cc.

References MillePedeFileConverter_cfg::e, getAntiDecayTable(), getDecayTable(), getMassTable(), GeV, geometryCSVtoXML::line, loaded, MeV, common_cff::pdgId, alignCSCRings::s, AlCaHLTBitMon_QueryRunRegistry::string, ToLower(), and ApeEstimator_cff::width.

32  {
33 
34  if(loaded) { return; }
35 #ifdef G4MULTITHREADED
36  G4MUTEXLOCK(&customParticleFactoryMutex);
37  if(loaded) { return; }
38 #endif
39 
40  // loading once
41  loaded = true;
42  std::ifstream configFile(filePath.c_str());
43 
45  edm::LogInfo("SimG4CoreCustomPhysics")
46  <<"CustomParticleFactory: Reading Custom Particle and G4DecayTable from \n"
47  << filePath;
48  // This should be compatible IMO to SLHA
49  G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
50  while(getline(configFile, line)){
51  line.erase(0, line.find_first_not_of(" \t")); // Remove leading whitespace.
52  if (line.length()==0 || line.at(0) == '#') { continue; } // Skip blank lines and comments.
53  // The mass table begins with a line containing "BLOCK MASS"
54  if (ToLower(line).find("block") < line.npos &&
55  ToLower(line).find("mass") < line.npos) {
56  edm::LogInfo("SimG4CoreCustomPhysics") <<"CustomParticleFactory: Retrieving mass table.";
57  getMassTable(&configFile);
58  }
59  if(line.find("DECAY")<line.npos) {
60  int pdgId;
61  double width;
62  std::string tmpString;
63  std::stringstream lineStream(line);
64  lineStream >> tmpString >> pdgId >> width;
65  // assume SLHA format, e.g.: DECAY 1000021 5.50675438E+00 # gluino decays
66  edm::LogInfo("SimG4CoreCustomPhysics")
67  <<"CustomParticleFactory: entry to G4DecayTable: pdgID, width "
68  << pdgId << ", " << width;
69  G4ParticleDefinition *aParticle = theParticleTable->FindParticle(pdgId);
70  if (!aParticle || width == 0.0) { continue; }
71  G4DecayTable* aDecayTable = getDecayTable(&configFile, pdgId);
72  aParticle->SetDecayTable(aDecayTable);
73  aParticle->SetPDGStable(false);
74  aParticle->SetPDGLifeTime(1.0/(width*CLHEP::GeV)*6.582122e-22*CLHEP::MeV*CLHEP::s);
75  G4ParticleDefinition *aAntiParticle = theParticleTable->FindAntiParticle(pdgId);
76  if(aAntiParticle && aAntiParticle->GetPDGEncoding() != pdgId){
77  aAntiParticle->SetDecayTable(getAntiDecayTable(pdgId,aDecayTable));
78  aAntiParticle->SetPDGStable(false);
79  aAntiParticle->SetPDGLifeTime(1.0/(width*CLHEP::GeV)*6.582122e-22*CLHEP::MeV*CLHEP::s);
80  }
81  }
82  }
83 #ifdef G4MULTITHREADED
84  G4MUTEXUNLOCK(&customParticleFactoryMutex);
85 #endif
86 }
const double GeV
Definition: MathUtil.h:16
const double MeV
std::string ToLower(std::string str)
void getMassTable(std::ifstream *configFile)
G4DecayTable * getDecayTable(std::ifstream *configFile, int pdgId)
G4DecayTable * getAntiDecayTable(int pdgId, G4DecayTable *theDecayTable)
std::string CustomParticleFactory::ToLower ( std::string  str)
private

Definition at line 413 of file CustomParticleFactory.cc.

References mps_fire::i, create_public_lumi_plots::loc, and str.

Referenced by getDecayTable(), getMassTable(), and loadCustomParticles().

413  {
414  std::locale loc;
415  for (std::string::size_type i=0; i<str.length(); ++i)
416  str.at(i) = std::tolower(str.at(i),loc);
417  return str;
418 }
uint16_t size_type
#define str(s)

Member Data Documentation

bool CustomParticleFactory::loaded = false
staticprivate

Definition at line 30 of file CustomParticleFactory.h.

Referenced by loadCustomParticles().

std::vector< G4ParticleDefinition * > CustomParticleFactory::m_particles
staticprivate

Definition at line 31 of file CustomParticleFactory.h.

Referenced by addCustomParticle(), and GetCustomParticles().