#include <G4HEProtonInelastic.hh>
Inheritance diagram for G4HEProtonInelastic:
Public Member Functions | |
G4HEProtonInelastic () | |
~G4HEProtonInelastic () | |
virtual void | ModelDescription (std::ostream &) const |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
G4int | GetNumberOfSecondaries () |
void | FirstIntInCasProton (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight) |
Data Fields | |
G4int | vecLength |
Definition at line 52 of file G4HEProtonInelastic.hh.
G4HEProtonInelastic::G4HEProtonInelastic | ( | ) | [inline] |
Definition at line 55 of file G4HEProtonInelastic.hh.
References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, vecLength, and G4HEInelastic::verboseLevel.
00055 : G4HEInelastic("G4HEProtonInelastic") 00056 { 00057 vecLength = 0; 00058 theMinEnergy = 45*CLHEP::GeV; 00059 theMaxEnergy = 10*CLHEP::TeV; 00060 MAXPART = 2048; 00061 verboseLevel = 0; 00062 G4cout << "WARNING: model G4HEProtonInelastic is being deprecated and will\n" 00063 << "disappear in Geant4 version 10.0" << G4endl; 00064 }
G4HEProtonInelastic::~G4HEProtonInelastic | ( | ) | [inline] |
G4HadFinalState * G4HEProtonInelastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | targetNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 57 of file G4HEProtonInelastic.cc.
References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasProton(), G4cout, G4endl, G4UniformRand, G4Nucleus::GetA_asInt(), G4HEVector::getCode(), G4HEVector::getEnergy(), G4HEVector::getKineticEnergy(), G4HEVector::getMass(), G4HEVector::getName(), G4Nucleus::GetZ_asInt(), G4HEInelastic::HighEnergyCascading(), G4HEInelastic::HighEnergyClusterProduction(), G4HEInelastic::MAXPART, G4HEInelastic::MediumEnergyCascading(), G4HEInelastic::MediumEnergyClusterProduction(), G4HEInelastic::NuclearExcitation(), G4HEInelastic::NuclearInelasticity(), G4HEInelastic::QuasiElasticScattering(), G4HEVector::setDefinition(), G4HadFinalState::SetStatusChange(), stopAndKill, G4HEInelastic::StrangeParticlePairProduction(), G4HadronicInteraction::theParticleChange, vecLength, and G4HEInelastic::verboseLevel.
00059 { 00060 G4HEVector* pv = new G4HEVector[MAXPART]; 00061 const G4HadProjectile* aParticle = &aTrack; 00062 const G4double A = targetNucleus.GetA_asInt(); 00063 const G4double Z = targetNucleus.GetZ_asInt(); 00064 G4HEVector incidentParticle(aParticle); 00065 00066 G4double atomicNumber = Z; 00067 G4double atomicWeight = A; 00068 00069 if (verboseLevel > 1) 00070 G4cout << "Z , A = " << atomicNumber << " " << atomicWeight << G4endl; 00071 00072 G4int incidentCode = incidentParticle.getCode(); 00073 G4double incidentMass = incidentParticle.getMass(); 00074 G4double incidentTotalEnergy = incidentParticle.getEnergy(); 00075 00076 // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00077 // DHW 19 May 2011: variable set but not used 00078 00079 G4double incidentKineticEnergy = incidentParticle.getKineticEnergy(); 00080 00081 if (incidentKineticEnergy < 1.) 00082 G4cout << "GHEProtonInelastic: incident energy < 1 GeV" << G4endl; 00083 00084 if (verboseLevel > 1) { 00085 G4cout << "G4HEProtonInelastic::ApplyYourself" << G4endl; 00086 G4cout << "incident particle " << incidentParticle.getName() << " " 00087 << "mass " << incidentMass << " " 00088 << "kinetic energy " << incidentKineticEnergy 00089 << G4endl; 00090 G4cout << "target material with (A,Z) = (" 00091 << atomicWeight << "," << atomicNumber << ")" << G4endl; 00092 } 00093 00094 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 00095 atomicWeight, atomicNumber); 00096 if (verboseLevel > 1) 00097 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 00098 00099 incidentKineticEnergy -= inelasticity; 00100 00101 G4double excitationEnergyGNP = 0.; 00102 G4double excitationEnergyDTA = 0.; 00103 00104 G4double excitation = NuclearExcitation(incidentKineticEnergy, 00105 atomicWeight, atomicNumber, 00106 excitationEnergyGNP, 00107 excitationEnergyDTA); 00108 00109 if (verboseLevel > 1) 00110 G4cout << "nuclear excitation = " << excitation << " " 00111 << excitationEnergyGNP << " " << excitationEnergyDTA << G4endl; 00112 00113 incidentKineticEnergy -= excitation; 00114 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 00115 // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 00116 // *(incidentTotalEnergy+incidentMass)); 00117 // DHW 19 May 2011: variable set but not used 00118 00119 G4HEVector targetParticle; 00120 if (G4UniformRand() < atomicNumber/atomicWeight) { 00121 targetParticle.setDefinition("Proton"); 00122 } else { 00123 targetParticle.setDefinition("Neutron"); 00124 } 00125 00126 G4double targetMass = targetParticle.getMass(); 00127 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass 00128 + targetMass*targetMass 00129 + 2.0*targetMass*incidentTotalEnergy); 00130 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 00131 00132 // In the original Gheisha code, the inElastic flag was defined as follows: 00133 // G4bool inElastic = InElasticCrossSectionInFirstInt 00134 // (availableEnergy, incidentCode, incidentTotalMomentum); 00135 // For unknown reasons, it was replaced by the following code in Geant??? 00136 00137 G4bool inElastic = true; 00138 // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false; 00139 00140 vecLength = 0; 00141 00142 if (verboseLevel > 1) 00143 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 00144 << incidentCode << G4endl; 00145 00146 G4bool successful = false; 00147 00148 FirstIntInCasProton(inElastic, availableEnergy, pv, vecLength, 00149 incidentParticle, targetParticle, atomicWeight); 00150 00151 if (verboseLevel > 1) 00152 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 00153 00154 if ((vecLength > 0) && (availableEnergy > 1.)) 00155 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 00156 pv, vecLength, 00157 incidentParticle, targetParticle); 00158 00159 HighEnergyCascading(successful, pv, vecLength, 00160 excitationEnergyGNP, excitationEnergyDTA, 00161 incidentParticle, targetParticle, 00162 atomicWeight, atomicNumber); 00163 if (!successful) 00164 HighEnergyClusterProduction(successful, pv, vecLength, 00165 excitationEnergyGNP, excitationEnergyDTA, 00166 incidentParticle, targetParticle, 00167 atomicWeight, atomicNumber); 00168 if (!successful) 00169 MediumEnergyCascading(successful, pv, vecLength, 00170 excitationEnergyGNP, excitationEnergyDTA, 00171 incidentParticle, targetParticle, 00172 atomicWeight, atomicNumber); 00173 00174 if (!successful) 00175 MediumEnergyClusterProduction(successful, pv, vecLength, 00176 excitationEnergyGNP, excitationEnergyDTA, 00177 incidentParticle, targetParticle, 00178 atomicWeight, atomicNumber); 00179 if (!successful) 00180 QuasiElasticScattering(successful, pv, vecLength, 00181 excitationEnergyGNP, excitationEnergyDTA, 00182 incidentParticle, targetParticle, 00183 atomicWeight, atomicNumber); 00184 if (!successful) 00185 ElasticScattering(successful, pv, vecLength, 00186 incidentParticle, 00187 atomicWeight, atomicNumber); 00188 00189 if (!successful) 00190 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 00191 << G4endl; 00192 00193 FillParticleChange(pv, vecLength); 00194 delete [] pv; 00195 theParticleChange.SetStatusChange(stopAndKill); 00196 return &theParticleChange; 00197 }
void G4HEProtonInelastic::FirstIntInCasProton | ( | G4bool & | inElastic, | |
const G4double | availableEnergy, | |||
G4HEVector | pv[], | |||
G4int & | vecLen, | |||
const G4HEVector & | incidentParticle, | |||
const G4HEVector & | targetParticle, | |||
const G4double | atomicWeight | |||
) |
Definition at line 201 of file G4HEProtonInelastic.cc.
References G4HEInelastic::Amax(), G4HEInelastic::Amin(), G4HEInelastic::Factorial(), G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), G4HEInelastic::Imax(), CLHEP::detail::n, G4HEInelastic::Neutron, neutronCode, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, sqr(), and G4HEInelastic::verboseLevel.
Referenced by ApplyYourself().
00216 { 00217 static const G4double expxu = 82.; // upper bound for arg. of exp 00218 static const G4double expxl = -expxu; // lower bound for arg. of exp 00219 00220 static const G4double protb = 0.7; 00221 static const G4double neutb = 0.35; 00222 static const G4double c = 1.25; 00223 00224 static const G4int numMul = 1200; 00225 static const G4int numSec = 60; 00226 00227 G4int neutronCode = Neutron.getCode(); 00228 G4int protonCode = Proton.getCode(); 00229 G4double pionMass = PionPlus.getMass(); 00230 00231 G4int targetCode = targetParticle.getCode(); 00232 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00233 00234 static G4bool first = true; 00235 static G4double protmul[numMul], protnorm[numSec]; // proton constants 00236 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 00237 00238 // misc. local variables 00239 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0 00240 00241 G4int i, counter, nt, npos, nneg, nzero; 00242 00243 if (first) { 00244 // compute normalization constants, this will only be done once 00245 first = false; 00246 for (i=0; i<numMul; i++) protmul[i] = 0.0; 00247 for (i=0; i<numSec; i++) protnorm[i] = 0.0; 00248 counter = -1; 00249 for (npos=0; npos<(numSec/3); npos++) { 00250 for (nneg=Imax(0,npos-2); nneg<=npos; nneg++) { 00251 for (nzero=0; nzero<numSec/3; nzero++) { 00252 if (++counter < numMul) { 00253 nt = npos+nneg+nzero; 00254 if ( (nt>0) && (nt<=numSec) ) { 00255 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c) 00256 /(Factorial(2-npos+nneg)*Factorial(npos-nneg)) ; 00257 protnorm[nt-1] += protmul[counter]; 00258 } 00259 } 00260 } 00261 } 00262 } 00263 00264 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 00265 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 00266 counter = -1; 00267 for (npos=0; npos<numSec/3; npos++) { 00268 for (nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++) { 00269 for (nzero=0; nzero<numSec/3; nzero++) { 00270 if (++counter < numMul) { 00271 nt = npos+nneg+nzero; 00272 if ( (nt>0) && (nt<=numSec) ) { 00273 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c) 00274 /(Factorial(1-npos+nneg)*Factorial(1+npos-nneg)); 00275 neutnorm[nt-1] += neutmul[counter]; 00276 } 00277 } 00278 } 00279 } 00280 } 00281 for (i=0; i<numSec; i++) { 00282 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 00283 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 00284 } 00285 } // end of initialization 00286 00287 00288 // initialize the first two places 00289 // the same as beam and target 00290 pv[0] = incidentParticle; 00291 pv[1] = targetParticle; 00292 vecLen = 2; 00293 00294 if (!inElastic) { // quasi-elastic scattering, no pions produced 00295 if (targetCode == neutronCode) { 00296 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.}; 00297 G4int iplab = G4int( Amin( 9.0, incidentTotalMomentum*2.5 ) ); 00298 if (G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) { 00299 // charge exchange pi+ n -> pi0 p 00300 pv[0] = PionZero; 00301 pv[1] = Proton; 00302 } 00303 } 00304 return; 00305 } else if (availableEnergy <= pionMass) return; 00306 00307 // inelastic scattering 00308 npos = 0, nneg = 0, nzero = 0; 00309 G4double eab = availableEnergy; 00310 G4int ieab = G4int( eab*5.0 ); 00311 00312 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98}; 00313 if ( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) ) { 00314 // suppress high multiplicity events at low momentum 00315 // only one additional pion will be produced 00316 G4double w0, wp, wm, wt, ran; 00317 if (targetCode == protonCode) { // target is a proton 00318 w0 = - sqr(1.+protb)/(2.*c*c); 00319 wp = w0 = std::exp(w0); 00320 if (G4UniformRand() < w0/(w0+wp) ) { 00321 npos = 0; 00322 nneg = 0; 00323 nzero = 1; 00324 } else { 00325 npos = 1; 00326 nneg = 0; 00327 nzero = 0; } 00328 } else { // target is a neutron 00329 w0 = -sqr(1.+neutb)/(2.*c*c); 00330 w0 = std::exp(w0); 00331 wp = w0/2.; 00332 wm = -sqr(-1.+neutb)/(2.*c*c); 00333 wm = std::exp(wm)/2.; 00334 wt = w0+wp+wm; 00335 wp = w0+wp; 00336 ran = G4UniformRand(); 00337 if( ran < w0/wt) 00338 { npos = 0; nneg = 0; nzero = 1; } 00339 else if( ran < wp/wt) 00340 { npos = 1; nneg = 0; nzero = 0; } 00341 else 00342 { npos = 0; nneg = 1; nzero = 0; } 00343 } 00344 } else { 00345 // number of total particles vs. centre of mass Energy - 2*proton mass 00346 00347 G4double aleab = std::log(availableEnergy); 00348 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00349 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00350 00351 // normalization constant for kno-distribution. 00352 // calculate first the sum of all constants, check for numerical problems. 00353 G4double test, dum, anpn = 0.0; 00354 00355 for (nt=1; nt<=numSec; nt++) { 00356 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00357 dum = pi*nt/(2.0*n*n); 00358 if (std::fabs(dum) < 1.0) { 00359 if( test >= 1.0e-10 )anpn += dum*test; 00360 } else { 00361 anpn += dum*test; 00362 } 00363 } 00364 00365 G4double ran = G4UniformRand(); 00366 G4double excs = 0.0; 00367 if( targetCode == protonCode ) 00368 { 00369 counter = -1; 00370 for (npos=0; npos<numSec/3; npos++) { 00371 for (nneg=Imax(0,npos-2); nneg<=npos; nneg++) { 00372 for (nzero=0; nzero<numSec/3; nzero++) { 00373 if (++counter < numMul) { 00374 nt = npos+nneg+nzero; 00375 if ( (nt>0) && (nt<=numSec) ) { 00376 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00377 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 00378 if (std::fabs(dum) < 1.0) { 00379 if( test >= 1.0e-10 )excs += dum*test; 00380 } else { 00381 excs += dum*test; 00382 } 00383 if (ran < excs) goto outOfLoop; //------------------> 00384 } 00385 } 00386 } 00387 } 00388 } 00389 00390 // 3 previous loops continued to the end 00391 inElastic = false; // quasi-elastic scattering 00392 return; 00393 } 00394 else 00395 { // target must be a neutron 00396 counter = -1; 00397 for (npos=0; npos<numSec/3; npos++) { 00398 for (nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++) { 00399 for (nzero=0; nzero<numSec/3; nzero++) { 00400 if (++counter < numMul) { 00401 nt = npos+nneg+nzero; 00402 if ( (nt>=1) && (nt<=numSec) ) { 00403 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00404 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 00405 if (std::fabs(dum) < 1.0) { 00406 if( test >= 1.0e-10 )excs += dum*test; 00407 } else { 00408 excs += dum*test; 00409 } 00410 if (ran < excs) goto outOfLoop; // -------------> 00411 } 00412 } 00413 } 00414 } 00415 } 00416 // 3 previous loops continued to the end 00417 inElastic = false; // quasi-elastic scattering. 00418 return; 00419 } 00420 } 00421 outOfLoop: // <----------------------------------------------- 00422 00423 if( targetCode == protonCode) 00424 { 00425 if( npos == nneg) 00426 { 00427 } 00428 else if (npos == (1+nneg)) 00429 { 00430 if( G4UniformRand() < 0.5) 00431 { 00432 pv[1] = Neutron; 00433 } 00434 else 00435 { 00436 pv[0] = Neutron; 00437 } 00438 } 00439 else 00440 { 00441 pv[0] = Neutron; 00442 pv[1] = Neutron; 00443 } 00444 } 00445 else 00446 { 00447 if( npos == nneg) 00448 { 00449 if( G4UniformRand() < 0.25) 00450 { 00451 pv[0] = Neutron; 00452 pv[1] = Proton; 00453 } 00454 else 00455 { 00456 } 00457 } 00458 else if ( npos == (1+nneg)) 00459 { 00460 pv[0] = Neutron; 00461 } 00462 else 00463 { 00464 pv[1] = Proton; 00465 } 00466 } 00467 00468 00469 nt = npos + nneg + nzero; 00470 while ( nt > 0) 00471 { 00472 G4double ran = G4UniformRand(); 00473 if ( ran < (G4double)npos/nt) 00474 { 00475 if( npos > 0 ) 00476 { pv[vecLen++] = PionPlus; 00477 npos--; 00478 } 00479 } 00480 else if ( ran < (G4double)(npos+nneg)/nt) 00481 { 00482 if( nneg > 0 ) 00483 { 00484 pv[vecLen++] = PionMinus; 00485 nneg--; 00486 } 00487 } 00488 else 00489 { 00490 if( nzero > 0 ) 00491 { 00492 pv[vecLen++] = PionZero; 00493 nzero--; 00494 } 00495 } 00496 nt = npos + nneg + nzero; 00497 } 00498 if (verboseLevel > 1) 00499 { 00500 G4cout << "Particles produced: " ; 00501 G4cout << pv[0].getName() << " " ; 00502 G4cout << pv[1].getName() << " " ; 00503 for (i=2; i < vecLen; i++) 00504 { 00505 G4cout << pv[i].getName() << " " ; 00506 } 00507 G4cout << G4endl; 00508 } 00509 return; 00510 }
G4int G4HEProtonInelastic::GetNumberOfSecondaries | ( | ) | [inline] |
Definition at line 75 of file G4HEProtonInelastic.hh.
References vecLength.
00075 {return vecLength;}
void G4HEProtonInelastic::ModelDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 42 of file G4HEProtonInelastic.cc.
00043 { 00044 outFile << "G4HEProtonInelastic is one of the High Energy\n" 00045 << "Parameterized (HEP) models used to implement inelastic\n" 00046 << "proton scattering from nuclei. It is a re-engineered\n" 00047 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n" 00048 << "initial collision products into backward- and forward-going\n" 00049 << "clusters which are then decayed into final state hadrons.\n" 00050 << "The model does not conserve energy on an event-by-event\n" 00051 << "basis. It may be applied to protons with initial energies\n" 00052 << "above 20 GeV.\n"; 00053 }
Definition at line 70 of file G4HEProtonInelastic.hh.
Referenced by ApplyYourself(), G4HEProtonInelastic(), and GetNumberOfSecondaries().