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.

26 {}
CustomParticleFactory::~CustomParticleFactory ( )

Definition at line 28 of file CustomParticleFactory.cc.

28 {}

Member Function Documentation

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

Definition at line 90 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, TopDecayID::stable, str, AlCaHLTBitMon_QueryRunRegistry::string, and ApeEstimator_cff::width.

Referenced by getMassTable().

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)) {
104  pType = "rhadron";
105  }
106  if (CustomPDGParser::s_isSLepton(pdgCode)) {
107  pType = "sLepton";
108  }
109  if (CustomPDGParser::s_isMesonino(pdgCode)) {
110  pType = "mesonino";
111  }
112  if (CustomPDGParser::s_isSbaryon(pdgCode)) {
113  pType = "sbaryon";
114  }
115 
116  double massGeV = mass * CLHEP::GeV;
117  double width = 0.0;
118  double charge = CLHEP::eplus * CustomPDGParser::s_charge(pdgCode);
119  if (name.compare(0, 4, "~HIP") == 0) {
120  if ((name.compare(0, 7, "~HIPbar") == 0)) {
121  std::string str = name.substr(7);
122  charge = CLHEP::eplus * atoi(str.c_str()) / 3.;
123  } else {
124  std::string str = name.substr(4);
125  charge = -CLHEP::eplus * atoi(str.c_str()) / 3.;
126  }
127  }
128  if (name.compare(0, 9, "anti_~HIP") == 0) {
129  if ((name.compare(0, 12, "anti_~HIPbar") == 0)) {
130  std::string str = name.substr(12);
131  charge = -CLHEP::eplus * atoi(str.c_str()) / 3.;
132  } else {
133  std::string str = name.substr(9);
134  charge = CLHEP::eplus * atoi(str.c_str()) / 3.;
135  }
136  }
137  int spin = (int)CustomPDGParser::s_spin(pdgCode) - 1;
138  int parity = +1;
139  int conjugation = 0;
140  int isospin = 0;
141  int isospinZ = 0;
142  int gParity = 0;
143  int lepton = 0; //FIXME:
144  int baryon = 1; //FIXME:
145  bool stable = true;
146  double lifetime = -1;
147 
148  if (CustomPDGParser::s_isDphoton(pdgCode)) {
149  pType = "darkpho";
150  spin = 2;
151  parity = -1;
152  conjugation = -1;
153  isospin = 0;
154  isospinZ = 0;
155  gParity = 0;
156  lepton = 0;
157  baryon = 0;
158  stable = true;
159  lifetime = -1;
160  }
161 
162  G4DecayTable *decaytable = nullptr;
163  G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
164 
165  CustomParticle *particle = new CustomParticle(name,
166  massGeV,
167  width,
168  charge,
169  spin,
170  parity,
171  conjugation,
172  isospin,
173  isospinZ,
174  gParity,
175  pType,
176  lepton,
177  baryon,
178  pdgCode,
179  stable,
180  lifetime,
181  decaytable);
182 
183  if (pType == "rhadron" && name != "~g") {
184  G4String cloudname = name + "cloud";
185  G4String cloudtype = pType + "cloud";
186  spectator = theParticleTable->FindParticle(1000021);
187  spectatormass = spectator->GetPDGMass();
188  G4double cloudmass = massGeV - spectatormass;
189  CustomParticle *tmpParticle =
190  new CustomParticle(cloudname, cloudmass, 0.0, 0, 0, +1, 0, 0, 0, 0, cloudtype, 0, +1, 0, true, -1.0, nullptr);
191  particle->SetCloud(tmpParticle);
192  particle->SetSpectator(spectator);
193 
194  edm::LogInfo("SimG4CoreCustomPhysics")
195  << "CustomParticleFactory: " << name << " being assigned spectator" << spectator->GetParticleName()
196  << " and cloud " << cloudname << "\n"
197  << " Masses: " << mass << " Gev, " << spectatormass / CLHEP::GeV << " GeV and "
198  << cloudmass / CLHEP::GeV << " GeV.";
199  } else if (pType == "mesonino" || pType == "sbaryon") {
200  int sign = 1;
201  if (pdgCode < 0)
202  sign = -1;
203 
204  G4String cloudname = name + "cloud";
205  G4String cloudtype = pType + "cloud";
206  if (CustomPDGParser::s_isstopHadron(pdgCode)) {
207  spectator = theParticleTable->FindParticle(1000006 * sign);
208  } else {
209  if (CustomPDGParser::s_issbottomHadron(pdgCode)) {
210  spectator = theParticleTable->FindParticle(1000005 * sign);
211  } else {
212  spectator = nullptr;
213  edm::LogError("SimG4CoreCustomPhysics") << "CustomParticleFactory: Cannot find spectator parton";
214  }
215  }
216  if (spectator) {
217  spectatormass = spectator->GetPDGMass();
218  }
219  G4double cloudmass = massGeV - spectatormass;
220  CustomParticle *tmpParticle =
221  new CustomParticle(cloudname, cloudmass, 0.0, 0, 0, +1, 0, 0, 0, 0, cloudtype, 0, +1, 0, true, -1.0, nullptr);
222  particle->SetCloud(tmpParticle);
223  particle->SetSpectator(spectator);
224 
225  if (spectator)
226  edm::LogVerbatim("SimG4CoreCustomPhysics")
227  << "CustomParticleFactory: " << name << " being assigned spectator" << spectator->GetParticleName()
228  << " and cloud " << cloudname << "\n"
229  << " Masses: " << mass << " Gev, " << spectatormass / CLHEP::GeV << " GeV and "
230  << cloudmass / CLHEP::GeV << " GeV.";
231  } else {
232  particle->SetCloud(nullptr);
233  particle->SetSpectator(nullptr);
234  }
235  m_particles.push_back(particle);
236 }
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:10
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)
G4DecayTable * CustomParticleFactory::getAntiDecayTable ( int  pdgId,
G4DecayTable *  theDecayTable 
)
private

