#include <G4HEKaonPlusInelastic.hh>
Inheritance diagram for G4HEKaonPlusInelastic:
Public Member Functions | |
| G4HEKaonPlusInelastic () | |
| ~G4HEKaonPlusInelastic () | |
| virtual void | ModelDescription (std::ostream &) const |
| G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
| G4int | GetNumberOfSecondaries () |
| void | FirstIntInCasKaonPlus (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 G4HEKaonPlusInelastic.hh.
| G4HEKaonPlusInelastic::G4HEKaonPlusInelastic | ( | ) | [inline] |
Definition at line 55 of file G4HEKaonPlusInelastic.hh.
References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, vecLength, and G4HEInelastic::verboseLevel.
00055 : G4HEInelastic("G4HEKaonPlusInelastic") 00056 { 00057 vecLength = 0; 00058 theMinEnergy = 20*CLHEP::GeV; 00059 theMaxEnergy = 10*CLHEP::TeV; 00060 MAXPART = 2048; 00061 verboseLevel = 0; 00062 G4cout << "WARNING: model G4HEKaonPlusInelastic is being deprecated and will\n" 00063 << "disappear in Geant4 version 10.0" << G4endl; 00064 }
| G4HEKaonPlusInelastic::~G4HEKaonPlusInelastic | ( | ) | [inline] |
| G4HadFinalState * G4HEKaonPlusInelastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
| G4Nucleus & | targetNucleus | |||
| ) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 57 of file G4HEKaonPlusInelastic.cc.
References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasKaonPlus(), 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 << "GHEKaonPlusInelastic: incident energy < 1 GeV" << G4endl; 00080 00081 if (verboseLevel > 1) { 00082 G4cout << "G4HEKaonPlusInelastic::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 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 vecLength = 0; 00131 00132 if (verboseLevel > 1) 00133 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 00134 << incidentCode << G4endl; 00135 00136 G4bool successful = false; 00137 00138 FirstIntInCasKaonPlus(inElastic, availableEnergy, pv, vecLength, 00139 incidentParticle, targetParticle, atomicWeight); 00140 00141 if (verboseLevel > 1) 00142 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 00143 00144 if ((vecLength > 0) && (availableEnergy > 1.)) 00145 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 00146 pv, vecLength, 00147 incidentParticle, targetParticle); 00148 00149 HighEnergyCascading(successful, pv, vecLength, 00150 excitationEnergyGNP, excitationEnergyDTA, 00151 incidentParticle, targetParticle, 00152 atomicWeight, atomicNumber); 00153 if (!successful) 00154 HighEnergyClusterProduction(successful, pv, vecLength, 00155 excitationEnergyGNP, excitationEnergyDTA, 00156 incidentParticle, targetParticle, 00157 atomicWeight, atomicNumber); 00158 if (!successful) 00159 MediumEnergyCascading(successful, pv, vecLength, 00160 excitationEnergyGNP, excitationEnergyDTA, 00161 incidentParticle, targetParticle, 00162 atomicWeight, atomicNumber); 00163 00164 if (!successful) 00165 MediumEnergyClusterProduction(successful, pv, vecLength, 00166 excitationEnergyGNP, excitationEnergyDTA, 00167 incidentParticle, targetParticle, 00168 atomicWeight, atomicNumber); 00169 if (!successful) 00170 QuasiElasticScattering(successful, pv, vecLength, 00171 excitationEnergyGNP, excitationEnergyDTA, 00172 incidentParticle, targetParticle, 00173 atomicWeight, atomicNumber); 00174 if (!successful) 00175 ElasticScattering(successful, pv, vecLength, 00176 incidentParticle, 00177 atomicWeight, atomicNumber); 00178 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 G4HEKaonPlusInelastic::FirstIntInCasKaonPlus | ( | G4bool & | inElastic, | |
| const G4double | availableEnergy, | |||
| G4HEVector | pv[], | |||
| G4int & | vecLen, | |||
| const G4HEVector & | incidentParticle, | |||
| const G4HEVector & | targetParticle, | |||
| const G4double | atomicWeight | |||
| ) |
Definition at line 191 of file G4HEKaonPlusInelastic.cc.
References G4HEInelastic::Amax(), G4HEInelastic::Amin(), G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), G4HEInelastic::Imax(), G4HEInelastic::KaonZero, CLHEP::detail::n, G4HEInelastic::Neutron, neutronCode, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, sqr(), and G4HEInelastic::verboseLevel.
Referenced by ApplyYourself().
00206 { 00207 static const G4double expxu = 82.; // upper bound for arg. of exp 00208 static const G4double expxl = -expxu; // lower bound for arg. of exp 00209 00210 static const G4double protb = 0.7; 00211 static const G4double neutb = 0.7; 00212 static const G4double c = 1.25; 00213 00214 static const G4int numMul = 1200; 00215 static const G4int numSec = 60; 00216 00217 G4int neutronCode = Neutron.getCode(); 00218 G4int protonCode = Proton.getCode(); 00219 00220 G4int targetCode = targetParticle.getCode(); 00221 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00222 00223 static G4bool first = true; 00224 static G4double protmul[numMul], protnorm[numSec]; // proton constants 00225 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 00226 00227 // misc. local variables 00228 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0 00229 00230 G4int i, counter, nt, npos, nneg, nzero; 00231 00232 if( first ) 00233 { // compute normalization constants, this will only be done once 00234 first = false; 00235 for( i=0; i<numMul; i++ )protmul[i] = 0.0; 00236 for( i=0; i<numSec; i++ )protnorm[i] = 0.0; 00237 counter = -1; 00238 for( npos=0; npos<(numSec/3); npos++ ) 00239 { 00240 for( nneg=Imax(0,npos-2); nneg<=npos; nneg++ ) 00241 { 00242 for( nzero=0; nzero<numSec/3; nzero++ ) 00243 { 00244 if( ++counter < numMul ) 00245 { 00246 nt = npos+nneg+nzero; 00247 if( (nt>0) && (nt<=numSec) ) 00248 { 00249 protmul[counter] = 00250 pmltpc(npos,nneg,nzero,nt,protb,c) ; 00251 protnorm[nt-1] += protmul[counter]; 00252 } 00253 } 00254 } 00255 } 00256 } 00257 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 00258 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 00259 counter = -1; 00260 for( npos=0; npos<numSec/3; npos++ ) 00261 { 00262 for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ ) 00263 { 00264 for( nzero=0; nzero<numSec/3; nzero++ ) 00265 { 00266 if( ++counter < numMul ) 00267 { 00268 nt = npos+nneg+nzero; 00269 if( (nt>0) && (nt<=numSec) ) 00270 { 00271 neutmul[counter] = 00272 pmltpc(npos,nneg,nzero,nt,neutb,c); 00273 neutnorm[nt-1] += neutmul[counter]; 00274 } 00275 } 00276 } 00277 } 00278 } 00279 for( i=0; i<numSec; i++ ) 00280 { 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 00287 // initialize the first two places 00288 // the same as beam and target 00289 pv[0] = incidentParticle; 00290 pv[1] = targetParticle; 00291 vecLen = 2; 00292 00293 if( !inElastic ) 00294 { // quasi-elastic scattering, no pions produced 00295 if( targetCode == neutronCode ) 00296 { 00297 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07}; 00298 G4int iplab = G4int( Amin( 9.0, incidentTotalMomentum*5. ) ); 00299 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 00300 { // charge exchange K+ n -> K0 p 00301 pv[0] = KaonZero; 00302 pv[1] = Proton; 00303 } 00304 } 00305 return; 00306 } 00307 else if (availableEnergy <= PionPlus.getMass()) 00308 return; 00309 00310 // inelastic scattering 00311 00312 npos = 0, nneg = 0, nzero = 0; 00313 G4double eab = availableEnergy; 00314 G4int ieab = G4int( eab*5.0 ); 00315 00316 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98}; 00317 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) ) 00318 { 00319 // suppress high multiplicity events at low momentum 00320 // only one additional pion will be produced 00321 G4double w0, wp, wm, wt, ran; 00322 if( targetCode == protonCode ) // target is a proton 00323 { 00324 w0 = - sqr(1.+protb)/(2.*c*c); 00325 wp = w0 = std::exp(w0); 00326 wp *= 2.; 00327 if( G4UniformRand() < w0/(w0+wp) ) 00328 { npos = 0; nneg = 0; nzero = 1; } 00329 else 00330 { npos = 1; nneg = 0; nzero = 0; } 00331 } 00332 else 00333 { // target is a neutron 00334 w0 = -sqr(1.+neutb)/(2.*c*c); 00335 wp = w0 = std::exp(w0); 00336 wm = -sqr(-1.+neutb)/(2.*c*c); 00337 wm = std::exp(wm); 00338 wt = w0+wp+wm; 00339 wp = w0+wp; 00340 ran = G4UniformRand(); 00341 if( ran < w0/wt) 00342 { npos = 0; nneg = 0; nzero = 1; } 00343 else if( ran < wp/wt) 00344 { npos = 1; nneg = 0; nzero = 0; } 00345 else 00346 { npos = 0; nneg = 1; nzero = 0; } 00347 } 00348 } 00349 else 00350 { 00351 // number of total particles vs. centre of mass Energy - 2*proton mass 00352 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-2); nneg<=npos; 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 if (ran < excs) goto outOfLoop; //-----------------------> 00392 } 00393 } 00394 } 00395 } 00396 } 00397 00398 // 3 previous loops continued to the end 00399 inElastic = false; // quasi-elastic scattering 00400 return; 00401 } 00402 else 00403 { // target must be a neutron 00404 counter = -1; 00405 for( npos=0; npos<numSec/3; npos++ ) 00406 { 00407 for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ ) 00408 { 00409 for (nzero=0; nzero<numSec/3; nzero++) { 00410 if (++counter < numMul) { 00411 nt = npos+nneg+nzero; 00412 if ( (nt>=1) && (nt<=numSec) ) { 00413 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00414 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 00415 if (std::fabs(dum) < 1.0) { 00416 if( test >= 1.0e-10 )excs += dum*test; 00417 } else { 00418 excs += dum*test; 00419 } 00420 if (ran < excs) goto outOfLoop; // --------------------------> 00421 } 00422 } 00423 } 00424 } 00425 } 00426 // 3 previous loops continued to the end 00427 inElastic = false; // quasi-elastic scattering. 00428 return; 00429 } 00430 } 00431 outOfLoop: // <------------------------------------------------------------------------ 00432 00433 if( targetCode == protonCode) 00434 { 00435 if( npos == nneg) 00436 { 00437 } 00438 else if (npos == (1+nneg)) 00439 { 00440 if( G4UniformRand() < 0.5) 00441 { 00442 pv[1] = Neutron; 00443 } 00444 else 00445 { 00446 pv[0] = KaonZero; 00447 } 00448 } 00449 else 00450 { 00451 pv[0] = KaonZero; 00452 pv[1] = Neutron; 00453 } 00454 } 00455 else 00456 { 00457 if( npos == nneg) 00458 { 00459 if( G4UniformRand() < 0.25) 00460 { 00461 pv[0] = KaonZero; 00462 pv[1] = Proton; 00463 } 00464 else 00465 { 00466 } 00467 } 00468 else if ( npos == (1+nneg)) 00469 { 00470 pv[0] = KaonZero; 00471 } 00472 else 00473 { 00474 pv[1] = Proton; 00475 } 00476 } 00477 00478 00479 nt = npos + nneg + nzero; 00480 while ( nt > 0) 00481 { 00482 G4double ran = G4UniformRand(); 00483 if ( ran < (G4double)npos/nt) 00484 { 00485 if( npos > 0 ) 00486 { pv[vecLen++] = PionPlus; 00487 npos--; 00488 } 00489 } 00490 else if ( ran < (G4double)(npos+nneg)/nt) 00491 { 00492 if( nneg > 0 ) 00493 { 00494 pv[vecLen++] = PionMinus; 00495 nneg--; 00496 } 00497 } 00498 else 00499 { 00500 if( nzero > 0 ) 00501 { 00502 pv[vecLen++] = PionZero; 00503 nzero--; 00504 } 00505 } 00506 nt = npos + nneg + nzero; 00507 } 00508 if (verboseLevel > 1) 00509 { 00510 G4cout << "Particles produced: " ; 00511 G4cout << pv[0].getName() << " " ; 00512 G4cout << pv[1].getName() << " " ; 00513 for (i=2; i < vecLen; i++) 00514 { 00515 G4cout << pv[i].getName() << " " ; 00516 } 00517 G4cout << G4endl; 00518 } 00519 return; 00520 }
| G4int G4HEKaonPlusInelastic::GetNumberOfSecondaries | ( | ) | [inline] |
Definition at line 75 of file G4HEKaonPlusInelastic.hh.
References vecLength.
00075 {return vecLength;}
| void G4HEKaonPlusInelastic::ModelDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 42 of file G4HEKaonPlusInelastic.cc.
00043 { 00044 outFile << "G4HEKaonPlusInelastic is one of the High Energy\n" 00045 << "Parameterized (HEP) models used to implement inelastic\n" 00046 << "K+ 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 K+ with initial energies\n" 00052 << "above 20 GeV.\n"; 00053 }
Definition at line 70 of file G4HEKaonPlusInelastic.hh.
Referenced by ApplyYourself(), G4HEKaonPlusInelastic(), and GetNumberOfSecondaries().
1.4.7