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