#include <G4HEAntiProtonInelastic.hh>
Inheritance diagram for G4HEAntiProtonInelastic:
Public Member Functions | |
G4HEAntiProtonInelastic () | |
~G4HEAntiProtonInelastic () | |
virtual void | ModelDescription (std::ostream &) const |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
G4int | GetNumberOfSecondaries () |
void | FirstIntInCasAntiProton (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 G4HEAntiProtonInelastic.hh.
G4HEAntiProtonInelastic::G4HEAntiProtonInelastic | ( | ) | [inline] |
Definition at line 55 of file G4HEAntiProtonInelastic.hh.
References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, vecLength, and G4HEInelastic::verboseLevel.
00055 : G4HEInelastic("G4HEAntiProtonInelastic") 00056 { 00057 vecLength = 0; 00058 theMinEnergy = 20*CLHEP::GeV; 00059 theMaxEnergy = 10*CLHEP::TeV; 00060 MAXPART = 2048; 00061 verboseLevel = 0; 00062 G4cout << "WARNING: model G4HEAntiProtonInelastic is being deprecated and will\n" 00063 << "disappear in Geant4 version 10.0" << G4endl; 00064 }
G4HEAntiProtonInelastic::~G4HEAntiProtonInelastic | ( | ) | [inline] |
G4HadFinalState * G4HEAntiProtonInelastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | targetNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 58 of file G4HEAntiProtonInelastic.cc.
References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasAntiProton(), 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.
00060 { 00061 G4HEVector* pv = new G4HEVector[MAXPART]; 00062 const G4HadProjectile* aParticle = &aTrack; 00063 const G4double atomicWeight = targetNucleus.GetA_asInt(); 00064 const G4double atomicNumber = targetNucleus.GetZ_asInt(); 00065 G4HEVector incidentParticle(aParticle); 00066 00067 G4int incidentCode = incidentParticle.getCode(); 00068 G4double incidentMass = incidentParticle.getMass(); 00069 G4double incidentTotalEnergy = incidentParticle.getEnergy(); 00070 00071 // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00072 // DHW 19 may 2011: variable set but not used 00073 00074 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass; 00075 00076 if (incidentKineticEnergy < 1.) 00077 G4cout << "GHEAntiProtonInelastic: incident energy < 1 GeV" << G4endl; 00078 00079 if (verboseLevel > 1) { 00080 G4cout << "G4HEAntiProtonInelastic::ApplyYourself" << G4endl; 00081 G4cout << "incident particle " << incidentParticle.getName() 00082 << "mass " << incidentMass 00083 << "kinetic energy " << incidentKineticEnergy 00084 << G4endl; 00085 G4cout << "target material with (A,Z) = (" 00086 << atomicWeight << "," << atomicNumber << ")" << G4endl; 00087 } 00088 00089 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 00090 atomicWeight, atomicNumber); 00091 if (verboseLevel > 1) 00092 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 00093 00094 incidentKineticEnergy -= inelasticity; 00095 00096 G4double excitationEnergyGNP = 0.; 00097 G4double excitationEnergyDTA = 0.; 00098 00099 G4double excitation = NuclearExcitation(incidentKineticEnergy, 00100 atomicWeight, atomicNumber, 00101 excitationEnergyGNP, 00102 excitationEnergyDTA); 00103 if (verboseLevel > 1) 00104 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 00105 << excitationEnergyDTA << G4endl; 00106 00107 00108 incidentKineticEnergy -= excitation; 00109 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 00110 // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 00111 // *(incidentTotalEnergy+incidentMass)); 00112 // DHW 19 may 2011: variable set but not used 00113 00114 G4HEVector targetParticle; 00115 if (G4UniformRand() < atomicNumber/atomicWeight) { 00116 targetParticle.setDefinition("Proton"); 00117 } else { 00118 targetParticle.setDefinition("Neutron"); 00119 } 00120 00121 G4double targetMass = targetParticle.getMass(); 00122 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass 00123 + targetMass*targetMass 00124 + 2.0*targetMass*incidentTotalEnergy); 00125 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 00126 00127 G4bool inElastic = true; 00128 vecLength = 0; 00129 00130 if (verboseLevel > 1) 00131 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 00132 << incidentCode << G4endl; 00133 00134 G4bool successful = false; 00135 00136 FirstIntInCasAntiProton(inElastic, availableEnergy, pv, vecLength, 00137 incidentParticle, targetParticle, atomicWeight); 00138 00139 if (verboseLevel > 1) 00140 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 00141 00142 if ((vecLength > 0) && (availableEnergy > 1.)) 00143 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 00144 pv, vecLength, 00145 incidentParticle, targetParticle); 00146 HighEnergyCascading(successful, pv, vecLength, 00147 excitationEnergyGNP, excitationEnergyDTA, 00148 incidentParticle, targetParticle, 00149 atomicWeight, atomicNumber); 00150 if (!successful) 00151 HighEnergyClusterProduction(successful, pv, vecLength, 00152 excitationEnergyGNP, excitationEnergyDTA, 00153 incidentParticle, targetParticle, 00154 atomicWeight, atomicNumber); 00155 if (!successful) 00156 MediumEnergyCascading(successful, pv, vecLength, 00157 excitationEnergyGNP, excitationEnergyDTA, 00158 incidentParticle, targetParticle, 00159 atomicWeight, atomicNumber); 00160 00161 if (!successful) 00162 MediumEnergyClusterProduction(successful, pv, vecLength, 00163 excitationEnergyGNP, excitationEnergyDTA, 00164 incidentParticle, targetParticle, 00165 atomicWeight, atomicNumber); 00166 if (!successful) 00167 QuasiElasticScattering(successful, pv, vecLength, 00168 excitationEnergyGNP, excitationEnergyDTA, 00169 incidentParticle, targetParticle, 00170 atomicWeight, atomicNumber); 00171 if (!successful) 00172 ElasticScattering(successful, pv, vecLength, 00173 incidentParticle, 00174 atomicWeight, atomicNumber); 00175 00176 if (!successful) 00177 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 00178 << G4endl; 00179 00180 FillParticleChange(pv, vecLength); 00181 delete [] pv; 00182 theParticleChange.SetStatusChange(stopAndKill); 00183 return & theParticleChange; 00184 }
void G4HEAntiProtonInelastic::FirstIntInCasAntiProton | ( | G4bool & | inElastic, | |
const G4double | availableEnergy, | |||
G4HEVector | pv[], | |||
G4int & | vecLen, | |||
const G4HEVector & | incidentParticle, | |||
const G4HEVector & | targetParticle, | |||
const G4double | atomicWeight | |||
) |
Definition at line 188 of file G4HEAntiProtonInelastic.cc.
References G4HEInelastic::Amax(), G4HEInelastic::Amin(), G4HEInelastic::AntiNeutron, G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), G4HEInelastic::Imax(), G4HEInelastic::Imin(), 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().
00203 { 00204 static const G4double expxu = 82; // upper bound for arg. of exp 00205 static const G4double expxl = -expxu; // lower bound for arg. of exp 00206 00207 static const G4double protb = 0.7; 00208 static const G4double neutb = 0.7; 00209 static const G4double c = 1.25; 00210 00211 static const G4int numMul = 1200; 00212 static const G4int numMulAn = 400; 00213 static const G4int numSec = 60; 00214 00215 G4int neutronCode = Neutron.getCode(); 00216 G4int protonCode = Proton.getCode(); 00217 00218 G4int targetCode = targetParticle.getCode(); 00219 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00220 00221 static G4bool first = true; 00222 static G4double protmul[numMul], protnorm[numSec]; // proton constants 00223 static G4double protmulAn[numMulAn],protnormAn[numSec]; 00224 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 00225 static G4double neutmulAn[numMulAn],neutnormAn[numSec]; 00226 00227 // misc. local variables 00228 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0 00229 00230 G4int i, counter, nt, npos, nneg, nzero; 00231 00232 if( first ) 00233 { // compute normalization constants, this will only be done once 00234 first = false; 00235 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0; 00236 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0; 00237 counter = -1; 00238 for( npos=0; npos<(numSec/3); npos++ ) 00239 { 00240 for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ ) 00241 { 00242 for( nzero=0; nzero<numSec/3; nzero++ ) 00243 { 00244 if( ++counter < numMul ) 00245 { 00246 nt = npos+nneg+nzero; 00247 if( (nt>0) && (nt<=numSec) ) 00248 { 00249 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c); 00250 protnorm[nt-1] += protmul[counter]; 00251 } 00252 } 00253 } 00254 } 00255 } 00256 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 00257 00258 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 00259 counter = -1; 00260 for( npos=0; npos<numSec/3; npos++ ) 00261 { 00262 for( nneg=npos; nneg<=(npos+2); nneg++ ) 00263 { 00264 for( nzero=0; nzero<numSec/3; nzero++ ) 00265 { 00266 if( ++counter < numMul ) 00267 { 00268 nt = npos+nneg+nzero; 00269 if( (nt>0) && (nt<=numSec) ) 00270 { 00271 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c); 00272 neutnorm[nt-1] += neutmul[counter]; 00273 } 00274 } 00275 } 00276 } 00277 } 00278 for( i=0; i<numSec; i++ ) 00279 { 00280 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 00281 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 00282 } 00283 // annihilation 00284 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0; 00285 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0; 00286 counter = -1; 00287 for( npos=1; npos<(numSec/3); npos++ ) 00288 { 00289 nneg = npos; 00290 for( nzero=0; nzero<numSec/3; nzero++ ) 00291 { 00292 if( ++counter < numMulAn ) 00293 { 00294 nt = npos+nneg+nzero; 00295 if( (nt>0) && (nt<=numSec) ) 00296 { 00297 protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c); 00298 protnormAn[nt-1] += protmulAn[counter]; 00299 } 00300 } 00301 } 00302 } 00303 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0; 00304 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0; 00305 counter = -1; 00306 for( npos=1; npos<numSec/3; npos++ ) 00307 { 00308 nneg = npos+1; 00309 for( nzero=0; nzero<numSec/3; nzero++ ) 00310 { 00311 if( ++counter < numMulAn ) 00312 { 00313 nt = npos+nneg+nzero; 00314 if( (nt>0) && (nt<=numSec) ) 00315 { 00316 neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c); 00317 neutnormAn[nt-1] += neutmulAn[counter]; 00318 } 00319 } 00320 } 00321 } 00322 for( i=0; i<numSec; i++ ) 00323 { 00324 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i]; 00325 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i]; 00326 } 00327 } // end of initialization 00328 00329 00330 // initialize the first two places 00331 // the same as beam and target 00332 pv[0] = incidentParticle; 00333 pv[1] = targetParticle; 00334 vecLen = 2; 00335 00336 if( !inElastic ) 00337 { // pb p --> nb n 00338 if( targetCode == protonCode ) 00339 { 00340 G4double cech[] = {0.14, 0.170, 0.180, 0.180, 0.180, 0.170, 0.170, 0.160, 0.155, 0.145, 00341 0.11, 0.082, 0.065, 0.050, 0.041, 0.035, 0.028, 0.024, 0.010, 0.000}; 00342 00343 G4int iplab = G4int( incidentTotalMomentum*10.); 00344 if (iplab > 9) iplab = Imin(19, G4int( incidentTotalMomentum) + 9); 00345 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 00346 { // charge exchange pi+ n -> pi0 p 00347 pv[0] = AntiNeutron; 00348 pv[1] = Neutron; 00349 } 00350 } 00351 return; 00352 } 00353 else if (availableEnergy <= PionPlus.getMass()) 00354 return; 00355 00356 // inelastic scattering 00357 00358 npos = 0; nneg = 0; nzero = 0; 00359 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.90, 00360 0.60, 0.52, 0.47, 0.44, 0.41, 0.39, 0.37, 0.35, 0.34, 0.24, 00361 0.19, 0.15, 0.12, 0.10, 0.09, 0.07, 0.06, 0.05, 0.00}; 00362 G4int iplab = G4int( incidentTotalMomentum*10.); 00363 if ( iplab > 9) iplab = 9 + G4int( incidentTotalMomentum); 00364 if ( iplab > 18) iplab = 18 + G4int( incidentTotalMomentum*10.); 00365 iplab = Imin(28, iplab); 00366 00367 if ( G4UniformRand() > anhl[iplab] ) 00368 { 00369 00370 G4double eab = availableEnergy; 00371 G4int ieab = G4int( eab*5.0 ); 00372 00373 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98}; 00374 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) ) 00375 { 00376 // suppress high multiplicity events at low momentum 00377 // only one additional pion will be produced 00378 G4double w0, wp, wm, wt, ran; 00379 if( targetCode == neutronCode ) // target is a neutron 00380 { 00381 w0 = - sqr(1.+neutb)/(2.*c*c); 00382 w0 = std::exp(w0); 00383 wm = - sqr(-1.+neutb)/(2.*c*c); 00384 wm = std::exp(wm); 00385 if( G4UniformRand() < w0/(w0+wm) ) 00386 { npos = 0; nneg = 0; nzero = 1; } 00387 else 00388 { npos = 0; nneg = 1; nzero = 0; } 00389 } 00390 else 00391 { // target is a proton 00392 w0 = -sqr(1.+protb)/(2.*c*c); 00393 w0 = std::exp(w0); 00394 wp = w0; 00395 wm = -sqr(-1.+protb)/(2.*c*c); 00396 wm = std::exp(wm); 00397 wt = w0+wp+wm; 00398 wp = w0+wp; 00399 ran = G4UniformRand(); 00400 if( ran < w0/wt) 00401 { npos = 0; nneg = 0; nzero = 1; } 00402 else if( ran < wp/wt) 00403 { npos = 1; nneg = 0; nzero = 0; } 00404 else 00405 { npos = 0; nneg = 1; nzero = 0; } 00406 } 00407 } 00408 else 00409 { 00410 // number of total particles vs. centre of mass Energy - 2*proton mass 00411 00412 G4double aleab = std::log(availableEnergy); 00413 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00414 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00415 00416 // normalization constant for kno-distribution. 00417 // calculate first the sum of all constants, check for numerical problems. 00418 G4double test, dum, anpn = 0.0; 00419 00420 for (nt=1; nt<=numSec; nt++) { 00421 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00422 dum = pi*nt/(2.0*n*n); 00423 if (std::fabs(dum) < 1.0) { 00424 if( test >= 1.0e-10 )anpn += dum*test; 00425 } else { 00426 anpn += dum*test; 00427 } 00428 } 00429 00430 G4double ran = G4UniformRand(); 00431 G4double excs = 0.0; 00432 if( targetCode == protonCode ) 00433 { 00434 counter = -1; 00435 for( npos=0; npos<numSec/3; npos++ ) 00436 { 00437 for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ ) 00438 { 00439 for( nzero=0; nzero<numSec/3; nzero++ ) 00440 { 00441 if( ++counter < numMul ) 00442 { 00443 nt = npos+nneg+nzero; 00444 if ( (nt>0) && (nt<=numSec) ) { 00445 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00446 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 00447 if (std::fabs(dum) < 1.0) { 00448 if( test >= 1.0e-10 )excs += dum*test; 00449 } else { 00450 excs += dum*test; 00451 } 00452 00453 if (ran < excs) goto outOfLoop; //-----------------------> 00454 } 00455 } 00456 } 00457 } 00458 } 00459 00460 // 3 previous loops continued to the end 00461 inElastic = false; // quasi-elastic scattering 00462 return; 00463 } 00464 else 00465 { // target must be a neutron 00466 counter = -1; 00467 for( npos=0; npos<numSec/3; npos++ ) 00468 { 00469 for( nneg=npos; nneg<=(npos+2); nneg++ ) 00470 { 00471 for( nzero=0; nzero<numSec/3; nzero++ ) 00472 { 00473 if( ++counter < numMul ) 00474 { 00475 nt = npos+nneg+nzero; 00476 if ( (nt>=1) && (nt<=numSec) ) { 00477 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00478 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 00479 if (std::fabs(dum) < 1.0) { 00480 if( test >= 1.0e-10 )excs += dum*test; 00481 } else { 00482 excs += dum*test; 00483 } 00484 00485 if (ran < excs) goto outOfLoop; // --------------------------> 00486 } 00487 } 00488 } 00489 } 00490 } 00491 // 3 previous loops continued to the end 00492 inElastic = false; // quasi-elastic scattering. 00493 return; 00494 } 00495 } 00496 outOfLoop: // <------------------------------------------------------------------------ 00497 00498 if( targetCode == neutronCode) 00499 { 00500 if( npos == nneg) 00501 { 00502 } 00503 else if (npos == (nneg-1)) 00504 { 00505 if( G4UniformRand() < 0.5) 00506 { 00507 pv[1] = Proton; 00508 } 00509 else 00510 { 00511 pv[0] = AntiNeutron; 00512 } 00513 } 00514 else 00515 { 00516 pv[0] = AntiNeutron; 00517 pv[1] = Proton; 00518 } 00519 } 00520 else 00521 { 00522 if( npos == nneg) 00523 { 00524 if( G4UniformRand() < 0.25) 00525 { 00526 pv[0] = AntiNeutron; 00527 pv[1] = Neutron; 00528 } 00529 else 00530 { 00531 } 00532 } 00533 else if ( npos == (1+nneg)) 00534 { 00535 pv[1] = Neutron; 00536 } 00537 else 00538 { 00539 pv[0] = AntiNeutron; 00540 } 00541 } 00542 00543 } 00544 else // annihilation 00545 { 00546 if ( availableEnergy > 2. * PionPlus.getMass() ) 00547 { 00548 00549 G4double aleab = std::log(availableEnergy); 00550 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00551 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00552 00553 // normalization constant for kno-distribution. 00554 // calculate first the sum of all constants, check for numerical problems. 00555 G4double test, dum, anpn = 0.0; 00556 00557 for (nt=2; nt<=numSec; nt++) { 00558 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00559 dum = pi*nt/(2.0*n*n); 00560 if (std::fabs(dum) < 1.0) { 00561 if( test >= 1.0e-10 )anpn += dum*test; 00562 } else { 00563 anpn += dum*test; 00564 } 00565 } 00566 00567 G4double ran = G4UniformRand(); 00568 G4double excs = 0.0; 00569 if( targetCode == protonCode ) 00570 { 00571 counter = -1; 00572 for( npos=1; npos<numSec/3; npos++ ) 00573 { 00574 nneg = npos; 00575 for( nzero=0; nzero<numSec/3; nzero++ ) 00576 { 00577 if( ++counter < numMulAn ) 00578 { 00579 nt = npos+nneg+nzero; 00580 if ( (nt>0) && (nt<=numSec) ) { 00581 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00582 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n); 00583 if (std::fabs(dum) < 1.0) { 00584 if( test >= 1.0e-10 )excs += dum*test; 00585 } else { 00586 excs += dum*test; 00587 } 00588 00589 if (ran < excs) goto outOfLoopAn; //-----------------------> 00590 } 00591 } 00592 } 00593 } 00594 // 3 previous loops continued to the end 00595 inElastic = false; // quasi-elastic scattering 00596 return; 00597 } 00598 else 00599 { // target must be a neutron 00600 counter = -1; 00601 for( npos=1; npos<numSec/3; npos++ ) 00602 { 00603 nneg = npos+1; 00604 for( nzero=0; nzero<numSec/3; nzero++ ) 00605 { 00606 if( ++counter < numMulAn ) 00607 { 00608 nt = npos+nneg+nzero; 00609 if ( (nt>=1) && (nt<=numSec) ) { 00610 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00611 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n); 00612 if (std::fabs(dum) < 1.0) { 00613 if( test >= 1.0e-10 )excs += dum*test; 00614 } else { 00615 excs += dum*test; 00616 } 00617 00618 if (ran < excs) goto outOfLoopAn; // --------------------------> 00619 } 00620 } 00621 } 00622 } 00623 inElastic = false; // quasi-elastic scattering. 00624 return; 00625 } 00626 outOfLoopAn: // <------------------------------------------------------------------ 00627 vecLen = 0; 00628 } 00629 } 00630 00631 nt = npos + nneg + nzero; 00632 while ( nt > 0) 00633 { 00634 G4double ran = G4UniformRand(); 00635 if ( ran < (G4double)npos/nt) 00636 { 00637 if( npos > 0 ) 00638 { pv[vecLen++] = PionPlus; 00639 npos--; 00640 } 00641 } 00642 else if ( ran < (G4double)(npos+nneg)/nt) 00643 { 00644 if( nneg > 0 ) 00645 { 00646 pv[vecLen++] = PionMinus; 00647 nneg--; 00648 } 00649 } 00650 else 00651 { 00652 if( nzero > 0 ) 00653 { 00654 pv[vecLen++] = PionZero; 00655 nzero--; 00656 } 00657 } 00658 nt = npos + nneg + nzero; 00659 } 00660 if (verboseLevel > 1) 00661 { 00662 G4cout << "Particles produced: " ; 00663 G4cout << pv[0].getName() << " " ; 00664 G4cout << pv[1].getName() << " " ; 00665 for (i=2; i < vecLen; i++) 00666 { 00667 G4cout << pv[i].getName() << " " ; 00668 } 00669 G4cout << G4endl; 00670 } 00671 return; 00672 }
G4int G4HEAntiProtonInelastic::GetNumberOfSecondaries | ( | ) | [inline] |
Definition at line 75 of file G4HEAntiProtonInelastic.hh.
References vecLength.
00075 {return vecLength;}
void G4HEAntiProtonInelastic::ModelDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 43 of file G4HEAntiProtonInelastic.cc.
00044 { 00045 outFile << "G4HEAntiProtonInelastic is one of the High Energy\n" 00046 << "Parameterized (HEP) models used to implement inelastic\n" 00047 << "anti-proton scattering from nuclei. It is a re-engineered\n" 00048 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n" 00049 << "initial collision products into backward- and forward-going\n" 00050 << "clusters which are then decayed into final state hadrons.\n" 00051 << "The model does not conserve energy on an event-by-event\n" 00052 << "basis. It may be applied to anti-protons with initial\n" 00053 << "energies above 20 GeV.\n"; 00054 }
Definition at line 70 of file G4HEAntiProtonInelastic.hh.
Referenced by ApplyYourself(), G4HEAntiProtonInelastic(), and GetNumberOfSecondaries().