CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
G4ProcessHelper Class Reference

#include <G4ProcessHelper.h>

Public Member Functions

G4bool ApplicabilityTester (const G4ParticleDefinition &aPart)
 
 G4ProcessHelper (const edm::ParameterSet &p, CustomParticleFactory *ptr)
 
 G4ProcessHelper (const G4ProcessHelper &)=delete
 
ReactionProduct GetFinalState (const G4Track &aTrack, G4ParticleDefinition *&aTarget)
 
G4double GetInclusiveCrossSection (const G4DynamicParticle *aParticle, const G4Element *anElement)
 
G4ProcessHelperoperator= (const G4ProcessHelper &)=delete
 
 ~G4ProcessHelper ()
 

Private Member Functions

G4double PhaseSpace (const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle)
 
G4double Pom (const double boost)
 
G4bool ReactionGivesBaryon (const ReactionProduct &aReaction)
 
G4bool ReactionIsPossible (const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle)
 
G4double ReactionProductMass (const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle)
 
void ReadAndParse (const G4String &str, std::vector< G4String > &tokens, const G4String &delimiters=" ")
 
G4double Regge (const double boost)
 

Private Attributes

double amplitude
 
G4double checkfraction
 
double ek_0
 
CustomParticleFactoryfParticleFactory
 
double gamma
 
TProfile * h_q_gamma
 
TProfile * h_q_p
 
TH1D * h_sqrts
 
TProfile * h_xsec_cms
 
TProfile * h_xsec_lab
 
double hadronlifetime
 
std::map< const
G4ParticleDefinition *, G4bool > 
known_particles
 
double mixing
 
G4int n_22
 
G4int n_23
 
ReactionMap nReactionMap
 
std::map< G4String, G4double > parameters
 
G4ParticleTable * particleTable
 
ReactionMap pReactionMap
 
bool reggemodel
 
bool resonant
 
double suppressionfactor
 
HistoHelper * theHistoHelper
 
G4ParticleDefinition * theNeutron
 
G4ParticleDefinition * theProton
 
G4ParticleDefinition * theRbaryoncloud
 
ReactionMaptheReactionMap
 
G4ParticleDefinition * theRmesoncloud
 
G4ParticleDefinition * theTarget
 

Detailed Description

Definition at line 26 of file G4ProcessHelper.h.

Constructor & Destructor Documentation

G4ProcessHelper::G4ProcessHelper ( const edm::ParameterSet p,
CustomParticleFactory ptr 
)

Definition at line 19 of file G4ProcessHelper.cc.

References personalPlayback::fp, edm::FileInPath::fullPath(), edm::ParameterSet::getParameter(), GeV, mps_fire::i, geometryCSVtoXML::line, alignCSCRings::s, AlCaHLTBitMon_QueryRunRegistry::string, and filterCSVwithJSON::target.