Definition at line 379 of file CustomParticleFactory.cc.

References mps_fire::i, dqmiolumiharvest::j, min(), Skims_PA_cff::name, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by loadCustomParticles().

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, theDecayChannel->GetBR(), nd, name[0], name[1], name[2], name[3]);
401  decaytable->Insert(aDecayChannel);
402  edm::LogInfo("SimG4CoreCustomPhysics")
403  << "CustomParticleFactory: inserted decay channel "
404  << " for pdgID= " << -pdgId << " " << parentName << " BR= " << theDecayChannel->GetBR()
405  << " Daughters: " << name[0] << " " << name[1] << " " << name[2] << " " << name[3];
406  }
407  return decaytable;
408 }
T min(T a, T b)
Definition: MathUtil.h:58
part
Definition: HCALResponse.h:20
const std::vector< G4ParticleDefinition * > & CustomParticleFactory::GetCustomParticles ( )

Definition at line 417 of file CustomParticleFactory.cc.

References m_particles.

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

Definition at line 318 of file CustomParticleFactory.cc.

References beamvalidation::br, mps_fire::i, mps_splice::line, LogDebug, Skims_PA_cff::name, EgammaValidation_cff::pdgId, AlCaHLTBitMon_QueryRunRegistry::string, and ToLower().

Referenced by loadCustomParticles().

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

Definition at line 238 of file CustomParticleFactory.cc.

References funct::abs(), addCustomParticle(), mps_splice::line, EgHLTOffHistBins_cfi::mass, SiStripPI::max, Skims_PA_cff::name, EgammaValidation_cff::pdgId, CustomPDGParser::s_isRHadron(), CustomPDGParser::s_isstopHadron(), Validation_hcalonly_cfi::sign, AlCaHLTBitMon_QueryRunRegistry::string, createJobs::tmp, and ToLower().

Referenced by loadCustomParticles().

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

Definition at line 30 of file CustomParticleFactory.cc.

References MillePedeFileConverter_cfg::e, dt4ml_dqm_sourceclient-live_cfg::filePath, getAntiDecayTable(), getDecayTable(), getMassTable(), GeV, mps_splice::line, loaded, MeV, EgammaValidation_cff::pdgId, alignCSCRings::s, AlCaHLTBitMon_QueryRunRegistry::string, ToLower(), and ApeEstimator_cff::width.

30  {
31  if (loaded) {
32  return;
33  }
34 #ifdef G4MULTITHREADED
35  G4MUTEXLOCK(&customParticleFactoryMutex);
36  if (loaded) {
37  return;
38  }
39 #endif
40 
41  // loading once
42  loaded = true;
43  std::ifstream configFile(filePath.c_str());
44 
46  edm::LogInfo("SimG4CoreCustomPhysics") << "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) == '#') {
53  continue;
54  } // Skip blank lines and comments.
55  // The mass table begins with a line containing "BLOCK MASS"
56  if (ToLower(line).find("block") < line.npos && ToLower(line).find("mass") < line.npos) {
57  edm::LogInfo("SimG4CoreCustomPhysics") << "CustomParticleFactory: Retrieving mass table.";
58  getMassTable(&configFile);
59  }
60  if (line.find("DECAY") < line.npos) {
61  int pdgId;
62  double width;
63  std::string tmpString;
64  std::stringstream lineStream(line);
65  lineStream >> tmpString >> pdgId >> width;
66  // assume SLHA format, e.g.: DECAY 1000021 5.50675438E+00 # gluino decays
67  edm::LogInfo("SimG4CoreCustomPhysics")
68  << "CustomParticleFactory: entry to G4DecayTable: pdgID, width " << pdgId << ", " << width;
69  G4ParticleDefinition *aParticle = theParticleTable->FindParticle(pdgId);
70  if (!aParticle || width == 0.0) {
71  continue;
72  }
73  G4DecayTable *aDecayTable = getDecayTable(&configFile, pdgId);
74  aParticle->SetDecayTable(aDecayTable);
75  aParticle->SetPDGStable(false);
76  aParticle->SetPDGLifeTime(1.0 / (width * CLHEP::GeV) * 6.582122e-22 * CLHEP::MeV * CLHEP::s);
77  G4ParticleDefinition *aAntiParticle = theParticleTable->FindAntiParticle(pdgId);
78  if (aAntiParticle && aAntiParticle->GetPDGEncoding() != pdgId) {
79  aAntiParticle->SetDecayTable(getAntiDecayTable(pdgId, aDecayTable));
80  aAntiParticle->SetPDGStable(false);
81  aAntiParticle->SetPDGLifeTime(1.0 / (width * CLHEP::GeV) * 6.582122e-22 * CLHEP::MeV * CLHEP::s);
82  }
83  }
84  }
85 #ifdef G4MULTITHREADED
86  G4MUTEXUNLOCK(&customParticleFactoryMutex);
87 #endif
88 }
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 410 of file CustomParticleFactory.cc.

References mps_fire::i, and str.

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

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

Member Data Documentation

bool CustomParticleFactory::loaded = false
staticprivate

Definition at line 27 of file CustomParticleFactory.h.

Referenced by loadCustomParticles().

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

Definition at line 28 of file CustomParticleFactory.h.

Referenced by addCustomParticle(), and GetCustomParticles().