#include <G4HEAntiSigmaMinusInelastic.hh>
Inheritance diagram for G4HEAntiSigmaMinusInelastic:
Public Member Functions | |
| G4HEAntiSigmaMinusInelastic () | |
| ~G4HEAntiSigmaMinusInelastic () | |
| virtual void | ModelDescription (std::ostream &) const |
| G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
| G4int | GetNumberOfSecondaries () |
| void | FirstIntInCasAntiSigmaMinus (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 G4HEAntiSigmaMinusInelastic.hh.
| G4HEAntiSigmaMinusInelastic::G4HEAntiSigmaMinusInelastic | ( | ) | [inline] |
Definition at line 56 of file G4HEAntiSigmaMinusInelastic.hh.
References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, vecLength, and G4HEInelastic::verboseLevel.
00056 : G4HEInelastic("G4HEAntiSigmaMinusInelastic") 00057 { 00058 vecLength = 0; 00059 theMinEnergy = 20*CLHEP::GeV; 00060 theMaxEnergy = 10*CLHEP::TeV; 00061 MAXPART = 2048; 00062 verboseLevel = 0; 00063 G4cout << "WARNING: model G4HEAntiSigmaMinusInelastic is being deprecated and will\n" 00064 << "disappear in Geant4 version 10.0" << G4endl; 00065 }
| G4HEAntiSigmaMinusInelastic::~G4HEAntiSigmaMinusInelastic | ( | ) | [inline] |
| G4HadFinalState * G4HEAntiSigmaMinusInelastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
| G4Nucleus & | targetNucleus | |||
| ) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 58 of file G4HEAntiSigmaMinusInelastic.cc.
References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasAntiSigmaMinus(), 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 << "GHEAntiSigmaMinusInelastic: incident energy < 1 GeV" << G4endl; 00078 00079 if (verboseLevel > 1) { 00080 G4cout << "G4HEAntiSigmaMinusInelastic::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 FirstIntInCasAntiSigmaMinus(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 HighEnergyCascading(successful, pv, vecLength, 00146 excitationEnergyGNP, excitationEnergyDTA, 00147 incidentParticle, targetParticle, 00148 atomicWeight, atomicNumber); 00149 if (!successful) 00150 HighEnergyClusterProduction(successful, pv, vecLength, 00151 excitationEnergyGNP, excitationEnergyDTA, 00152 incidentParticle, targetParticle, 00153 atomicWeight, atomicNumber); 00154 if (!successful) 00155 MediumEnergyCascading(successful, pv, vecLength, 00156 excitationEnergyGNP, excitationEnergyDTA, 00157 incidentParticle, targetParticle, 00158 atomicWeight, atomicNumber); 00159 00160 if (!successful) 00161 MediumEnergyClusterProduction(successful, pv, vecLength, 00162 excitationEnergyGNP, excitationEnergyDTA, 00163 incidentParticle, targetParticle, 00164 atomicWeight, atomicNumber); 00165 if (!successful) 00166 QuasiElasticScattering(successful, pv, vecLength, 00167 excitationEnergyGNP, excitationEnergyDTA, 00168 incidentParticle, targetParticle, 00169 atomicWeight, atomicNumber); 00170 if (!successful) 00171 ElasticScattering(successful, pv, vecLength, 00172 incidentParticle, 00173 atomicWeight, atomicNumber); 00174 00175 if (!successful) 00176 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 00177 << G4endl; 00178 FillParticleChange(pv, vecLength); 00179 delete [] pv; 00180 theParticleChange.SetStatusChange(stopAndKill); 00181 return &theParticleChange; 00182 }
| void G4HEAntiSigmaMinusInelastic::FirstIntInCasAntiSigmaMinus | ( | G4bool & | inElastic, | |
| const G4double | availableEnergy, | |||
| G4HEVector | pv[], | |||
| G4int & | vecLen, | |||
| const G4HEVector & | incidentParticle, | |||
| const G4HEVector & | targetParticle, | |||
| const G4double | atomicWeight | |||
| ) |
Definition at line 186 of file G4HEAntiSigmaMinusInelastic.cc.
References G4HEInelastic::AntiLambda, G4HEInelastic::AntiSigmaMinus, G4HEInelastic::AntiSigmaZero, 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, and G4HEInelastic::verboseLevel.
Referenced by ApplyYourself().
00201 { 00202 static const G4double expxu = 82; // upper bound for arg. of exp 00203 static const G4double expxl = -expxu; // lower bound for arg. of exp 00204 00205 static const G4double protb = 0.7; 00206 static const G4double neutb = 0.7; 00207 static const G4double c = 1.25; 00208 00209 static const G4int numMul = 1200; 00210 static const G4int numMulAn = 400; 00211 static const G4int numSec = 60; 00212 00213 G4int neutronCode = Neutron.getCode(); 00214 G4int protonCode = Proton.getCode(); 00215 00216 G4int targetCode = targetParticle.getCode(); 00217 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00218 00219 static G4bool first = true; 00220 static G4double protmul[numMul], protnorm[numSec]; // proton constants 00221 static G4double protmulAn[numMulAn],protnormAn[numSec]; 00222 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 00223 static G4double neutmulAn[numMulAn],neutnormAn[numSec]; 00224 00225 // misc. local variables 00226 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0 00227 00228 G4int i, counter, nt, npos, nneg, nzero; 00229 00230 if( first ) 00231 { // compute normalization constants, this will only be done once 00232 first = false; 00233 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0; 00234 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0; 00235 counter = -1; 00236 for( npos=0; npos<(numSec/3); npos++ ) 00237 { 00238 for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ ) 00239 { 00240 for( nzero=0; nzero<numSec/3; nzero++ ) 00241 { 00242 if( ++counter < numMul ) 00243 { 00244 nt = npos+nneg+nzero; 00245 if( (nt>0) && (nt<=numSec) ) 00246 { 00247 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c); 00248 protnorm[nt-1] += protmul[counter]; 00249 } 00250 } 00251 } 00252 } 00253 } 00254 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 00255 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 00256 counter = -1; 00257 for( npos=0; npos<numSec/3; npos++ ) 00258 { 00259 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ ) 00260 { 00261 for( nzero=0; nzero<numSec/3; nzero++ ) 00262 { 00263 if( ++counter < numMul ) 00264 { 00265 nt = npos+nneg+nzero; 00266 if( (nt>0) && (nt<=numSec) ) 00267 { 00268 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c); 00269 neutnorm[nt-1] += neutmul[counter]; 00270 } 00271 } 00272 } 00273 } 00274 } 00275 for( i=0; i<numSec; i++ ) 00276 { 00277 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 00278 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 00279 } 00280 // annihilation 00281 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0; 00282 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0; 00283 counter = -1; 00284 for( npos=1; npos<(numSec/3); npos++ ) 00285 { 00286 nneg = std::max(0,npos-2); 00287 for( nzero=0; nzero<numSec/3; nzero++ ) 00288 { 00289 if( ++counter < numMulAn ) 00290 { 00291 nt = npos+nneg+nzero; 00292 if( (nt>1) && (nt<=numSec) ) 00293 { 00294 protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c); 00295 protnormAn[nt-1] += protmulAn[counter]; 00296 } 00297 } 00298 } 00299 } 00300 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0; 00301 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0; 00302 counter = -1; 00303 for( npos=0; npos<numSec/3; npos++ ) 00304 { 00305 nneg = npos-1; 00306 for( nzero=0; nzero<numSec/3; nzero++ ) 00307 { 00308 if( ++counter < numMulAn ) 00309 { 00310 nt = npos+nneg+nzero; 00311 if( (nt>1) && (nt<=numSec) ) 00312 { 00313 neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c); 00314 neutnormAn[nt-1] += neutmulAn[counter]; 00315 } 00316 } 00317 } 00318 } 00319 for( i=0; i<numSec; i++ ) 00320 { 00321 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i]; 00322 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i]; 00323 } 00324 } // end of initialization 00325 00326 00327 // initialize the first two places 00328 // the same as beam and target 00329 pv[0] = incidentParticle; 00330 pv[1] = targetParticle; 00331 vecLen = 2; 00332 00333 if( !inElastic ) 00334 { // some two-body reactions 00335 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.}; 00336 00337 G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5 )); 00338 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 00339 { 00340 G4double ran = G4UniformRand(); 00341 00342 if ( targetCode == neutronCode) 00343 { 00344 if(ran < 0.2) 00345 { 00346 pv[0] = AntiSigmaZero; 00347 pv[1] = Proton; 00348 } 00349 else if (ran < 0.4) 00350 { 00351 pv[0] = AntiLambda; 00352 pv[1] = Proton; 00353 } 00354 else if (ran < 0.6) 00355 { 00356 pv[0] = Proton; 00357 pv[1] = AntiLambda; 00358 } 00359 else if (ran < 0.8) 00360 { 00361 pv[0] = Proton; 00362 pv[1] = AntiSigmaZero; 00363 } 00364 else 00365 { 00366 pv[0] = Neutron; 00367 pv[1] = AntiSigmaMinus; 00368 } 00369 } 00370 else 00371 { 00372 pv[0] = Proton; 00373 pv[1] = AntiSigmaMinus; 00374 } 00375 } 00376 return; 00377 } 00378 else if (availableEnergy <= PionPlus.getMass()) 00379 return; 00380 00381 // inelastic scattering 00382 00383 npos = 0; nneg = 0; nzero = 0; 00384 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88, 00385 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40, 00386 0.39, 0.36, 0.33, 0.10, 0.01}; 00387 G4int iplab = G4int( incidentTotalMomentum*10.); 00388 if ( iplab > 9) iplab = 10 + G4int( (incidentTotalMomentum -1.)*5. ); 00389 if ( iplab > 14) iplab = 15 + G4int( incidentTotalMomentum -2. ); 00390 if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.); 00391 iplab = std::min(24, iplab); 00392 00393 if (G4UniformRand() > anhl[iplab]) { // non- annihilation channels 00394 00395 // number of total particles vs. centre of mass Energy - 2*proton mass 00396 G4double aleab = std::log(availableEnergy); 00397 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00398 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00399 00400 // normalization constant for kno-distribution. 00401 // calculate first the sum of all constants, check for numerical problems. 00402 G4double test, dum, anpn = 0.0; 00403 00404 for (nt=1; nt<=numSec; nt++) { 00405 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00406 dum = pi*nt/(2.0*n*n); 00407 if (std::fabs(dum) < 1.0) { 00408 if( test >= 1.0e-10 )anpn += dum*test; 00409 } else { 00410 anpn += dum*test; 00411 } 00412 } 00413 00414 G4double ran = G4UniformRand(); 00415 G4double excs = 0.0; 00416 if (targetCode == protonCode) { 00417 counter = -1; 00418 for (npos = 0; npos < numSec/3; npos++) { 00419 for (nneg = std::max(0,npos-2); nneg <= npos; nneg++) { 00420 for (nzero = 0; nzero < numSec/3; nzero++) { 00421 if (++counter < numMul) { 00422 nt = npos+nneg+nzero; 00423 if ((nt > 0) && (nt <= numSec) ) { 00424 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00425 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 00426 if (std::fabs(dum) < 1.0) { 00427 if( test >= 1.0e-10 )excs += dum*test; 00428 } else { 00429 excs += dum*test; 00430 } 00431 00432 if (ran < excs) goto outOfLoop; //-----------------------> 00433 } 00434 } 00435 } 00436 } 00437 } 00438 // 3 previous loops continued to the end 00439 inElastic = false; // quasi-elastic scattering 00440 return; 00441 } else { // target must be a neutron 00442 counter = -1; 00443 for( npos=0; npos<numSec/3; npos++ ) 00444 { 00445 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ ) 00446 { 00447 for( nzero=0; nzero<numSec/3; nzero++ ) 00448 { 00449 if( ++counter < numMul ) 00450 { 00451 nt = npos+nneg+nzero; 00452 if ( (nt>0) && (nt<=numSec) ) { 00453 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00454 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 00455 if (std::fabs(dum) < 1.0) { 00456 if( test >= 1.0e-10 )excs += dum*test; 00457 } else { 00458 excs += dum*test; 00459 } 00460 00461 if (ran < excs) goto outOfLoop; // --------------------------> 00462 } 00463 } 00464 } 00465 } 00466 } 00467 // 3 previous loops continued to the end 00468 inElastic = false; // quasi-elastic scattering. 00469 return; 00470 } 00471 00472 outOfLoop: // <------------------------------------------------------------------------ 00473 00474 ran = G4UniformRand(); 00475 00476 if( targetCode == protonCode) 00477 { 00478 if( npos == nneg) 00479 { 00480 } 00481 else if (npos == (nneg+1)) 00482 { 00483 if( ran < 0.50) 00484 { 00485 pv[1] = Neutron; 00486 } 00487 else if (ran < 0.75) 00488 { 00489 pv[0] = AntiSigmaZero; 00490 } 00491 else 00492 { 00493 pv[0] = AntiLambda; 00494 } 00495 } 00496 else 00497 { 00498 if (ran < 0.5) 00499 { 00500 pv[0] = AntiSigmaZero; 00501 pv[1] = Neutron; 00502 } 00503 else 00504 { 00505 pv[0] = AntiLambda; 00506 pv[1] = Neutron; 00507 } 00508 } 00509 } 00510 else 00511 { 00512 if( npos == nneg) 00513 { 00514 if (ran < 0.5) 00515 { 00516 } 00517 else if(ran < 0.75) 00518 { 00519 pv[0] = AntiSigmaZero; 00520 pv[1] = Proton; 00521 } 00522 else 00523 { 00524 pv[0] = AntiLambda; 00525 pv[1] = Proton; 00526 } 00527 } 00528 else if ( npos == (nneg+1)) 00529 { 00530 if (ran < 0.5) 00531 { 00532 pv[0] = AntiSigmaZero; 00533 } 00534 else 00535 { 00536 pv[0] = AntiLambda; 00537 } 00538 } 00539 else 00540 { 00541 pv[1] = Proton; 00542 } 00543 } 00544 00545 } 00546 else // annihilation 00547 { 00548 if ( availableEnergy > 2. * PionPlus.getMass() ) 00549 { 00550 00551 G4double aleab = std::log(availableEnergy); 00552 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00553 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00554 00555 // normalization constant for kno-distribution. 00556 // calculate first the sum of all constants, check for numerical problems. 00557 G4double test, dum, anpn = 0.0; 00558 00559 for (nt=2; nt<=numSec; nt++) { 00560 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00561 dum = pi*nt/(2.0*n*n); 00562 if (std::fabs(dum) < 1.0) { 00563 if( test >= 1.0e-10 )anpn += dum*test; 00564 } else { 00565 anpn += dum*test; 00566 } 00567 } 00568 00569 G4double ran = G4UniformRand(); 00570 G4double excs = 0.0; 00571 if( targetCode == protonCode ) 00572 { 00573 counter = -1; 00574 for( npos=2; npos<numSec/3; npos++ ) 00575 { 00576 nneg = npos-2; 00577 for( nzero=0; nzero<numSec/3; nzero++ ) 00578 { 00579 if( ++counter < numMulAn ) 00580 { 00581 nt = npos+nneg+nzero; 00582 if ( (nt>1) && (nt<=numSec) ) { 00583 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00584 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n); 00585 if (std::fabs(dum) < 1.0) { 00586 if( test >= 1.0e-10 )excs += dum*test; 00587 } else { 00588 excs += dum*test; 00589 } 00590 00591 if (ran < excs) goto outOfLoopAn; //-----------------------> 00592 } 00593 } 00594 } 00595 } 00596 // 3 previous loops continued to the end 00597 inElastic = false; // quasi-elastic scattering 00598 return; 00599 } 00600 else 00601 { // target must be a neutron 00602 counter = -1; 00603 for( npos=1; npos<numSec/3; npos++ ) 00604 { 00605 nneg = npos-1; 00606 for( nzero=0; nzero<numSec/3; nzero++ ) 00607 { 00608 if (++counter < numMulAn) { 00609 nt = npos+nneg+nzero; 00610 if ( (nt>1) && (nt<=numSec) ) { 00611 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00612 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n); 00613 if (std::fabs(dum) < 1.0) { 00614 if( test >= 1.0e-10 )excs += dum*test; 00615 } else { 00616 excs += dum*test; 00617 } 00618 00619 if (ran < excs) goto outOfLoopAn; // --------------------------> 00620 } 00621 } 00622 } 00623 } 00624 inElastic = false; // quasi-elastic scattering. 00625 return; 00626 } 00627 outOfLoopAn: // <------------------------------------------------------------------ 00628 vecLen = 0; 00629 } 00630 } 00631 00632 nt = npos + nneg + nzero; 00633 while ( nt > 0) 00634 { 00635 G4double ran = G4UniformRand(); 00636 if ( ran < (G4double)npos/nt) 00637 { 00638 if( npos > 0 ) 00639 { pv[vecLen++] = PionPlus; 00640 npos--; 00641 } 00642 } 00643 else if ( ran < (G4double)(npos+nneg)/nt) 00644 { 00645 if( nneg > 0 ) 00646 { 00647 pv[vecLen++] = PionMinus; 00648 nneg--; 00649 } 00650 } 00651 else 00652 { 00653 if( nzero > 0 ) 00654 { 00655 pv[vecLen++] = PionZero; 00656 nzero--; 00657 } 00658 } 00659 nt = npos + nneg + nzero; 00660 } 00661 if (verboseLevel > 1) 00662 { 00663 G4cout << "Particles produced: " ; 00664 G4cout << pv[0].getName() << " " ; 00665 G4cout << pv[1].getName() << " " ; 00666 for (i=2; i < vecLen; i++) 00667 { 00668 G4cout << pv[i].getName() << " " ; 00669 } 00670 G4cout << G4endl; 00671 } 00672 return; 00673 }
| G4int G4HEAntiSigmaMinusInelastic::GetNumberOfSecondaries | ( | ) | [inline] |
Definition at line 76 of file G4HEAntiSigmaMinusInelastic.hh.
References vecLength.
00076 {return vecLength;}
| void G4HEAntiSigmaMinusInelastic::ModelDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 43 of file G4HEAntiSigmaMinusInelastic.cc.
00044 { 00045 outFile << "G4HEAntiSigmaMinusInelastic 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 G4HEAntiSigmaMinusInelastic.hh.
Referenced by ApplyYourself(), G4HEAntiSigmaMinusInelastic(), and GetNumberOfSecondaries().
1.4.7