#include <G4HEKaonMinusInelastic.hh>
Inheritance diagram for G4HEKaonMinusInelastic:
Public Member Functions | |
G4HEKaonMinusInelastic () | |
~G4HEKaonMinusInelastic () | |
virtual void | ModelDescription (std::ostream &) const |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
G4int | GetNumberOfSecondaries () |
void | FirstIntInCasKaonMinus (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle) |
Data Fields | |
G4int | vecLength |
Definition at line 52 of file G4HEKaonMinusInelastic.hh.
G4HEKaonMinusInelastic::G4HEKaonMinusInelastic | ( | ) | [inline] |
Definition at line 55 of file G4HEKaonMinusInelastic.hh.
References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, vecLength, and G4HEInelastic::verboseLevel.
00055 : G4HEInelastic("G4HEKaonMinusInelastic") 00056 { 00057 vecLength = 0; 00058 theMinEnergy = 20*CLHEP::GeV; 00059 theMaxEnergy = 10*CLHEP::TeV; 00060 MAXPART = 2048; 00061 verboseLevel = 0; 00062 G4cout << "WARNING: model G4HEKaonMinusInelastic is being deprecated and will\n" 00063 << "disappear in Geant4 version 10.0" << G4endl; 00064 }
G4HEKaonMinusInelastic::~G4HEKaonMinusInelastic | ( | ) | [inline] |
G4HadFinalState * G4HEKaonMinusInelastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | targetNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 58 of file G4HEKaonMinusInelastic.cc.
References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasKaonMinus(), 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 A = targetNucleus.GetA_asInt(); 00064 const G4double Z = targetNucleus.GetZ_asInt(); 00065 G4HEVector incidentParticle(aParticle); 00066 00067 G4double atomicNumber = Z; 00068 G4double atomicWeight = A; 00069 00070 G4int incidentCode = incidentParticle.getCode(); 00071 G4double incidentMass = incidentParticle.getMass(); 00072 G4double incidentTotalEnergy = incidentParticle.getEnergy(); 00073 00074 // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00075 // DHW 19 May 2011: variable set but not used 00076 00077 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass; 00078 00079 if (incidentKineticEnergy < 1.) 00080 G4cout << "GHEKaonMinusInelastic: incident energy < 1 GeV" << G4endl; 00081 00082 if (verboseLevel > 1) { 00083 G4cout << "G4HEKaonMinusInelastic::ApplyYourself" << G4endl; 00084 G4cout << "incident particle " << incidentParticle.getName() 00085 << "mass " << incidentMass 00086 << "kinetic energy " << incidentKineticEnergy 00087 << G4endl; 00088 G4cout << "target material with (A,Z) = (" 00089 << atomicWeight << "," << atomicNumber << ")" << G4endl; 00090 } 00091 00092 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 00093 atomicWeight, atomicNumber); 00094 if (verboseLevel > 1) 00095 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 00096 00097 incidentKineticEnergy -= inelasticity; 00098 00099 G4double excitationEnergyGNP = 0.; 00100 G4double excitationEnergyDTA = 0.; 00101 00102 G4double excitation = NuclearExcitation(incidentKineticEnergy, 00103 atomicWeight, atomicNumber, 00104 excitationEnergyGNP, 00105 excitationEnergyDTA); 00106 if (verboseLevel > 1) 00107 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 00108 << excitationEnergyDTA << G4endl; 00109 00110 incidentKineticEnergy -= excitation; 00111 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 00112 // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 00113 // *(incidentTotalEnergy+incidentMass)); 00114 // DHW 19 May 2011: variable set but not used 00115 00116 G4HEVector targetParticle; 00117 if (G4UniformRand() < atomicNumber/atomicWeight) { 00118 targetParticle.setDefinition("Proton"); 00119 } else { 00120 targetParticle.setDefinition("Neutron"); 00121 } 00122 00123 G4double targetMass = targetParticle.getMass(); 00124 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass 00125 + targetMass*targetMass 00126 + 2.0*targetMass*incidentTotalEnergy); 00127 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 00128 00129 G4bool inElastic = true; 00130 00131 vecLength = 0; 00132 00133 if (verboseLevel > 1) 00134 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 00135 << incidentCode << G4endl; 00136 00137 G4bool successful = false; 00138 00139 FirstIntInCasKaonMinus(inElastic, availableEnergy, pv, vecLength, 00140 incidentParticle, targetParticle); 00141 00142 if (verboseLevel > 1) 00143 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 00144 00145 if ((vecLength > 0) && (availableEnergy > 1.)) 00146 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 00147 pv, vecLength, 00148 incidentParticle, targetParticle); 00149 00150 HighEnergyCascading(successful, pv, vecLength, 00151 excitationEnergyGNP, excitationEnergyDTA, 00152 incidentParticle, targetParticle, 00153 atomicWeight, atomicNumber); 00154 if (!successful) 00155 HighEnergyClusterProduction(successful, pv, vecLength, 00156 excitationEnergyGNP, excitationEnergyDTA, 00157 incidentParticle, targetParticle, 00158 atomicWeight, atomicNumber); 00159 if (!successful) 00160 MediumEnergyCascading(successful, pv, vecLength, 00161 excitationEnergyGNP, excitationEnergyDTA, 00162 incidentParticle, targetParticle, 00163 atomicWeight, atomicNumber); 00164 00165 if (!successful) 00166 MediumEnergyClusterProduction(successful, pv, vecLength, 00167 excitationEnergyGNP, excitationEnergyDTA, 00168 incidentParticle, targetParticle, 00169 atomicWeight, atomicNumber); 00170 if (!successful) 00171 QuasiElasticScattering(successful, pv, vecLength, 00172 excitationEnergyGNP, excitationEnergyDTA, 00173 incidentParticle, targetParticle, 00174 atomicWeight, atomicNumber); 00175 if (!successful) 00176 ElasticScattering(successful, pv, vecLength, 00177 incidentParticle, 00178 atomicWeight, atomicNumber); 00179 if (!successful) 00180 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 00181 << G4endl; 00182 00183 FillParticleChange(pv, vecLength); 00184 delete [] pv; 00185 theParticleChange.SetStatusChange(stopAndKill); 00186 return &theParticleChange; 00187 }
void G4HEKaonMinusInelastic::FirstIntInCasKaonMinus | ( | G4bool & | inElastic, | |
const G4double | availableEnergy, | |||
G4HEVector | pv[], | |||
G4int & | vecLen, | |||
const G4HEVector & | incidentParticle, | |||
const G4HEVector & | targetParticle | |||
) |
Definition at line 191 of file G4HEKaonMinusInelastic.cc.
References G4HEInelastic::Amax(), G4HEInelastic::Amin(), G4HEInelastic::AntiKaonZero, G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), G4HEInelastic::Imax(), G4HEInelastic::Imin(), G4HEInelastic::KaonMinus, G4HEInelastic::KaonZero, G4HEInelastic::Lambda, CLHEP::detail::n, G4HEInelastic::Neutron, neutronCode, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, G4HEInelastic::SigmaMinus, G4HEInelastic::SigmaPlus, G4HEInelastic::SigmaZero, and G4HEInelastic::verboseLevel.
Referenced by ApplyYourself().
00205 { 00206 static const G4double expxu = 82.; // upper bound for arg. of exp 00207 static const G4double expxl = -expxu; // lower bound for arg. of exp 00208 00209 static const G4double protb = 0.7; 00210 static const G4double neutb = 0.7; 00211 static const G4double c = 1.25; 00212 00213 static const G4int numMul = 1200; 00214 static const G4int numSec = 60; 00215 00216 G4int neutronCode = Neutron.getCode(); 00217 G4int protonCode = Proton.getCode(); 00218 G4int kaonMinusCode = KaonMinus.getCode(); 00219 G4int kaonZeroCode = KaonZero.getCode(); 00220 G4int antiKaonZeroCode = AntiKaonZero.getCode(); 00221 00222 G4int targetCode = targetParticle.getCode(); 00223 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00224 00225 static G4bool first = true; 00226 static G4double protmul[numMul], protnorm[numSec]; // proton constants 00227 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 00228 00229 // misc. local variables 00230 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0 00231 00232 G4int i, counter, nt, npos, nneg, nzero; 00233 00234 if (first) { 00235 // compute normalization constants, this will only be done once 00236 first = false; 00237 for (i = 0; i < numMul; i++) protmul[i] = 0.0; 00238 for (i = 0; i < numSec; i++) protnorm[i] = 0.0; 00239 counter = -1; 00240 for (npos = 0; npos < (numSec/3); npos++) { 00241 for( nneg=Imax(0,npos-1); nneg<=npos+1; nneg++ ) 00242 { 00243 for( nzero=0; nzero<numSec/3; nzero++ ) 00244 { 00245 if( ++counter < numMul ) 00246 { 00247 nt = npos+nneg+nzero; 00248 if( (nt>0) && (nt<=numSec) ) 00249 { 00250 protmul[counter] = 00251 pmltpc(npos,nneg,nzero,nt,protb,c) ; 00252 protnorm[nt-1] += protmul[counter]; 00253 } 00254 } 00255 } 00256 } 00257 } 00258 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 00259 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 00260 counter = -1; 00261 for (npos = 0; npos < numSec/3; npos++) { 00262 for (nneg = npos; nneg <= (npos+2); nneg++) { 00263 for( nzero=0; nzero<numSec/3; nzero++ ) 00264 { 00265 if( ++counter < numMul ) 00266 { 00267 nt = npos+nneg+nzero; 00268 if( (nt>0) && (nt<=numSec) ) 00269 { 00270 neutmul[counter] = 00271 pmltpc(npos,nneg,nzero,nt,neutb,c); 00272 neutnorm[nt-1] += neutmul[counter]; 00273 } 00274 } 00275 } 00276 } 00277 } 00278 00279 for (i = 0; i < numSec; i++) { 00280 if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i]; 00281 if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i]; 00282 } 00283 } // end of initialization 00284 00285 pv[0] = incidentParticle; // initialize the first two places 00286 pv[1] = targetParticle; // the same as beam and target 00287 vecLen = 2; 00288 00289 if (!inElastic || (availableEnergy <= PionPlus.getMass())) return; 00290 00291 // inelastic scattering 00292 npos = 0, nneg = 0, nzero = 0; 00293 G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15}; 00294 G4int iplab = G4int( incidentTotalMomentum*5.); 00295 if ( (iplab < 10) && (G4UniformRand() < cech[iplab])) { 00296 G4int ipl = Imin(19, G4int( incidentTotalMomentum*5.)); 00297 G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45, 00298 0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33}; 00299 if (G4UniformRand() < cnk0[ipl]) { 00300 if (targetCode == protonCode) { 00301 pv[0] = AntiKaonZero; 00302 pv[1] = Neutron; 00303 return; 00304 } else { 00305 return; 00306 } 00307 } 00308 G4double ran = G4UniformRand(); 00309 if (targetCode == protonCode) { // target is a proton 00310 if( ran < 0.25 ) 00311 { 00312 pv[0] = PionMinus; 00313 pv[1] = SigmaPlus; 00314 } 00315 else if (ran < 0.50) 00316 { 00317 pv[0] = PionZero; 00318 pv[1] = SigmaZero; 00319 } 00320 else if (ran < 0.75) 00321 { 00322 pv[0] = PionPlus; 00323 pv[1] = SigmaMinus; 00324 } 00325 else 00326 { 00327 pv[0] = PionZero; 00328 pv[1] = Lambda; 00329 } 00330 } else { // target is a neutron 00331 if( ran < 0.25 ) 00332 { 00333 } 00334 else if (ran < 0.50) 00335 { 00336 pv[0] = PionMinus; 00337 pv[1] = SigmaZero; 00338 } 00339 else if (ran < 0.75) 00340 { 00341 pv[0] = PionZero; 00342 pv[1] = SigmaMinus; 00343 } 00344 else 00345 { 00346 pv[0] = PionMinus; 00347 pv[1] = Lambda; 00348 } 00349 } 00350 return; 00351 } else { 00352 // number of total particles vs. centre of mass Energy - 2*proton mass 00353 G4double aleab = std::log(availableEnergy); 00354 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00355 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00356 00357 // normalization constant for kno-distribution. 00358 // calculate first the sum of all constants, check for numerical problems. 00359 G4double test, dum, anpn = 0.0; 00360 00361 for (nt=1; nt<=numSec; nt++) { 00362 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00363 dum = pi*nt/(2.0*n*n); 00364 if (std::fabs(dum) < 1.0) { 00365 if( test >= 1.0e-10 )anpn += dum*test; 00366 } else { 00367 anpn += dum*test; 00368 } 00369 } 00370 00371 G4double ran = G4UniformRand(); 00372 G4double excs = 0.0; 00373 if( targetCode == protonCode ) 00374 { 00375 counter = -1; 00376 for( npos=0; npos<numSec/3; npos++ ) 00377 { 00378 for( nneg=Imax(0,npos-1); nneg<=npos+1; nneg++ ) 00379 { 00380 for (nzero=0; nzero<numSec/3; nzero++) { 00381 if (++counter < numMul) { 00382 nt = npos+nneg+nzero; 00383 if ( (nt>0) && (nt<=numSec) ) { 00384 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00385 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 00386 if (std::fabs(dum) < 1.0) { 00387 if( test >= 1.0e-10 )excs += dum*test; 00388 } else { 00389 excs += dum*test; 00390 } 00391 00392 if (ran < excs) goto outOfLoop; //-----------------------> 00393 } 00394 } 00395 } 00396 } 00397 } 00398 00399 // 3 previous loops continued to the end 00400 inElastic = false; // quasi-elastic scattering 00401 return; 00402 } 00403 else 00404 { // target must be a neutron 00405 counter = -1; 00406 for( npos=0; npos<numSec/3; npos++ ) 00407 { 00408 for( nneg=npos; nneg<=(npos+2); nneg++ ) 00409 { 00410 for (nzero=0; nzero<numSec/3; nzero++) { 00411 if (++counter < numMul) { 00412 nt = npos+nneg+nzero; 00413 if ( (nt>=1) && (nt<=numSec) ) { 00414 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00415 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 00416 if (std::fabs(dum) < 1.0) { 00417 if( test >= 1.0e-10 )excs += dum*test; 00418 } else { 00419 excs += dum*test; 00420 } 00421 if (ran < excs) goto outOfLoop; // --------------------------> 00422 } 00423 } 00424 } 00425 } 00426 } 00427 // 3 previous loops continued to the end 00428 inElastic = false; // quasi-elastic scattering. 00429 return; 00430 } 00431 } 00432 outOfLoop: // <--------------------------------------------- 00433 00434 if( targetCode == protonCode) 00435 { 00436 if( npos == (1+nneg)) 00437 { 00438 pv[1] = Neutron; 00439 } 00440 else if (npos == nneg) 00441 { 00442 if( G4UniformRand() < 0.75) 00443 { 00444 } 00445 else 00446 { 00447 pv[0] = AntiKaonZero; 00448 pv[1] = Neutron; 00449 } 00450 } 00451 else 00452 { 00453 pv[0] = AntiKaonZero; 00454 } 00455 } 00456 else 00457 { 00458 if( npos == (nneg-1)) 00459 { 00460 if( G4UniformRand() < 0.5) 00461 { 00462 pv[1] = Proton; 00463 } 00464 else 00465 { 00466 pv[0] = AntiKaonZero; 00467 } 00468 } 00469 else if ( npos == nneg) 00470 { 00471 } 00472 else 00473 { 00474 pv[0] = AntiKaonZero; 00475 pv[1] = Proton; 00476 } 00477 } 00478 00479 if( G4UniformRand() < 0.5 ) 00480 { 00481 if( ( (pv[0].getCode() == kaonMinusCode) 00482 && (pv[1].getCode() == neutronCode) ) 00483 || ( (pv[0].getCode() == kaonZeroCode) 00484 && (pv[1].getCode() == protonCode) ) 00485 || ( (pv[0].getCode() == antiKaonZeroCode) 00486 && (pv[1].getCode() == protonCode) ) ) 00487 { 00488 G4double ran = G4UniformRand(); 00489 if( pv[1].getCode() == protonCode) 00490 { 00491 if(ran < 0.68) 00492 { 00493 pv[0] = PionPlus; 00494 pv[1] = Lambda; 00495 } 00496 else if (ran < 0.84) 00497 { 00498 pv[0] = PionZero; 00499 pv[1] = SigmaPlus; 00500 } 00501 else 00502 { 00503 pv[0] = PionPlus; 00504 pv[1] = SigmaZero; 00505 } 00506 } 00507 else 00508 { 00509 if(ran < 0.68) 00510 { 00511 pv[0] = PionMinus; 00512 pv[1] = Lambda; 00513 } 00514 else if (ran < 0.84) 00515 { 00516 pv[0] = PionMinus; 00517 pv[1] = SigmaZero; 00518 } 00519 else 00520 { 00521 pv[0] = PionZero; 00522 pv[1] = SigmaMinus; 00523 } 00524 } 00525 } 00526 else 00527 { 00528 G4double ran = G4UniformRand(); 00529 if (ran < 0.67) 00530 { 00531 pv[0] = PionZero; 00532 pv[1] = Lambda; 00533 } 00534 else if (ran < 0.78) 00535 { 00536 pv[0] = PionMinus; 00537 pv[1] = SigmaPlus; 00538 } 00539 else if (ran < 0.89) 00540 { 00541 pv[0] = PionZero; 00542 pv[1] = SigmaZero; 00543 } 00544 else 00545 { 00546 pv[0] = PionPlus; 00547 pv[1] = SigmaMinus; 00548 } 00549 } 00550 } 00551 00552 00553 nt = npos + nneg + nzero; 00554 while ( nt > 0) 00555 { 00556 G4double ran = G4UniformRand(); 00557 if ( ran < (G4double)npos/nt) 00558 { 00559 if( npos > 0 ) 00560 { pv[vecLen++] = PionPlus; 00561 npos--; 00562 } 00563 } 00564 else if ( ran < (G4double)(npos+nneg)/nt) 00565 { 00566 if( nneg > 0 ) 00567 { 00568 pv[vecLen++] = PionMinus; 00569 nneg--; 00570 } 00571 } 00572 else 00573 { 00574 if( nzero > 0 ) 00575 { 00576 pv[vecLen++] = PionZero; 00577 nzero--; 00578 } 00579 } 00580 nt = npos + nneg + nzero; 00581 } 00582 if (verboseLevel > 1) 00583 { 00584 G4cout << "Particles produced: " ; 00585 G4cout << pv[0].getName() << " " ; 00586 G4cout << pv[1].getName() << " " ; 00587 for (i=2; i < vecLen; i++) 00588 { 00589 G4cout << pv[i].getName() << " " ; 00590 } 00591 G4cout << G4endl; 00592 } 00593 return; 00594 }
G4int G4HEKaonMinusInelastic::GetNumberOfSecondaries | ( | ) | [inline] |
Definition at line 75 of file G4HEKaonMinusInelastic.hh.
References vecLength.
00075 {return vecLength;}
void G4HEKaonMinusInelastic::ModelDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 43 of file G4HEKaonMinusInelastic.cc.
00044 { 00045 outFile << "G4HEKaonMinusInelastic is one of the High Energy\n" 00046 << "Parameterized (HEP) models used to implement inelastic\n" 00047 << "K- 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 K- with initial energies\n" 00053 << "above 20 GeV.\n"; 00054 }
Definition at line 70 of file G4HEKaonMinusInelastic.hh.
Referenced by ApplyYourself(), G4HEKaonMinusInelastic(), and GetNumberOfSecondaries().