19  {
20  fParticleFactory = ptr;
21 
22  particleTable = G4ParticleTable::GetParticleTable();
23 
24  theProton = particleTable->FindParticle("proton");
25  theNeutron = particleTable->FindParticle("neutron");
26 
27  G4String line;
28 
29  edm::FileInPath fp = p.getParameter<edm::FileInPath>("processesDef");
30  std::string processDefFilePath = fp.fullPath();
31  std::ifstream process_stream(processDefFilePath.c_str());
32 
33  resonant = p.getParameter<bool>("resonant");
34  ek_0 = p.getParameter<double>("resonanceEnergy") * GeV;
35  gamma = p.getParameter<double>("gamma") * GeV;
36  amplitude = p.getParameter<double>("amplitude") * millibarn;
37  suppressionfactor = p.getParameter<double>("reggeSuppression");
38  hadronlifetime = p.getParameter<double>("hadronLifeTime");
39  reggemodel = p.getParameter<bool>("reggeModel");
40  mixing = p.getParameter<double>("mixing");
41 
42  edm::LogInfo("SimG4CoreCustomPhysics") << "ProcessHelper: Read in physics parameters:"
43  << "\n Resonant = " << resonant << "\n ResonanceEnergy = " << ek_0 / GeV
44  << " GeV"
45  << "\n Gamma = " << gamma / GeV << " GeV"
46  << "\n Amplitude = " << amplitude / millibarn << " millibarn"
47  << "ReggeSuppression = " << 100 * suppressionfactor << " %"
48  << "HadronLifeTime = " << hadronlifetime << " s"
49  << "ReggeModel = " << reggemodel << "Mixing = " << mixing * 100 << " %";
50 
51  checkfraction = 0;
52  n_22 = 0;
53  n_23 = 0;
54 
55  while (getline(process_stream, line)) {
56  std::vector<G4String> tokens;
57  //Getting a line
58  ReadAndParse(line, tokens, "#");
59  //Important info
60  G4String incident = tokens[0];
61 
62  G4ParticleDefinition* incidentDef = particleTable->FindParticle(incident);
63  //particleTable->DumpTable();
64  G4int incidentPDG = incidentDef->GetPDGEncoding();
65  known_particles[incidentDef] = true;
66 
67  G4String target = tokens[1];
68  edm::LogInfo("SimG4CoreCustomPhysics") << "ProcessHelper: Incident " << incident << "; Target " << target;
69 
70  // Making a ReactionProduct
71  ReactionProduct prod;
72  for (size_t i = 2; i != tokens.size(); i++) {
73  G4String part = tokens[i];
74  if (particleTable->contains(part)) {
75  prod.push_back(particleTable->FindParticle(part)->GetPDGEncoding());
76  } else {
77  G4Exception("G4ProcessHelper",
78  "UnkownParticle",
79  FatalException,
80  "Initialization: The reaction product list contained an unknown particle");
81  }
82  }
83  if (target == "proton") {
84  pReactionMap[incidentPDG].push_back(prod);
85  } else if (target == "neutron") {
86  nReactionMap[incidentPDG].push_back(prod);
87  } else {
88  G4Exception("G4ProcessHelper",
89  "IllegalTarget",
90  FatalException,
91  "Initialization: The reaction product list contained an illegal target particle");
92  }
93  }
94 
95  process_stream.close();
96 
97  for (auto part : fParticleFactory->getCustomParticles()) {
98  CustomParticle* particle = dynamic_cast<CustomParticle*>(part);
99  if (particle) {
100  edm::LogInfo("SimG4CoreCustomPhysics") << "ProcessHelper: Lifetime of " << part->GetParticleName() << " set to "
101  << particle->GetPDGLifeTime() / s << " s;"
102  << " isStable: " << particle->GetPDGStable();
103  }
104  }
105 }
const double GeV
Definition: MathUtil.h:16
G4ParticleTable * particleTable
G4ParticleDefinition * theProton
G4ParticleDefinition * theNeutron
std::vector< G4int > ReactionProduct
void ReadAndParse(const G4String &str, std::vector< G4String > &tokens, const G4String &delimiters=" ")
double suppressionfactor
Log< level::Info, false > LogInfo
ReactionMap pReactionMap
ReactionMap nReactionMap
std::map< const G4ParticleDefinition *, G4bool > known_particles
part
Definition: HCALResponse.h:20
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
const std::vector< G4ParticleDefinition * > & getCustomParticles()
CustomParticleFactory * fParticleFactory
std::string fullPath() const
Definition: FileInPath.cc:161
G4double checkfraction
G4ProcessHelper::~G4ProcessHelper ( )

Definition at line 107 of file G4ProcessHelper.cc.

107 {}
G4ProcessHelper::G4ProcessHelper ( const G4ProcessHelper )
delete

Member Function Documentation

G4bool G4ProcessHelper::ApplicabilityTester ( const G4ParticleDefinition &  aPart)

Definition at line 109 of file G4ProcessHelper.cc.

Referenced by FullModelHadronicProcess::IsApplicable().

109  {
110  const G4ParticleDefinition* aP = &aPart;
111  if (known_particles[aP])
112  return true;
113  return false;
114 }
std::map< const G4ParticleDefinition *, G4bool > known_particles
ReactionProduct G4ProcessHelper::GetFinalState ( const G4Track &  aTrack,
G4ParticleDefinition *&  aTarget 
)

Definition at line 187 of file G4ProcessHelper.cc.

References mps_fire::i, CustomPDGParser::s_isMesonino(), CustomPDGParser::s_isRBaryon(), CustomPDGParser::s_isRMeson(), CustomPDGParser::s_isSbaryon(), PixelPluginsPhase0_cfi::select, and findQualityFiles::size.

