#include <G4HEAntiSigmaPlusInelastic.hh>
Inheritance diagram for G4HEAntiSigmaPlusInelastic:
Public Member Functions | |
G4HEAntiSigmaPlusInelastic () | |
~G4HEAntiSigmaPlusInelastic () | |
virtual void | ModelDescription (std::ostream &) const |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
G4int | GetNumberOfSecondaries () |
void | FirstIntInCasAntiSigmaPlus (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 53 of file G4HEAntiSigmaPlusInelastic.hh.
G4HEAntiSigmaPlusInelastic::G4HEAntiSigmaPlusInelastic | ( | ) | [inline] |
Definition at line 56 of file G4HEAntiSigmaPlusInelastic.hh.
References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, vecLength, and G4HEInelastic::verboseLevel.
00056 : G4HEInelastic("G4HEAntiSigmaPlusInelastic") 00057 { 00058 vecLength = 0; 00059 theMinEnergy = 20*CLHEP::GeV; 00060 theMaxEnergy = 10*CLHEP::TeV; 00061 MAXPART = 2048; 00062 verboseLevel = 0; 00063 G4cout << "WARNING: model G4HEAntiSigmaPlusInelastic is being deprecated and will\n" 00064 << "disappear in Geant4 version 10.0" << G4endl; 00065 }
G4HEAntiSigmaPlusInelastic::~G4HEAntiSigmaPlusInelastic | ( | ) | [inline] |
G4HadFinalState * G4HEAntiSigmaPlusInelastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | targetNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 58 of file G4HEAntiSigmaPlusInelastic.cc.
References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasAntiSigmaPlus(), 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 << "GHEAntiSigmaPlusInelastic: incident energy < 1 GeV" << G4endl; 00078 00079 if (verboseLevel > 1) { 00080 G4cout << "G4HEAntiSigmaPlusInelastic::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 incidentKineticEnergy -= excitation; 00108 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 00109 // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 00110 // *(incidentTotalEnergy+incidentMass)); 00111 // DHW 19 May 2011: variable set but not used 00112 00113 G4HEVector targetParticle; 00114 if (G4UniformRand() < atomicNumber/atomicWeight) { 00115 targetParticle.setDefinition("Proton"); 00116 } else { 00117 targetParticle.setDefinition("Neutron"); 00118 } 00119 00120 G4double targetMass = targetParticle.getMass(); 00121 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass 00122 + targetMass*targetMass 00123 + 2.0*targetMass*incidentTotalEnergy); 00124 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 00125 00126 G4bool inElastic = true; 00127 vecLength = 0; 00128 00129 if (verboseLevel > 1) 00130 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 00131 << incidentCode << G4endl; 00132 00133 G4bool successful = false; 00134 00135 FirstIntInCasAntiSigmaPlus(inElastic, availableEnergy, pv, vecLength, 00136 incidentParticle, targetParticle, atomicWeight); 00137 00138 if (verboseLevel > 1) 00139 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 00140 00141 if ((vecLength > 0) && (availableEnergy > 1.)) 00142 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 00143 pv, vecLength, 00144 incidentParticle, targetParticle); 00145 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 G4HEAntiSigmaPlusInelastic::FirstIntInCasAntiSigmaPlus | ( | G4bool & | inElastic, | |
const G4double | availableEnergy, | |||
G4HEVector | pv[], | |||
G4int & | vecLen, | |||
const G4HEVector & | incidentParticle, | |||
const G4HEVector & | targetParticle, | |||
const G4double | atomicWeight | |||
) |
Definition at line 188 of file G4HEAntiSigmaPlusInelastic.cc.
References G4HEInelastic::AntiLambda, G4HEInelastic::AntiSigmaPlus, G4HEInelastic::AntiSigmaZero, G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), CLHEP::detail::n, G4HEInelastic::Neutron, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, 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 protonCode = Proton.getCode(); 00216 00217 G4int targetCode = targetParticle.getCode(); 00218 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00219 00220 static G4bool first = true; 00221 static G4double protmul[numMul], protnorm[numSec]; // proton constants 00222 static G4double protmulAn[numMulAn],protnormAn[numSec]; 00223 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 00224 static G4double neutmulAn[numMulAn],neutnormAn[numSec]; 00225 00226 // misc. local variables 00227 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0 00228 00229 G4int i, counter, nt, npos, nneg, nzero; 00230 00231 if (first) { 00232 // compute normalization constants, this will only be done once 00233 first = false; 00234 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0; 00235 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0; 00236 counter = -1; 00237 for (npos = 0; npos < (numSec/3); npos++) { 00238 for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) { 00239 for (nzero = 0; nzero < numSec/3; nzero++) { 00240 if (++counter < numMul) { 00241 nt = npos+nneg+nzero; 00242 if ((nt>0) && (nt<=numSec) ) { 00243 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c); 00244 protnorm[nt-1] += protmul[counter]; 00245 } 00246 } 00247 } 00248 } 00249 } 00250 00251 for (i = 0; i < numMul; i++) neutmul[i] = 0.0; 00252 for (i = 0; i < numSec; i++) neutnorm[i] = 0.0; 00253 counter = -1; 00254 for (npos = 0; npos < numSec/3; npos++) { 00255 for (nneg = npos; nneg <= (npos+2); nneg++) { 00256 for (nzero = 0; nzero < numSec/3; nzero++) { 00257 if (++counter < numMul) { 00258 nt = npos+nneg+nzero; 00259 if ((nt>0) && (nt<=numSec) ) { 00260 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c); 00261 neutnorm[nt-1] += neutmul[counter]; 00262 } 00263 } 00264 } 00265 } 00266 } 00267 00268 for (i = 0; i < numSec; i++) { 00269 if (protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 00270 if (neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 00271 } 00272 00273 // annihilation 00274 for (i = 0; i < numMulAn; i++) protmulAn[i] = 0.0; 00275 for (i = 0; i < numSec; i++) protnormAn[i] = 0.0; 00276 counter = -1; 00277 for (npos = 1; npos < (numSec/3); npos++) { 00278 nneg = npos; 00279 for (nzero = 0; nzero < numSec/3; nzero++) { 00280 if (++counter < numMulAn) { 00281 nt = npos+nneg+nzero; 00282 if ((nt>1) && (nt<=numSec) ) { 00283 protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c); 00284 protnormAn[nt-1] += protmulAn[counter]; 00285 } 00286 } 00287 } 00288 } 00289 00290 for (i = 0; i < numMulAn; i++) neutmulAn[i] = 0.0; 00291 for (i = 0; i < numSec; i++) neutnormAn[i] = 0.0; 00292 counter = -1; 00293 for (npos = 0; npos < numSec/3; npos++) { 00294 nneg = npos+1; 00295 for (nzero = 0; nzero < numSec/3; nzero++) { 00296 if (++counter < numMulAn) { 00297 nt = npos+nneg+nzero; 00298 if ((nt>1) && (nt<=numSec) ) { 00299 neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c); 00300 neutnormAn[nt-1] += neutmulAn[counter]; 00301 } 00302 } 00303 } 00304 } 00305 for (i = 0; i < numSec; i++) { 00306 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i]; 00307 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i]; 00308 } 00309 } // end of initialization 00310 00311 00312 // initialize the first two places the same as beam and target 00313 pv[0] = incidentParticle; 00314 pv[1] = targetParticle; 00315 vecLen = 2; 00316 00317 if (!inElastic) { // some two-body reactions 00318 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.}; 00319 00320 G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5)); 00321 if (G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) { 00322 G4double ran = G4UniformRand(); 00323 00324 if (targetCode == protonCode) { 00325 if (ran < 0.2) { 00326 pv[0] = Proton; 00327 pv[1] = AntiSigmaPlus; 00328 } else if (ran < 0.4) { 00329 pv[0] = AntiLambda; 00330 pv[1] = Neutron; 00331 } else if (ran < 0.6) { 00332 pv[0] = Neutron; 00333 pv[1] = AntiLambda; 00334 } else if (ran < 0.8) { 00335 pv[0] = Neutron; 00336 pv[1] = AntiSigmaZero; 00337 } else { 00338 pv[0] = AntiSigmaZero; 00339 pv[1] = Neutron; 00340 } 00341 } else { 00342 pv[0] = Neutron; 00343 pv[1] = AntiSigmaPlus; 00344 } 00345 } 00346 00347 return; 00348 } 00349 else if (availableEnergy <= PionPlus.getMass()) return; 00350 00351 // inelastic scattering 00352 npos = 0; nneg = 0; nzero = 0; 00353 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88, 00354 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40, 00355 0.39, 0.36, 0.33, 0.10, 0.01}; 00356 G4int iplab = G4int( incidentTotalMomentum*10.); 00357 if ( iplab > 9) iplab = 10 + G4int( (incidentTotalMomentum -1.)*5. ); 00358 if ( iplab > 14) iplab = 15 + G4int( incidentTotalMomentum -2. ); 00359 if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.); 00360 iplab = std::min(24, iplab); 00361 00362 if ( G4UniformRand() > anhl[iplab] ) { // non- annihilation channels 00363 00364 // number of total particles vs. centre of mass Energy - 2*proton mass 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( std::min( expxu, std::max( 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 { 00390 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ ) 00391 { 00392 for( nzero=0; nzero<numSec/3; nzero++ ) 00393 { 00394 if( ++counter < numMul ) 00395 { 00396 nt = npos+nneg+nzero; 00397 if ( (nt>0) && (nt<=numSec) ) { 00398 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00399 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 00400 if (std::fabs(dum) < 1.0) { 00401 if( test >= 1.0e-10 )excs += dum*test; 00402 } else { 00403 excs += dum*test; 00404 } 00405 00406 if (ran < excs) goto outOfLoop; //-----------------------> 00407 } 00408 } 00409 } 00410 } 00411 } 00412 00413 // 3 previous loops continued to the end 00414 inElastic = false; // quasi-elastic scattering 00415 return; 00416 } 00417 else 00418 { // target must be a neutron 00419 counter = -1; 00420 for( npos=0; npos<numSec/3; npos++ ) 00421 { 00422 for( nneg=npos; nneg<=(npos+2); nneg++ ) 00423 { 00424 for( nzero=0; nzero<numSec/3; nzero++ ) 00425 { 00426 if( ++counter < numMul ) 00427 { 00428 nt = npos+nneg+nzero; 00429 if ( (nt>0) && (nt<=numSec) ) { 00430 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00431 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 00432 if (std::fabs(dum) < 1.0) { 00433 if( test >= 1.0e-10 )excs += dum*test; 00434 } else { 00435 excs += dum*test; 00436 } 00437 00438 if (ran < excs) goto outOfLoop; // --------------------------> 00439 } 00440 } 00441 } 00442 } 00443 } 00444 // 3 previous loops continued to the end 00445 inElastic = false; // quasi-elastic scattering. 00446 return; 00447 } 00448 00449 outOfLoop: // <------------------------------------------------------------------------ 00450 00451 ran = G4UniformRand(); 00452 00453 if( targetCode == protonCode) 00454 { 00455 if( npos == nneg) 00456 { 00457 if (ran < 0.50) 00458 { 00459 } 00460 else if (ran < 0.75) 00461 { 00462 pv[0] = AntiSigmaZero; 00463 pv[1] = Neutron; 00464 } 00465 else 00466 { 00467 pv[0] = AntiLambda; 00468 pv[1] = Neutron; 00469 } 00470 } 00471 else if (npos == (nneg-1)) 00472 { 00473 if( ran < 0.50) 00474 { 00475 pv[0] = AntiLambda; 00476 } 00477 else 00478 { 00479 pv[0] = AntiSigmaZero; 00480 } 00481 } 00482 else 00483 { 00484 pv[1] = Neutron; 00485 } 00486 } 00487 else 00488 { 00489 if( npos == nneg) 00490 { 00491 } 00492 else if ( npos == (nneg-1)) 00493 { 00494 if (ran < 0.5) 00495 { 00496 pv[1] = Proton; 00497 } 00498 else if (ran < 0.75) 00499 { 00500 pv[0] = AntiLambda; 00501 } 00502 else 00503 { 00504 pv[0] = AntiSigmaZero; 00505 } 00506 } 00507 else 00508 { 00509 if (ran < 0.5) 00510 { 00511 pv[0] = AntiLambda; 00512 pv[1] = Proton; 00513 } 00514 else 00515 { 00516 pv[0] = AntiSigmaZero; 00517 pv[1] = Proton; 00518 } 00519 } 00520 } 00521 } 00522 else // annihilation 00523 { 00524 if ( availableEnergy > 2. * PionPlus.getMass() ) 00525 { 00526 00527 G4double aleab = std::log(availableEnergy); 00528 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00529 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00530 00531 // normalization constant for kno-distribution. 00532 // calculate first the sum of all constants, check for numerical problems. 00533 G4double test, dum, anpn = 0.0; 00534 00535 for (nt=2; nt<=numSec; nt++) { 00536 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00537 dum = pi*nt/(2.0*n*n); 00538 if (std::fabs(dum) < 1.0) { 00539 if( test >= 1.0e-10 )anpn += dum*test; 00540 } else { 00541 anpn += dum*test; 00542 } 00543 } 00544 00545 G4double ran = G4UniformRand(); 00546 G4double excs = 0.0; 00547 if( targetCode == protonCode ) 00548 { 00549 counter = -1; 00550 for( npos=1; npos<numSec/3; npos++ ) 00551 { 00552 nneg = npos; 00553 for( nzero=0; nzero<numSec/3; nzero++ ) 00554 { 00555 if( ++counter < numMulAn ) 00556 { 00557 nt = npos+nneg+nzero; 00558 if ( (nt>1) && (nt<=numSec) ) { 00559 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00560 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n); 00561 if (std::fabs(dum) < 1.0) { 00562 if( test >= 1.0e-10 )excs += dum*test; 00563 } else { 00564 excs += dum*test; 00565 } 00566 00567 if (ran < excs) goto outOfLoopAn; //-----------------------> 00568 } 00569 } 00570 } 00571 } 00572 // 3 previous loops continued to the end 00573 inElastic = false; // quasi-elastic scattering 00574 return; 00575 } 00576 else 00577 { // target must be a neutron 00578 counter = -1; 00579 for( npos=0; npos<numSec/3; npos++ ) 00580 { 00581 nneg = npos+1; 00582 for( nzero=0; nzero<numSec/3; nzero++ ) 00583 { 00584 if( ++counter < numMulAn ) 00585 { 00586 nt = npos+nneg+nzero; 00587 if ( (nt>1) && (nt<=numSec) ) { 00588 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00589 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n); 00590 if (std::fabs(dum) < 1.0) { 00591 if( test >= 1.0e-10 )excs += dum*test; 00592 } else { 00593 excs += dum*test; 00594 } 00595 00596 if (ran < excs) goto outOfLoopAn; // --------------------------> 00597 } 00598 } 00599 } 00600 } 00601 inElastic = false; // quasi-elastic scattering. 00602 return; 00603 } 00604 outOfLoopAn: // <---------------------------------------- 00605 vecLen = 0; 00606 } 00607 } 00608 00609 nt = npos + nneg + nzero; 00610 while ( nt > 0) 00611 { 00612 G4double ran = G4UniformRand(); 00613 if ( ran < (G4double)npos/nt) 00614 { 00615 if( npos > 0 ) 00616 { pv[vecLen++] = PionPlus; 00617 npos--; 00618 } 00619 } 00620 else if ( ran < (G4double)(npos+nneg)/nt) 00621 { 00622 if( nneg > 0 ) 00623 { 00624 pv[vecLen++] = PionMinus; 00625 nneg--; 00626 } 00627 } 00628 else 00629 { 00630 if( nzero > 0 ) 00631 { 00632 pv[vecLen++] = PionZero; 00633 nzero--; 00634 } 00635 } 00636 nt = npos + nneg + nzero; 00637 } 00638 if (verboseLevel > 1) 00639 { 00640 G4cout << "Particles produced: " ; 00641 G4cout << pv[0].getName() << " " ; 00642 G4cout << pv[1].getName() << " " ; 00643 for (i=2; i < vecLen; i++) 00644 { 00645 G4cout << pv[i].getName() << " " ; 00646 } 00647 G4cout << G4endl; 00648 } 00649 return; 00650 }
G4int G4HEAntiSigmaPlusInelastic::GetNumberOfSecondaries | ( | ) | [inline] |
Definition at line 76 of file G4HEAntiSigmaPlusInelastic.hh.
References vecLength.
00076 {return vecLength;}
void G4HEAntiSigmaPlusInelastic::ModelDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 43 of file G4HEAntiSigmaPlusInelastic.cc.
00044 { 00045 outFile << "G4HEAntiSigmaPlusInelastic is one of the High Energy\n" 00046 << "Parameterized (HEP) models used to implement inelastic\n" 00047 << "anti-Sigma+ 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-Sigma+ with initial\n" 00053 << "energies above 20 GeV.\n"; 00054 }
Definition at line 71 of file G4HEAntiSigmaPlusInelastic.hh.
Referenced by ApplyYourself(), G4HEAntiSigmaPlusInelastic(), and GetNumberOfSecondaries().