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