Referenced by FullModelHadronicProcess::PostStepDoIt().

187  {
188  const G4DynamicParticle* aDynamicParticle = aTrack.GetDynamicParticle();
189 
190  //-----------------------------------------------
191  // Choose n / p as target
192  // and get ReactionProductList pointer
193  //-----------------------------------------------
194 
195  G4Material* aMaterial = aTrack.GetMaterial();
196  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
197  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
198 
199  G4double NumberOfProtons = 0;
200  G4double NumberOfNucleons = 0;
201 
202  for (size_t elm = 0; elm < aMaterial->GetNumberOfElements(); elm++) {
203  //Summing number of protons per unit volume
204  NumberOfProtons += NbOfAtomsPerVolume[elm] * (*theElementVector)[elm]->GetZ();
205  //Summing nucleons (not neutrons)
206  NumberOfNucleons += NbOfAtomsPerVolume[elm] * (*theElementVector)[elm]->GetN();
207  }
208 
209  if (G4UniformRand() < NumberOfProtons / NumberOfNucleons) {
212  } else {
215  }
216  aTarget = theTarget;
217 
218  G4int theIncidentPDG = aDynamicParticle->GetDefinition()->GetPDGEncoding();
219 
220  if (reggemodel && CustomPDGParser::s_isMesonino(theIncidentPDG) && G4UniformRand() * mixing > 0.5 &&
221  aDynamicParticle->GetDefinition()->GetPDGCharge() == 0.) {
222  // G4cout<<"Oscillating..."<<G4endl;
223  theIncidentPDG *= -1;
224  }
225 
226  bool baryonise = false;
227 
228  if (reggemodel && G4UniformRand() > 0.9 &&
229  ((CustomPDGParser::s_isMesonino(theIncidentPDG) && theIncidentPDG > 0) ||
230  CustomPDGParser::s_isRMeson(theIncidentPDG)))
231  baryonise = true;
232 
233  //Making a pointer directly to the ReactionProductList we are looking at. Makes life easier :-)
234  ReactionProductList* aReactionProductList = &((*theReactionMap)[theIncidentPDG]);
235 
236  //-----------------------------------------------
237  // Count processes
238  // kinematic check
239  // compute number of 2 -> 2 and 2 -> 3 processes
240  //-----------------------------------------------
241 
242  G4int N22 = 0; //Number of 2 -> 2 processes
243  G4int N23 = 0; //Number of 2 -> 3 processes. Couldn't think of more informative names
244 
245  //This is the list to be populated
246  ReactionProductList theReactionProductList;
247  std::vector<bool> theChargeChangeList;
248 
249  for (ReactionProductList::iterator prod_it = aReactionProductList->begin(); prod_it != aReactionProductList->end();
250  prod_it++) {
251  G4int secondaries = prod_it->size();
252  // If the reaction is not possible we will not consider it
253  if (ReactionIsPossible(*prod_it, aDynamicParticle) &&
254  (!reggemodel || (baryonise && ReactionGivesBaryon(*prod_it)) ||
255  (!baryonise && !ReactionGivesBaryon(*prod_it)) || (CustomPDGParser::s_isSbaryon(theIncidentPDG)) ||
256  (CustomPDGParser::s_isRBaryon(theIncidentPDG)))) {
257  // The reaction is possible. Let's store and count it
258  theReactionProductList.push_back(*prod_it);
259  if (secondaries == 2) {
260  N22++;
261  } else if (secondaries == 3) {
262  N23++;
263  } else {
264  G4cerr << "ReactionProduct has unsupported number of secondaries: " << secondaries << G4endl;
265  }
266  } /*else {
267  edm::LogInfo("SimG4CoreCustomPhysics")<<"There was an impossible process"<<G4endl;
268  }*/
269  }
270  // edm::LogInfo("SimG4CoreCustomPhysics")<<"The size of the ReactionProductList is: "<<theReactionProductList.size()<<G4endl;
271 
272  if (theReactionProductList.empty())
273  G4Exception("G4ProcessHelper",
274  "NoProcessPossible",
275  FatalException,
276  "GetFinalState: No process could be selected from the given list.");
277 
278  // For the Regge model no phase space considerations. We pick a process at random
279  if (reggemodel) {
280  int n_rps = theReactionProductList.size();
281  int select = (int)(G4UniformRand() * n_rps);
282  // G4cout<<"Possible: "<<n_rps<<", chosen: "<<select<<G4endl;
283  return theReactionProductList[select];
284  }
285 
286  // Fill a probability map. Remember total probability
287  // 2->2 is 0.15*1/n_22 2->3 uses phase space
288  G4double p22 = 0.15;
289  G4double p23 = 1 - p22; // :-)
290 
291  std::vector<G4double> Probabilities;
292  std::vector<G4bool> TwotoThreeFlag;
293 
294  G4double CumulatedProbability = 0;
295 
296  // To each ReactionProduct we assign a cumulated probability and a flag
297  // discerning between 2 -> 2 and 2 -> 3
298  for (unsigned int i = 0; i != theReactionProductList.size(); i++) {
299  if (theReactionProductList[i].size() == 2) {
300  CumulatedProbability += p22 / N22;
301  TwotoThreeFlag.push_back(false);
302  } else {
303  CumulatedProbability += p23 / N23;
304  TwotoThreeFlag.push_back(true);
305  }
306  Probabilities.push_back(CumulatedProbability);
307  // edm::LogInfo("SimG4CoreCustomPhysics")<<"Pushing back cumulated probability: "<<CumulatedProbability<<G4endl;
308  }
309 
310  //Renormalising probabilities
311  // edm::LogInfo("SimG4CoreCustomPhysics")<<"Probs: ";
312  for (std::vector<G4double>::iterator it = Probabilities.begin(); it != Probabilities.end(); it++) {
313  *it /= CumulatedProbability;
314  // edm::LogInfo("SimG4CoreCustomPhysics")<<*it<<" ";
315  }
316  // edm::LogInfo("SimG4CoreCustomPhysics")<<G4endl;
317 
318  // Choosing ReactionProduct
319 
320  G4bool selected = false;
321  G4int tries = 0;
322  // ReactionProductList::iterator prod_it;
323 
324  //Keep looping over the list until we have a choice, or until we have tried 100 times
325  unsigned int i;
326  while (!selected && tries < 100) {
327  i = 0;
328  G4double dice = G4UniformRand();
329  // edm::LogInfo("SimG4CoreCustomPhysics")<<"What's the dice?"<<dice<<G4endl;
330  while (dice > Probabilities[i] && i < theReactionProductList.size()) {
331  // edm::LogInfo("SimG4CoreCustomPhysics")<<"i: "<<i<<G4endl;
332  i++;
333  }
334 
335  // edm::LogInfo("SimG4CoreCustomPhysics")<<"Chosen i: "<<i<<G4endl;
336 
337  if (!TwotoThreeFlag[i]) {
338  // 2 -> 2 processes are chosen immediately
339  selected = true;
340  } else {
341  // 2 -> 3 processes require a phase space lookup
342  if (PhaseSpace(theReactionProductList[i], aDynamicParticle) > G4UniformRand())
343  selected = true;
344  //selected = true;
345  }
346  // double suppressionfactor=0.5;
347  if (selected && particleTable->FindParticle(theReactionProductList[i][0])->GetPDGCharge() !=
348  aDynamicParticle->GetDefinition()->GetPDGCharge()) {
349  /*
350  edm::LogInfo("SimG4CoreCustomPhysics")<<"Incoming particle "<<aDynamicParticle->GetDefinition()->GetParticleName()
351  <<" has charge "<<aDynamicParticle->GetDefinition()->GetPDGCharge()<<G4endl;
352  edm::LogInfo("SimG4CoreCustomPhysics")<<"Suggested particle "<<particleTable->FindParticle(theReactionProductList[i][0])->GetParticleName()
353  <<" has charge "<<particleTable->FindParticle(theReactionProductList[i][0])->GetPDGCharge()<<G4endl;
354  */
355  if (G4UniformRand() < suppressionfactor)
356  selected = false;
357  }
358  tries++;
359  // edm::LogInfo("SimG4CoreCustomPhysics")<<"Tries: "<<tries<<G4endl;
360  }
361  if (tries >= 100)
362  G4cerr << "Could not select process!!!!" << G4endl;
363 
364  // edm::LogInfo("SimG4CoreCustomPhysics")<<"So far so good"<<G4endl;
365  // edm::LogInfo("SimG4CoreCustomPhysics")<<"Sec's: "<<theReactionProductList[i].size()<<G4endl;
366 
367  //Updating checkfraction:
368  if (theReactionProductList[i].size() == 2) {
369  n_22++;
370  } else {
371  n_23++;
372  }
373 
374  checkfraction = (1.0 * n_22) / (n_22 + n_23);
375  // edm::LogInfo("SimG4CoreCustomPhysics")<<"n_22: "<<n_22<<" n_23: "<<n_23<<" Checkfraction: "<<checkfraction<<G4endl;
376  // edm::LogInfo("SimG4CoreCustomPhysics") <<"Biig number: "<<n_22+n_23<<G4endl;
377  //Return the chosen ReactionProduct
378  return theReactionProductList[i];
379 }
G4ParticleTable * particleTable
G4ParticleDefinition * theProton
static bool s_isSbaryon(int pdg)
G4ParticleDefinition * theNeutron
static bool s_isRMeson(int pdg)
G4double PhaseSpace(const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle)
static bool s_isMesonino(int pdg)
double suppressionfactor
static bool s_isRBaryon(int pdg)
ReactionMap pReactionMap
ReactionMap nReactionMap
ReactionMap * theReactionMap
G4double checkfraction
G4ParticleDefinition * theTarget
G4bool ReactionIsPossible(const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle)
std::vector< ReactionProduct > ReactionProductList
tuple size
Write out results.
G4bool ReactionGivesBaryon(const ReactionProduct &aReaction)
G4double G4ProcessHelper::GetInclusiveCrossSection ( const G4DynamicParticle *  aParticle,
const G4Element *  anElement 
)

