5 #include "G4ProcessHelper.hh"
8 #include "G4ParticleTable.hh"
9 #include "G4DecayTable.hh"
10 #include "CLHEP/Random/RandFlat.h"
11 #include "CLHEP/Units/PhysicalConstants.h"
18 G4ProcessHelper::G4ProcessHelper()
22 G4cout <<
"G4ProcessHelper constructor: start" << G4endl;
23 particleTable = G4ParticleTable::GetParticleTable();
24 theProton = particleTable->FindParticle(
"proton");
25 theNeutron = particleTable->FindParticle(
"neutron");
30 std::ifstream process_stream (
"ProcessList.txt");
32 while(getline(process_stream,
line)){
33 std::vector<G4String>
tokens;
39 G4String incident =
tokens[0];
41 G4ParticleDefinition* incidentDef = particleTable->FindParticle(incident);
42 G4int incidentPDG = incidentDef->GetPDGEncoding();
43 known_particles[incidentDef]=
true;
50 for (
unsigned int i = 2;
i !=
tokens.size();
i++){
52 if (particleTable->contains(
part))
54 prod.push_back(particleTable->FindParticle(
part)->GetPDGEncoding());
56 G4cout<<
"Particle: "<<
part<<
" is unknown."<<G4endl;
57 G4Exception(
"G4ProcessHelper",
"UnkownParticle", FatalException,
58 "Initialization: The reaction product list contained an unknown particle");
62 pReactionMap[incidentPDG].push_back(prod);
63 }
else if (
target ==
"neutron") {
64 nReactionMap[incidentPDG].push_back(prod);
66 G4Exception(
"G4ProcessHelper",
"IllegalTarget", FatalException,
67 "Initialization: The reaction product list contained an illegal target particle");
71 process_stream.close();
72 G4cout <<
"Found " << pReactionMap.size() <<
" proton interactions and " << nReactionMap.size() <<
" neutron interactions in ProcessList.txt." << G4endl;
84 suppressionfactor =
parameters[
"ReggeSuppression"];
87 if(
parameters[
"ReggeModel"]!=0.) reggemodel=
true;
90 G4cout<<
"Read in physics parameters:"<<G4endl;
91 G4cout<<
"Resonant = "<< resonant <<G4endl;
92 G4cout<<
"ResonanceEnergy = "<<ek_0/
CLHEP::GeV<<
" GeV"<<G4endl;
93 G4cout<<
"XsecMultiplier = "<<xsecmultiplier<<G4endl;
96 G4cout<<
"ReggeSuppression = "<<100*suppressionfactor<<
" %"<<G4endl;
97 G4cout<<
"HadronLifeTime = "<<hadronlifetime;
98 if (doDecays) G4cout<<
" ns"<<G4endl;
99 else G4cout<<
" s"<<G4endl;
100 G4cout<<
"ReggeModel = "<< reggemodel <<G4endl;
101 G4cout<<
"Mixing = "<< mixing*100 <<
" %"<<G4endl;
102 G4cout<<
"DoDecays = "<< doDecays << G4endl;
104 if ((!doDecays && hadronlifetime>0.) ||
105 (doDecays && hadronlifetime<=0.) ){
106 G4cout <<
"WARNING: Inconsistent treatment of R-Hadron properties! Lifetime of " << hadronlifetime
107 <<
" and doDecays= " << doDecays << G4endl;
110 G4ParticleTable::G4PTblDicIterator* theParticleIterator;
111 theParticleIterator = particleTable->GetIterator();
113 theParticleIterator->reset();
114 while( (*theParticleIterator)() ){
116 std::string
name = theParticleIterator->value()->GetParticleName();
117 G4DecayTable*
table = theParticleIterator->value()->GetDecayTable();
122 G4cout<<
"Stable: "<<
particle->GetPDGStable()<<G4endl;
128 G4cout<<
"Forcing a decay for "<<
name<<G4endl;
130 G4cout<<
"Stable: "<<
particle->GetPDGStable()<<G4endl;
133 G4cout <<
"Done with particle " <<
name << G4endl;
135 theParticleIterator->reset();
136 G4cout <<
"G4ProcessHelper constructor: end" << G4endl;
141 void G4ProcessHelper::ReadInPhysicsParameters(std::map<G4String,G4double>&
parameters)
const
154 std::ifstream physics_stream (
"PhysicsConfiguration.txt");
157 while (getline(physics_stream,
line)) {
158 std::vector<G4String>
tokens;
165 physics_stream.close();
169 const G4ProcessHelper* G4ProcessHelper::Instance()
171 static const G4ProcessHelper
instance;
176 G4bool G4ProcessHelper::ApplicabilityTester(
const G4ParticleDefinition& aPart)
const {
178 return known_particles.at(&aPart);
180 catch (
const std::out_of_range&
e) {
185 G4double G4ProcessHelper::GetInclusiveCrossSection(
const G4DynamicParticle *aParticle,
186 const G4Element *anElement)
const {
190 G4int thePDGCode = aParticle->GetDefinition()->GetPDGEncoding();
191 double boost = (aParticle->GetKineticEnergy()+aParticle->GetMass())/aParticle->GetMass();
192 G4double theXsec = 0;
193 G4String
name = aParticle->GetDefinition()->GetParticleName();
216 double P = Pom(
boost);
241 double e_0 = ek_0 + aParticle->GetDefinition()->GetPDGMass();
243 e_0 = sqrt(aParticle->GetDefinition()->GetPDGMass()*aParticle->GetDefinition()->GetPDGMass()
244 + theProton->GetPDGMass()*theProton->GetPDGMass()
245 + 2.*e_0*theProton->GetPDGMass());
246 double sqrts=sqrt(aParticle->GetDefinition()->GetPDGMass()*aParticle->GetDefinition()->GetPDGMass()
247 + theProton->GetPDGMass()*theProton->GetPDGMass() + 2*aParticle->GetTotalEnergy()*theProton->GetPDGMass());
251 theXsec += res_result;
255 return theXsec *
pow(anElement->GetN(),0.7)*1.25 * xsecmultiplier;
259 ReactionProduct G4ProcessHelper::GetFinalState(
const G4Track& aTrack, G4ParticleDefinition*& aTarget)
const {
260 return GetFinalStateInternal(aTrack,aTarget,
false);
264 ReactionProduct G4ProcessHelper::GetFinalStateInternal(
const G4Track& aTrack,G4ParticleDefinition*& aTarget,
const bool baryonize_failed)
const {
266 const G4DynamicParticle* aDynamicParticle = aTrack.GetDynamicParticle();
273 G4Material* aMaterial = aTrack.GetMaterial();
274 const G4ElementVector* theElementVector = aMaterial->GetElementVector() ;
275 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
277 G4double NumberOfProtons=0;
278 G4double NumberOfNucleons=0;
280 for (
size_t elm=0 ; elm < aMaterial->GetNumberOfElements() ; elm++ )
283 NumberOfProtons += NbOfAtomsPerVolume[elm]*(*theElementVector)[elm]->GetZ();
285 NumberOfNucleons += NbOfAtomsPerVolume[elm]*(*theElementVector)[elm]->GetN();
288 const ReactionMap* reactionMap;
289 if(CLHEP::RandFlat::shoot()<NumberOfProtons/NumberOfNucleons){
290 reactionMap = &pReactionMap;
293 reactionMap = &nReactionMap;
294 aTarget = theNeutron;
297 G4int theIncidentPDG = aDynamicParticle->GetDefinition()->GetPDGEncoding();
301 &&CLHEP::RandFlat::shoot()*mixing>0.5
302 &&aDynamicParticle->GetDefinition()->GetPDGCharge()==0.
306 theIncidentPDG *= -1;
310 bool baryonise=
false;
314 && CLHEP::RandFlat::shoot()>0.9
325 const ReactionProductList& aReactionProductList = reactionMap->at(theIncidentPDG);
337 ReactionProductList theReactionProductList;
338 std::vector<bool> theChargeChangeList;
340 for (
const ReactionProduct& prod : aReactionProductList) {
341 G4int secondaries = prod.size();
347 if(ReactionIsPossible(prod,*aTarget,aDynamicParticle)
349 (baryonise&&ReactionGivesBaryon(prod))
351 (!baryonise&&!ReactionGivesBaryon(prod))
361 theReactionProductList.push_back(prod);
362 if (secondaries == 2){
364 }
else if (secondaries ==3) {
367 G4cerr <<
"ReactionProduct has unsupported number of secondaries: "<<secondaries<<G4endl;
372 if (theReactionProductList.size()==0 && baryonize_failed){
373 G4Exception(
"G4ProcessHelper",
"NoProcessPossible", FatalException,
374 "GetFinalState: No process could be selected from the given list.");
375 }
else if (theReactionProductList.size()==0 && !baryonize_failed) {
377 G4cout <<
"G4ProcessHelper::GetFinalStateInternal WARNING Could not select an appropriate process in first pass" << G4endl;
378 return GetFinalStateInternal(aTrack,aTarget,
true);
384 int n_rps = theReactionProductList.size();
385 int select = (
int)(CLHEP::RandFlat::shoot()*n_rps);
387 return theReactionProductList[
select];
393 G4double p23 = 1-p22;
395 std::vector<G4double> Probabilities;
396 std::vector<G4bool> TwotoThreeFlag;
398 G4double CumulatedProbability = 0;
402 for (
unsigned int i = 0;
i != theReactionProductList.size();
i++){
403 if (theReactionProductList[
i].
size() == 2) {
404 if(0==N22) {
throw std::runtime_error(
"G4ProcessHelper::GetFinalState: N22 is zero!"); }
405 CumulatedProbability += p22/N22;
406 TwotoThreeFlag.push_back(
false);
408 if(0==N23) {
throw std::runtime_error(
"G4ProcessHelper::GetFinalState: N23 is zero!"); }
409 CumulatedProbability += p23/N23;
410 TwotoThreeFlag.push_back(
true);
412 Probabilities.push_back(CumulatedProbability);
417 it != Probabilities.end();
420 *
it /= CumulatedProbability;
425 G4bool selected =
false;
431 while(!selected && tries < 100){
433 G4double dice = CLHEP::RandFlat::shoot();
434 while(
i<theReactionProductList.size() && dice>Probabilities[
i] ){
438 if(!TwotoThreeFlag[
i]) {
443 if (PhaseSpace(theReactionProductList[
i],*aTarget,aDynamicParticle)>CLHEP::RandFlat::shoot()) selected =
true;
447 if(selected &&
table->FindParticle(theReactionProductList[
i][0])->GetPDGCharge()!=aDynamicParticle->GetDefinition()->GetPDGCharge())
455 if(CLHEP::RandFlat::shoot()<suppressionfactor) selected =
false;
460 if(tries>=100) G4cerr<<
"Could not select process!!!!"<<G4endl;
463 return theReactionProductList[
i];
466 G4double G4ProcessHelper::ReactionProductMass(
const ReactionProduct& aReaction,
const G4ParticleDefinition& aTarget,
const G4DynamicParticle* aDynamicParticle)
const{
468 G4double E_incident = aDynamicParticle->GetTotalEnergy();
471 G4double m_1 = aDynamicParticle->GetDefinition()->GetPDGMass();
472 G4double m_2 = aTarget.GetPDGMass();
474 G4double sqrts = sqrt(m_1*m_1 + m_2*(m_2 + 2 * E_incident));
477 G4double M_after = 0;
478 for (ReactionProduct::const_iterator r_it = aReaction.begin(); r_it !=aReaction.end(); ++r_it){
481 M_after +=
table->FindParticle(*r_it)->GetPDGMass();
484 return sqrts - M_after;
487 G4bool G4ProcessHelper::ReactionIsPossible(
const ReactionProduct& aReaction,
const G4ParticleDefinition& aTarget,
const G4DynamicParticle* aDynamicParticle)
const{
488 if (ReactionProductMass(aReaction,aTarget,aDynamicParticle)>0)
return true;
492 G4bool G4ProcessHelper::ReactionGivesBaryon(
const ReactionProduct& aReaction)
const{
493 for (ReactionProduct::const_iterator
it = aReaction.begin();
it!=aReaction.end();++
it)
498 G4double G4ProcessHelper::PhaseSpace(
const ReactionProduct& aReaction,
const G4ParticleDefinition& aTarget,
const G4DynamicParticle* aDynamicParticle)
const {
499 G4double qValue = ReactionProductMass(aReaction,aTarget,aDynamicParticle);
505 void G4ProcessHelper::ReadAndParse(
const G4String&
str,
506 std::vector<G4String>&
tokens,
507 const G4String& delimiters)
const
510 G4String::size_type lastPos =
str.find_first_not_of(delimiters, 0);
511 if (lastPos==G4String::npos)
return;
514 G4String::size_type
pos =
str.find_first_of(delimiters, lastPos);
516 while (G4String::npos !=
pos || G4String::npos != lastPos)
519 G4String temp =
str.substr(lastPos,
pos - lastPos);
520 while(temp.c_str()[0] ==
' ') temp.erase(0,1);
521 while(temp[temp.size()-1] ==
' ') temp.erase(temp.size()-1,1);
527 lastPos =
str.find_first_not_of(delimiters,
pos);
528 if (lastPos==G4String::npos)
continue;
531 pos =
str.find_first_of(delimiters, lastPos);
535 double G4ProcessHelper::Regge(
const double boost)
const
539 double a=2.165635078566177;
540 double b=0.1467453738547229;
541 double c=-0.9607903711871166;
546 double G4ProcessHelper::Pom(
const double boost)
const
550 double a=4.138224000651535;
551 double b=1.50377557581421;
552 double c=-0.05449742257808247;
553 double d=0.0008221235048211401;