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