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