Definition at line 116 of file G4ProcessHelper.cc.

References mergeVDriftHistosByStation::name, funct::pow(), dttmaxenums::R, CustomPDGParser::s_containedQuarks(), CustomPDGParser::s_isMesonino(), CustomPDGParser::s_isRBaryon(), CustomPDGParser::s_isRGlueball(), CustomPDGParser::s_isRMeson(), CustomPDGParser::s_isSbaryon(), and mathSSE::sqrt().

Referenced by FullModelHadronicProcess::GetMicroscopicCrossSection().

116  {
117  //We really do need a dedicated class to handle the cross sections. They might not always be constant
118 
119  //Disassemble the PDG-code
120 
121  G4int thePDGCode = aParticle->GetDefinition()->GetPDGEncoding();
122  double boost = (aParticle->GetKineticEnergy() + aParticle->GetMass()) / aParticle->GetMass();
123  // G4cout<<"thePDGCode: "<<thePDGCode<<G4endl;
124  G4double theXsec = 0;
125  G4String name = aParticle->GetDefinition()->GetParticleName();
126  if (!reggemodel) {
127  //Flat cross section
128  if (CustomPDGParser::s_isRGlueball(thePDGCode)) {
129  theXsec = 24 * millibarn;
130  } else {
131  std::vector<G4int> nq = CustomPDGParser::s_containedQuarks(thePDGCode);
132  // edm::LogInfo("SimG4CoreCustomPhysics")<<"Number of quarks: "<<nq.size()<<G4endl;
133  for (std::vector<G4int>::iterator it = nq.begin(); it != nq.end(); it++) {
134  // edm::LogInfo("SimG4CoreCustomPhysics")<<"Quarkvector: "<<*it<<G4endl;
135  if (*it == 1 || *it == 2)
136  theXsec += 12 * millibarn;
137  if (*it == 3)
138  theXsec += 6 * millibarn;
139  }
140  }
141  } else { //reggemodel
142  double R = Regge(boost);
143  double P = Pom(boost);
144  if (thePDGCode > 0) {
145  if (CustomPDGParser::s_isMesonino(thePDGCode))
146  theXsec = (P + R) * millibarn;
147  if (CustomPDGParser::s_isSbaryon(thePDGCode))
148  theXsec = 2 * P * millibarn;
149  if (CustomPDGParser::s_isRMeson(thePDGCode) || CustomPDGParser::s_isRGlueball(thePDGCode))
150  theXsec = (R + 2 * P) * millibarn;
151  if (CustomPDGParser::s_isRBaryon(thePDGCode))
152  theXsec = 3 * P * millibarn;
153  } else {
154  if (CustomPDGParser::s_isMesonino(thePDGCode))
155  theXsec = P * millibarn;
156  if (CustomPDGParser::s_isSbaryon(thePDGCode))
157  theXsec = (2 * (P + R) + 30 / sqrt(boost)) * millibarn;
158  if (CustomPDGParser::s_isRMeson(thePDGCode) || CustomPDGParser::s_isRGlueball(thePDGCode))
159  theXsec = (R + 2 * P) * millibarn;
160  if (CustomPDGParser::s_isRBaryon(thePDGCode))
161  theXsec = 3 * P * millibarn;
162  }
163  }
164  //Adding resonance
165  if (resonant) {
166  double e_0 = ek_0 + aParticle->GetDefinition()->GetPDGMass(); //Now total energy
167 
168  e_0 = sqrt(aParticle->GetDefinition()->GetPDGMass() * aParticle->GetDefinition()->GetPDGMass() +
169  theProton->GetPDGMass() * theProton->GetPDGMass() + 2. * e_0 * theProton->GetPDGMass());
170  // edm::LogInfo("SimG4CoreCustomPhysics")<<e_0/GeV<<G4endl;
171  // edm::LogInfo("SimG4CoreCustomPhysics")<<ek_0/GeV<<" "<<aParticle->GetDefinition()->GetPDGMass()/GeV<<" "<<theProton->GetPDGMass()/GeV<<G4endl;
172  double sqrts = sqrt(aParticle->GetDefinition()->GetPDGMass() * aParticle->GetDefinition()->GetPDGMass() +
173  theProton->GetPDGMass() * theProton->GetPDGMass() +
174  2 * aParticle->GetTotalEnergy() * theProton->GetPDGMass());
175 
176  double res_result = amplitude * (gamma * gamma / 4.) /
177  ((sqrts - e_0) * (sqrts - e_0) + (gamma * gamma / 4.)); //Non-relativistic Breit Wigner
178 
179  theXsec += res_result;
180  // if(fabs(aParticle->GetKineticEnergy()/GeV-200)<10) std::cout<<sqrts/GeV<<" "<<theXsec/millibarn<<std::endl;
181  }
182 
183  // std::cout<<"Xsec/nucleon: "<<theXsec/millibarn<<"millibarn, Total Xsec: "<<theXsec * anElement->GetN()/millibarn<<" millibarn"<<std::endl;
184  return theXsec * pow(anElement->GetN(), 0.7) * 1.25; // * 0.523598775598299;
185 }
static bool s_isRGlueball(int pdg)
G4ParticleDefinition * theProton
static bool s_isSbaryon(int pdg)
static bool s_isRMeson(int pdg)
static bool s_isMesonino(int pdg)
T sqrt(T t)
Definition: SSEVec.h:19
static bool s_isRBaryon(int pdg)
static std::vector< int > s_containedQuarks(int pdg)
std::pair< OmniClusterRef, TrackingParticleRef > P
G4double Regge(const double boost)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
G4double Pom(const double boost)
G4ProcessHelper& G4ProcessHelper::operator= ( const G4ProcessHelper )
delete
G4double G4ProcessHelper::PhaseSpace ( const ReactionProduct aReaction,
const G4DynamicParticle *  aDynamicParticle 
)
private

