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