1 #include"G4ParticleTable.hh"
2 #include "Randomize.hh"
12 #include"SimG4Core/CustomPhysics/interface/G4ProcessHelper.hh"
18 particleTable = G4ParticleTable::GetParticleTable();
20 theProton = particleTable->FindParticle(
"proton");
21 theNeutron = particleTable->FindParticle(
"neutron");
26 std::string processDefFilePath = fp.
fullPath();
27 std::ifstream process_stream (processDefFilePath.c_str());
32 amplitude = p.
getParameter<
double>(
"amplitude")*millibarn;
33 suppressionfactor = p.
getParameter<
double>(
"reggeSuppression");
34 hadronlifetime = p.
getParameter<
double>(
"hadronLifeTime");
38 edm::LogInfo(
"CustomPhysics")<<
"Read in physics parameters:"<<G4endl;
39 edm::LogInfo(
"CustomPhysics")<<
"Resonant = "<< resonant <<G4endl;
40 edm::LogInfo(
"CustomPhysics")<<
"ResonanceEnergy = "<<ek_0/GeV<<
" GeV"<<G4endl;
41 edm::LogInfo(
"CustomPhysics")<<
"Gamma = "<<gamma/GeV<<
" GeV"<<G4endl;
42 edm::LogInfo(
"CustomPhysics")<<
"Amplitude = "<<amplitude/millibarn<<
" millibarn"<<G4endl;
43 edm::LogInfo(
"CustomPhysics")<<
"ReggeSuppression = "<<100*suppressionfactor<<
" %"<<G4endl;
44 edm::LogInfo(
"CustomPhysics")<<
"HadronLifeTime = "<<hadronlifetime<<
" s"<<G4endl;
45 edm::LogInfo(
"CustomPhysics")<<
"ReggeModel = "<< reggemodel <<G4endl;
46 edm::LogInfo(
"CustomPhysics")<<
"Mixing = "<< mixing*100 <<
" %"<<G4endl;
55 while(getline(process_stream,line)){
56 std::vector<G4String> tokens;
58 ReadAndParse(line,tokens,
"#");
60 G4String incident = tokens[0];
62 G4ParticleDefinition* incidentDef = particleTable->FindParticle(incident);
64 G4int incidentPDG = incidentDef->GetPDGEncoding();
65 known_particles[incidentDef]=
true;
67 G4String
target = tokens[1];
69 <<
" Target: "<<target<<G4endl;
73 for (
size_t i = 2;
i != tokens.size();
i++){
74 G4String
part = tokens[
i];
75 if (particleTable->contains(part))
77 prod.push_back(particleTable->FindParticle(part)->GetPDGEncoding());
79 edm::LogInfo(
"CustomPhysics")<<
"Particle: "<<part<<
" is unknown."<<G4endl;
80 G4Exception(
"G4ProcessHelper",
"UnkownParticle", FatalException,
81 "Initialization: The reaction product list contained an unknown particle");
84 if (target ==
"proton")
86 pReactionMap[incidentPDG].push_back(prod);
87 }
else if (target ==
"neutron") {
88 nReactionMap[incidentPDG].push_back(prod);
90 G4Exception(
"G4ProcessHelper",
"IllegalTarget", FatalException,
91 "Initialization: The reaction product list contained an illegal target particle");
95 process_stream.close();
97 G4ParticleTable::G4PTblDicIterator* theParticleIterator;
98 theParticleIterator = particleTable->GetIterator();
100 theParticleIterator->reset();
101 while( (*theParticleIterator)() ){
103 std::string
name = theParticleIterator->value()->GetParticleName();
104 G4DecayTable*
table = theParticleIterator->value()->GetDecayTable();
105 if(particle!=0&&table!=0&&name.find(
"cloud")>name.size()&&hadronlifetime > 0)
107 particle->SetPDGLifeTime(hadronlifetime*
s);
108 particle->SetPDGStable(
false);
109 edm::LogInfo(
"CustomPhysics")<<
"Lifetime of: "<<name<<
" set to: "<<particle->GetPDGLifeTime()/s<<
" s."<<G4endl;
110 edm::LogInfo(
"CustomPhysics")<<
"Stable: "<<particle->GetPDGStable()<<G4endl;
113 theParticleIterator->reset();
118 G4bool G4ProcessHelper::ApplicabilityTester(
const G4ParticleDefinition& aPart){
119 const G4ParticleDefinition* aP = &aPart;
120 if (known_particles[aP])
return true;
124 G4double G4ProcessHelper::GetInclusiveCrossSection(
const G4DynamicParticle *aParticle,
125 const G4Element *anElement){
132 G4int thePDGCode = aParticle->GetDefinition()->GetPDGEncoding();
133 double boost = (aParticle->GetKineticEnergy()+aParticle->GetMass())/aParticle->GetMass();
135 G4double theXsec = 0;
136 G4String name = aParticle->GetDefinition()->GetParticleName();
141 theXsec = 24 * millibarn;
145 for (std::vector<G4int>::iterator it = nq.begin();
150 if (*it == 1 || *it == 2) theXsec += 12 * millibarn;
151 if (*it == 3) theXsec += 6 * millibarn;
157 double R = Regge(boost);
158 double P = Pom(boost);
181 double e_0 = ek_0 + aParticle->GetDefinition()->GetPDGMass();
183 e_0 =
sqrt(aParticle->GetDefinition()->GetPDGMass()*aParticle->GetDefinition()->GetPDGMass()
184 + theProton->GetPDGMass()*theProton->GetPDGMass()
185 + 2.*e_0*theProton->GetPDGMass());
188 double sqrts=
sqrt(aParticle->GetDefinition()->GetPDGMass()*aParticle->GetDefinition()->GetPDGMass()
189 + theProton->GetPDGMass()*theProton->GetPDGMass() + 2*aParticle->GetTotalEnergy()*theProton->GetPDGMass());
191 double res_result = amplitude*(gamma*gamma/4.)/((sqrts-e_0)*(sqrts-e_0)+(gamma*gamma/4.));
193 theXsec += res_result;
200 return theXsec *
pow(anElement->GetN(),0.7)*1.25;
204 ReactionProduct G4ProcessHelper::GetFinalState(
const G4Track& aTrack, G4ParticleDefinition*& aTarget){
206 const G4DynamicParticle* aDynamicParticle = aTrack.GetDynamicParticle();
213 G4Material* aMaterial = aTrack.GetMaterial();
214 const G4ElementVector* theElementVector = aMaterial->GetElementVector() ;
215 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
217 G4double NumberOfProtons=0;
218 G4double NumberOfNucleons=0;
220 for (
size_t elm=0 ;
elm < aMaterial->GetNumberOfElements() ;
elm++ )
223 NumberOfProtons += NbOfAtomsPerVolume[
elm]*(*theElementVector)[
elm]->GetZ();
225 NumberOfNucleons += NbOfAtomsPerVolume[
elm]*(*theElementVector)[
elm]->GetN();
228 if(G4UniformRand()<NumberOfProtons/NumberOfNucleons)
230 theReactionMap = &pReactionMap;
231 theTarget = theProton;
233 theReactionMap = &nReactionMap;
234 theTarget = theNeutron;
238 G4int theIncidentPDG = aDynamicParticle->GetDefinition()->GetPDGEncoding();
242 && G4UniformRand()*mixing>0.5
243 &&aDynamicParticle->GetDefinition()->GetPDGCharge()==0.
247 theIncidentPDG *= -1;
251 bool baryonise=
false;
254 && G4UniformRand()>0.9
264 ReactionProductList* aReactionProductList = &((*theReactionMap)[theIncidentPDG]);
276 ReactionProductList theReactionProductList;
277 std::vector<bool> theChargeChangeList;
279 for (ReactionProductList::iterator prod_it = aReactionProductList->begin();
280 prod_it != aReactionProductList->end();
282 G4int secondaries = prod_it->size();
285 if(ReactionIsPossible(*prod_it,aDynamicParticle)
287 (baryonise&&ReactionGivesBaryon(*prod_it))
289 (!baryonise&&!ReactionGivesBaryon(*prod_it))
298 theReactionProductList.push_back(*prod_it);
299 if (secondaries == 2){
301 }
else if (secondaries ==3) {
304 G4cerr <<
"ReactionProduct has unsupported number of secondaries: "<<secondaries<<G4endl;
312 if (theReactionProductList.size()==0) G4Exception(
"G4ProcessHelper",
"NoProcessPossible", FatalException,
313 "GetFinalState: No process could be selected from the given list.");
318 int n_rps = theReactionProductList.size();
319 int select = (int)( G4UniformRand()*n_rps);
321 return theReactionProductList[
select];
327 G4double p23 = 1-p22;
329 std::vector<G4double> Probabilities;
330 std::vector<G4bool> TwotoThreeFlag;
332 G4double CumulatedProbability = 0;
336 for (
unsigned int i = 0;
i != theReactionProductList.size();
i++){
337 if (theReactionProductList[
i].
size() == 2) {
338 CumulatedProbability += p22/N22;
339 TwotoThreeFlag.push_back(
false);
341 CumulatedProbability += p23/N23;
342 TwotoThreeFlag.push_back(
true);
344 Probabilities.push_back(CumulatedProbability);
350 for (std::vector<G4double>::iterator it = Probabilities.begin();
351 it != Probabilities.end();
354 *it /= CumulatedProbability;
361 G4bool selected =
false;
367 while(!selected && tries < 100){
369 G4double dice = G4UniformRand();
371 while(dice>Probabilities[i] && i<theReactionProductList.size()){
378 if(!TwotoThreeFlag[i]) {
383 if (PhaseSpace(theReactionProductList[i],aDynamicParticle)>G4UniformRand()) selected =
true;
387 if(selected&&particleTable->FindParticle(theReactionProductList[i][0])->GetPDGCharge()!=aDynamicParticle->GetDefinition()->GetPDGCharge())
395 if(G4UniformRand()<suppressionfactor) selected =
false;
400 if(tries>=100) G4cerr<<
"Could not select process!!!!"<<G4endl;
406 if (theReactionProductList[i].
size()==2) {
412 checkfraction = (1.0*n_22)/(n_22+n_23);
416 return theReactionProductList[
i];
419 G4double G4ProcessHelper::ReactionProductMass(
const ReactionProduct& aReaction,
const G4DynamicParticle* aDynamicParticle){
421 G4double E_incident = aDynamicParticle->GetTotalEnergy();
424 G4double m_1 = aDynamicParticle->GetDefinition()->GetPDGMass();
425 G4double m_2 = theTarget->GetPDGMass();
427 G4double sqrts =
sqrt(m_1*m_1 + m_2*(m_2 + 2 * E_incident));
430 G4double M_after = 0;
431 for (ReactionProduct::const_iterator r_it = aReaction.begin(); r_it !=aReaction.end(); r_it++){
433 M_after += particleTable->FindParticle(*r_it)->GetPDGMass();
436 return sqrts - M_after;
439 G4bool G4ProcessHelper::ReactionIsPossible(
const ReactionProduct& aReaction,
const G4DynamicParticle* aDynamicParticle){
440 if (ReactionProductMass(aReaction,aDynamicParticle)>0)
return true;
444 G4bool G4ProcessHelper::ReactionGivesBaryon(
const ReactionProduct& aReaction){
445 for (ReactionProduct::const_iterator it = aReaction.begin();it!=aReaction.end();it++)
450 G4double G4ProcessHelper::PhaseSpace(
const ReactionProduct& aReaction,
const G4DynamicParticle* aDynamicParticle){
451 G4double qValue = ReactionProductMass(aReaction,aDynamicParticle);
452 G4double
phi =
sqrt(1+qValue/(2*0.139*GeV))*
pow(qValue/(1.1*GeV),3./2.);
453 return (phi/(1+phi));
456 void G4ProcessHelper::ReadAndParse(
const G4String& str,
457 std::vector<G4String>& tokens,
458 const G4String& delimiters)
465 while (G4String::npos != pos || G4String::npos != lastPos)
468 G4String
temp = str.substr(lastPos, pos - lastPos);
469 while(temp.c_str()[0] ==
' ') temp.erase(0,1);
470 while(temp[temp.size()-1] ==
' ') temp.erase(temp.size()-1,1);
472 tokens.push_back(temp);
474 lastPos = str.find_first_not_of(delimiters, pos);
476 pos = str.find_first_of(delimiters, lastPos);
480 double G4ProcessHelper::Regge(
const double boost)
482 double a=2.165635078566177;
483 double b=0.1467453738547229;
484 double c=-0.9607903711871166;
485 return 1.5*
exp(a+b/boost+c*
log(boost));
489 double G4ProcessHelper::Pom(
const double boost)
491 double a=4.138224000651535;
492 double b=1.50377557581421;
493 double c=-0.05449742257808247;
494 double d=0.0008221235048211401;
495 return a + b*
sqrt(boost) + c*boost + d*
pow(boost,1.5);
T getParameter(std::string const &) const
static bool s_isRGlueball(int pdg)
static bool s_isSbaryon(int pdg)
static bool s_isRMeson(int pdg)
static bool s_isMesonino(int pdg)
static bool s_isRBaryon(int pdg)
static std::vector< int > s_containedQuarks(int pdg)
std::string fullPath() const
tuple size
Write out results.
Power< A, B >::type pow(const A &a, const B &b)