Definition at line 416 of file G4ProcessHelper.cc.

References GeV, funct::pow(), and mathSSE::sqrt().

416  {
417  G4double qValue = ReactionProductMass(aReaction, aDynamicParticle);
418  G4double phi = sqrt(1 + qValue / (2 * 0.139 * GeV)) * pow(qValue / (1.1 * GeV), 3. / 2.);
419  return (phi / (1 + phi));
420 }
G4double ReactionProductMass(const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle)
const double GeV
Definition: MathUtil.h:16
T sqrt(T t)
Definition: SSEVec.h:19
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
double G4ProcessHelper::Pom ( const double  boost)
private

Definition at line 451 of file G4ProcessHelper.cc.

References a, b, c, ztail::d, funct::pow(), and mathSSE::sqrt().

451  {
452  double a = 4.138224000651535;
453  double b = 1.50377557581421;
454  double c = -0.05449742257808247;
455  double d = 0.0008221235048211401;
456  return a + b * sqrt(boost) + c * boost + d * pow(boost, 1.5);
457 }
const edm::EventSetup & c
tuple d
Definition: ztail.py:151
T sqrt(T t)
Definition: SSEVec.h:19
double b
Definition: hdecay.h:118
double a
Definition: hdecay.h:119
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
G4bool G4ProcessHelper::ReactionGivesBaryon ( const ReactionProduct aReaction)
private

