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