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