Definition at line 409 of file G4ProcessHelper.cc.

References CustomPDGParser::s_isRBaryon(), and CustomPDGParser::s_isSbaryon().

409  {
410  for (ReactionProduct::const_iterator it = aReaction.begin(); it != aReaction.end(); it++)
412  return true;
413  return false;
414 }
static bool s_isSbaryon(int pdg)
static bool s_isRBaryon(int pdg)
G4bool G4ProcessHelper::ReactionIsPossible ( const ReactionProduct aReaction,
const G4DynamicParticle *  aDynamicParticle 
)
private

Definition at line 402 of file G4ProcessHelper.cc.

403  {
404  if (ReactionProductMass(aReaction, aDynamicParticle) > 0)
405  return true;
406  return false;
407 }
G4double ReactionProductMass(const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle)
G4double G4ProcessHelper::ReactionProductMass ( const ReactionProduct aReaction,
const G4DynamicParticle *  aDynamicParticle 
)
private

Definition at line 381 of file G4ProcessHelper.cc.

References mathSSE::sqrt().

382  {
383  // Incident energy:
384  G4double E_incident = aDynamicParticle->GetTotalEnergy();
385  //edm::LogInfo("SimG4CoreCustomPhysics")<<"Total energy: "<<E_incident<<" Kinetic: "<<aDynamicParticle->GetKineticEnergy()<<G4endl;
386  // sqrt(s)= sqrt(m_1^2 + m_2^2 + 2 E_1 m_2)
387  G4double m_1 = aDynamicParticle->GetDefinition()->GetPDGMass();
388  G4double m_2 = theTarget->GetPDGMass();
389  //edm::LogInfo("SimG4CoreCustomPhysics")<<"M_R: "<<m_1/GeV<<" GeV, M_np: "<<m_2/GeV<<" GeV"<<G4endl;
390  G4double sqrts = sqrt(m_1 * m_1 + m_2 * (m_2 + 2 * E_incident));
391  //edm::LogInfo("SimG4CoreCustomPhysics")<<"sqrt(s) = "<<sqrts/GeV<<" GeV"<<G4endl;
392  // Sum of rest masses after reaction:
393  G4double M_after = 0;
394  for (ReactionProduct::const_iterator r_it = aReaction.begin(); r_it != aReaction.end(); r_it++) {
395  //edm::LogInfo("SimG4CoreCustomPhysics")<<"Mass contrib: "<<(particleTable->FindParticle(*r_it)->GetPDGMass())/MeV<<" MeV"<<G4endl;
396  M_after += particleTable->FindParticle(*r_it)->GetPDGMass();
397  }
398  //edm::LogInfo("SimG4CoreCustomPhysics")<<"Intending to return this ReactionProductMass: "<<(sqrts - M_after)/MeV<<" MeV"<<G4endl;
399  return sqrts - M_after;
400 }
G4ParticleTable * particleTable
T sqrt(T t)
Definition: SSEVec.h:19
G4ParticleDefinition * theTarget
void G4ProcessHelper::ReadAndParse ( const G4String &  str,
std::vector< G4String > &  tokens,
const G4String &  delimiters = " " 
)
private

Definition at line 422 of file G4ProcessHelper.cc.

References groupFilesInBlocks::temp.

422  {
423  // Skip delimiters at beginning.
424  G4String::size_type lastPos = str.find_first_not_of(delimiters, 0);
425  // Find first "non-delimiter".
426  G4String::size_type pos = str.find_first_of(delimiters, lastPos);
427 
428  while (G4String::npos != pos || G4String::npos != lastPos) {
429  //Skipping leading / trailing whitespaces
430  G4String temp = str.substr(lastPos, pos - lastPos);
431  while (temp.c_str()[0] == ' ')
432  temp.erase(0, 1);
433  while (temp[temp.size() - 1] == ' ')
434  temp.erase(temp.size() - 1, 1);
435  // Found a token, add it to the vector.
436  tokens.push_back(temp);
437  // Skip delimiters. Note the "not_of"
438  lastPos = str.find_first_not_of(delimiters, pos);
439  // Find next "non-delimiter"
440  pos = str.find_first_of(delimiters, lastPos);
441  }
442 }
uint16_t size_type
#define str(s)
double G4ProcessHelper::Regge ( const double  boost)
private

Definition at line 444 of file G4ProcessHelper.cc.

References a, b, c, funct::exp(), and log.

444  {
445  double a = 2.165635078566177;
446  double b = 0.1467453738547229;
447  double c = -0.9607903711871166;
448  return 1.5 * exp(a + b / boost + c * log(boost));
449 }
static std::vector< std::string > checklist log
const edm::EventSetup & c
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
double b
Definition: hdecay.h:118
double a
Definition: hdecay.h:119

Member Data Documentation

double G4ProcessHelper::amplitude
private

Definition at line 78 of file G4ProcessHelper.h.

G4double G4ProcessHelper::checkfraction
private

Definition at line 46 of file G4ProcessHelper.h.

double G4ProcessHelper::ek_0
private

Definition at line 76 of file G4ProcessHelper.h.

CustomParticleFactory* G4ProcessHelper::fParticleFactory
private

Definition at line 90 of file G4ProcessHelper.h.

double G4ProcessHelper::gamma
private

Definition at line 77 of file G4ProcessHelper.h.

TProfile* G4ProcessHelper::h_q_gamma
private

Definition at line 97 of file G4ProcessHelper.h.

TProfile* G4ProcessHelper::h_q_p
private

Definition at line 96 of file G4ProcessHelper.h.

TH1D* G4ProcessHelper::h_sqrts
private

Definition at line 95 of file G4ProcessHelper.h.

TProfile* G4ProcessHelper::h_xsec_cms
private

Definition at line 94 of file G4ProcessHelper.h.

TProfile* G4ProcessHelper::h_xsec_lab
private

Definition at line 93 of file G4ProcessHelper.h.

double G4ProcessHelper::hadronlifetime
private

Definition at line 81 of file G4ProcessHelper.h.

std::map<const G4ParticleDefinition*, G4bool> G4ProcessHelper::known_particles
private

Definition at line 69 of file G4ProcessHelper.h.

double G4ProcessHelper::mixing
private

Definition at line 82 of file G4ProcessHelper.h.

G4int G4ProcessHelper::n_22
private

Definition at line 47 of file G4ProcessHelper.h.

G4int G4ProcessHelper::n_23
private

Definition at line 48 of file G4ProcessHelper.h.

ReactionMap G4ProcessHelper::nReactionMap
private

Definition at line 88 of file G4ProcessHelper.h.

std::map<G4String, G4double> G4ProcessHelper::parameters
private

Definition at line 72 of file G4ProcessHelper.h.

G4ParticleTable* G4ProcessHelper::particleTable
private

Definition at line 91 of file G4ProcessHelper.h.

ReactionMap G4ProcessHelper::pReactionMap
private

Definition at line 85 of file G4ProcessHelper.h.

bool G4ProcessHelper::reggemodel
private

Definition at line 80 of file G4ProcessHelper.h.

bool G4ProcessHelper::resonant
private

Definition at line 75 of file G4ProcessHelper.h.

double G4ProcessHelper::suppressionfactor
private

Definition at line 79 of file G4ProcessHelper.h.

HistoHelper* G4ProcessHelper::theHistoHelper
private

Definition at line 92 of file G4ProcessHelper.h.

G4ParticleDefinition* G4ProcessHelper::theNeutron
private

Definition at line 52 of file G4ProcessHelper.h.

G4ParticleDefinition* G4ProcessHelper::theProton
private

Definition at line 51 of file G4ProcessHelper.h.

G4ParticleDefinition* G4ProcessHelper::theRbaryoncloud
private

Definition at line 54 of file G4ProcessHelper.h.

ReactionMap* G4ProcessHelper::theReactionMap
private

Definition at line 56 of file G4ProcessHelper.h.

G4ParticleDefinition* G4ProcessHelper::theRmesoncloud
private

Definition at line 53 of file G4ProcessHelper.h.

G4ParticleDefinition* G4ProcessHelper::theTarget
private

Definition at line 50 of file G4ProcessHelper.h.