#include <G4HEOmegaMinusInelastic.hh>
Inheritance diagram for G4HEOmegaMinusInelastic:
Public Member Functions | |
G4HEOmegaMinusInelastic () | |
~G4HEOmegaMinusInelastic () | |
virtual void | ModelDescription (std::ostream &) const |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
G4int | GetNumberOfSecondaries () |
void | FirstIntInCasOmegaMinus (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 G4HEOmegaMinusInelastic.hh.
G4HEOmegaMinusInelastic::G4HEOmegaMinusInelastic | ( | ) | [inline] |
Definition at line 55 of file G4HEOmegaMinusInelastic.hh.
References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, vecLength, and G4HEInelastic::verboseLevel.
00055 : G4HEInelastic("G4HEOmegaMinusInelastic") 00056 { 00057 vecLength = 0; 00058 theMinEnergy = 20*CLHEP::GeV; 00059 theMaxEnergy = 10*CLHEP::TeV; 00060 MAXPART = 2048; 00061 verboseLevel = 0; 00062 G4cout << "WARNING: model G4HEOmegaMinusInelastic is being deprecated and will\n" 00063 << "disappear in Geant4 version 10.0" << G4endl; 00064 }
G4HEOmegaMinusInelastic::~G4HEOmegaMinusInelastic | ( | ) | [inline] |
G4HadFinalState * G4HEOmegaMinusInelastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | targetNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 57 of file G4HEOmegaMinusInelastic.cc.
References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasOmegaMinus(), 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 << "GHEOmegaMinusInelastic: incident energy < 1 GeV" << G4endl; 00080 00081 if (verboseLevel > 1) { 00082 G4cout << "G4HEOmegaMinusInelastic::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 if (G4UniformRand() < atomicNumber/atomicWeight) { 00117 targetParticle.setDefinition("Proton"); 00118 } else { 00119 targetParticle.setDefinition("Neutron"); 00120 } 00121 00122 G4double targetMass = targetParticle.getMass(); 00123 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass 00124 + targetMass*targetMass 00125 + 2.0*targetMass*incidentTotalEnergy); 00126 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 00127 00128 G4bool inElastic = true; 00129 vecLength = 0; 00130 00131 if (verboseLevel > 1) 00132 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 00133 << incidentCode << G4endl; 00134 00135 G4bool successful = false; 00136 00137 FirstIntInCasOmegaMinus(inElastic, availableEnergy, pv, vecLength, 00138 incidentParticle, targetParticle, atomicWeight); 00139 00140 if (verboseLevel > 1) 00141 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 00142 00143 if ((vecLength > 0) && (availableEnergy > 1.)) 00144 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 00145 pv, vecLength, 00146 incidentParticle, targetParticle); 00147 00148 HighEnergyCascading(successful, pv, vecLength, 00149 excitationEnergyGNP, excitationEnergyDTA, 00150 incidentParticle, targetParticle, 00151 atomicWeight, atomicNumber); 00152 if (!successful) 00153 HighEnergyClusterProduction(successful, pv, vecLength, 00154 excitationEnergyGNP, excitationEnergyDTA, 00155 incidentParticle, targetParticle, 00156 atomicWeight, atomicNumber); 00157 if (!successful) 00158 MediumEnergyCascading(successful, pv, vecLength, 00159 excitationEnergyGNP, excitationEnergyDTA, 00160 incidentParticle, targetParticle, 00161 atomicWeight, atomicNumber); 00162 00163 if (!successful) 00164 MediumEnergyClusterProduction(successful, pv, vecLength, 00165 excitationEnergyGNP, excitationEnergyDTA, 00166 incidentParticle, targetParticle, 00167 atomicWeight, atomicNumber); 00168 if (!successful) 00169 QuasiElasticScattering(successful, pv, vecLength, 00170 excitationEnergyGNP, excitationEnergyDTA, 00171 incidentParticle, targetParticle, 00172 atomicWeight, atomicNumber); 00173 if (!successful) 00174 ElasticScattering(successful, pv, vecLength, 00175 incidentParticle, 00176 atomicWeight, atomicNumber); 00177 00178 if (!successful) 00179 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 00180 << G4endl; 00181 00182 FillParticleChange(pv, vecLength); 00183 delete [] pv; 00184 theParticleChange.SetStatusChange(stopAndKill); 00185 return &theParticleChange; 00186 }
void G4HEOmegaMinusInelastic::FirstIntInCasOmegaMinus | ( | G4bool & | inElastic, | |
const G4double | availableEnergy, | |||
G4HEVector | pv[], | |||
G4int & | vecLen, | |||
const G4HEVector & | incidentParticle, | |||
const G4HEVector & | targetParticle, | |||
const G4double | atomicWeight | |||
) |
Definition at line 190 of file G4HEOmegaMinusInelastic.cc.
References G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), G4HEInelastic::Lambda, CLHEP::detail::n, G4HEInelastic::Neutron, G4HEInelastic::OmegaMinus, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, G4HEInelastic::SigmaMinus, G4HEInelastic::SigmaZero, G4HEInelastic::verboseLevel, G4HEInelastic::XiMinus, and G4HEInelastic::XiZero.
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 protonCode = Proton.getCode(); 00217 00218 G4int targetCode = targetParticle.getCode(); 00219 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00220 00221 static G4bool first = true; 00222 static G4double protmul[numMul], protnorm[numSec]; // proton constants 00223 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 00224 00225 // misc. local variables 00226 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0 00227 00228 G4int i, counter, nt, npos, nneg, nzero; 00229 00230 if (first) { 00231 // compute normalization constants, this will only be done once 00232 first = false; 00233 for( i=0; i<numMul; i++ )protmul[i] = 0.0; 00234 for( i=0; i<numSec; i++ )protnorm[i] = 0.0; 00235 counter = -1; 00236 for( npos=0; npos<(numSec/3); npos++ ) 00237 { 00238 for( nneg=std::max(0,npos-1); nneg<=npos; nneg++ ) 00239 { 00240 for( nzero=0; nzero<numSec/3; nzero++ ) 00241 { 00242 if( ++counter < numMul ) 00243 { 00244 nt = npos+nneg+nzero; 00245 if( (nt>0) && (nt<=numSec) ) 00246 { 00247 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c); 00248 protnorm[nt-1] += protmul[counter]; 00249 } 00250 } 00251 } 00252 } 00253 } 00254 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 00255 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 00256 counter = -1; 00257 for( npos=0; npos<numSec/3; npos++ ) 00258 { 00259 for( nneg=npos; nneg<=(npos+1); nneg++ ) 00260 { 00261 for( nzero=0; nzero<numSec/3; nzero++ ) 00262 { 00263 if( ++counter < numMul ) 00264 { 00265 nt = npos+nneg+nzero; 00266 if( (nt>0) && (nt<=numSec) ) 00267 { 00268 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c); 00269 neutnorm[nt-1] += neutmul[counter]; 00270 } 00271 } 00272 } 00273 } 00274 } 00275 for( i=0; i<numSec; i++ ) 00276 { 00277 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 00278 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 00279 } 00280 } // end of initialization 00281 00282 00283 // initialize the first two places 00284 // the same as beam and target 00285 pv[0] = incidentParticle; 00286 pv[1] = targetParticle; 00287 vecLen = 2; 00288 00289 if( !inElastic ) 00290 { // quasi-elastic scattering, no pions produced 00291 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.}; 00292 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) ); 00293 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 00294 { 00295 G4double ran = G4UniformRand(); 00296 if( targetCode == protonCode) 00297 { 00298 if (ran < 0.2) 00299 { 00300 pv[0] = XiZero; 00301 pv[1] = SigmaZero; 00302 } 00303 else if (ran < 0.4) 00304 { 00305 pv[0] = SigmaZero; 00306 pv[1] = XiZero; 00307 } 00308 else if (ran < 0.6) 00309 { 00310 pv[0] = XiZero; 00311 pv[1] = Lambda; 00312 } 00313 else if (ran < 0.8) 00314 { 00315 pv[0] = Lambda; 00316 pv[1] = XiZero; 00317 } 00318 else 00319 { 00320 pv[0] = Proton; 00321 pv[1] = OmegaMinus; 00322 } 00323 } 00324 else 00325 { 00326 if (ran < 0.2) 00327 { 00328 pv[0] = Neutron; 00329 pv[1] = OmegaMinus; 00330 } 00331 else if (ran < 0.4) 00332 { 00333 pv[0] = XiZero; 00334 pv[1] = SigmaMinus; 00335 } 00336 else if (ran < 0.6) 00337 { 00338 pv[0] = SigmaMinus; 00339 pv[1] = XiZero; 00340 } 00341 else if (ran < 0.8) 00342 { 00343 pv[0] = XiMinus; 00344 pv[1] = Lambda; 00345 } 00346 else 00347 { 00348 pv[0] = Lambda; 00349 pv[1] = XiMinus; 00350 } 00351 } 00352 } 00353 return; 00354 } 00355 else if (availableEnergy <= PionPlus.getMass()) 00356 return; 00357 00358 // inelastic scattering 00359 npos = 0; nneg = 0; nzero = 0; 00360 00361 // number of total particles vs. centre of mass Energy - 2*proton mass 00362 G4double aleab = std::log(availableEnergy); 00363 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00364 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00365 00366 // normalization constant for kno-distribution. 00367 // calculate first the sum of all constants, check for numerical problems. 00368 G4double test, dum, anpn = 0.0; 00369 00370 for (nt=1; nt<=numSec; nt++) { 00371 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00372 dum = pi*nt/(2.0*n*n); 00373 if (std::fabs(dum) < 1.0) { 00374 if( test >= 1.0e-10 )anpn += dum*test; 00375 } else { 00376 anpn += dum*test; 00377 } 00378 } 00379 00380 G4double ran = G4UniformRand(); 00381 G4double excs = 0.0; 00382 if (targetCode == protonCode) { 00383 counter = -1; 00384 for (npos = 0; npos < numSec/3; npos++) { 00385 for( nneg=std::max(0,npos-1); nneg<=npos; nneg++ ) 00386 { 00387 for( nzero=0; nzero<numSec/3; nzero++ ) 00388 { 00389 if( ++counter < numMul ) 00390 { 00391 nt = npos+nneg+nzero; 00392 if ( (nt>0) && (nt<=numSec) ) { 00393 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00394 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 00395 if (std::fabs(dum) < 1.0) { 00396 if( test >= 1.0e-10 )excs += dum*test; 00397 } else { 00398 excs += dum*test; 00399 } 00400 if (ran < excs) goto outOfLoop; //-----------------------> 00401 } 00402 } 00403 } 00404 } 00405 } 00406 00407 // 3 previous loops continued to the end 00408 inElastic = false; // quasi-elastic scattering 00409 return; 00410 } 00411 else 00412 { // target must be a neutron 00413 counter = -1; 00414 for( npos=0; npos<numSec/3; npos++ ) 00415 { 00416 for( nneg=npos; nneg<=(npos+1); nneg++ ) 00417 { 00418 for( nzero=0; nzero<numSec/3; nzero++ ) 00419 { 00420 if( ++counter < numMul ) 00421 { 00422 nt = npos+nneg+nzero; 00423 if ( (nt>=1) && (nt<=numSec) ) { 00424 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00425 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 00426 if (std::fabs(dum) < 1.0) { 00427 if( test >= 1.0e-10 )excs += dum*test; 00428 } else { 00429 excs += dum*test; 00430 } 00431 if (ran < excs) goto outOfLoop; // --------------------------> 00432 } 00433 } 00434 } 00435 } 00436 } 00437 // 3 previous loops continued to the end 00438 inElastic = false; // quasi-elastic scattering. 00439 return; 00440 } 00441 00442 outOfLoop: // <------------------------------------------------------------------------ 00443 00444 // in the following we do not consider strangeness transfer at high 00445 // multiplicity events. YK combinations are added in 00446 // StrangeParticlePairProduction 00447 ran = G4UniformRand(); 00448 if (targetCode == protonCode) { 00449 if( npos == nneg) 00450 { 00451 } 00452 else 00453 { 00454 pv[1] = Neutron; 00455 } 00456 } else { 00457 if (npos == nneg) 00458 { 00459 } 00460 else 00461 { 00462 pv[1] = Proton; 00463 } 00464 } 00465 00466 nt = npos + nneg + nzero; 00467 while (nt > 0) { 00468 G4double rnd = G4UniformRand(); 00469 if (rnd < (G4double)npos/nt) { 00470 if (npos > 0) { 00471 pv[vecLen++] = PionPlus; 00472 npos--; 00473 } 00474 } else if (rnd < (G4double)(npos+nneg)/nt) { 00475 if (nneg > 0) { 00476 pv[vecLen++] = PionMinus; 00477 nneg--; 00478 } 00479 } else { 00480 if (nzero > 0) { 00481 pv[vecLen++] = PionZero; 00482 nzero--; 00483 } 00484 } 00485 nt = npos + nneg + nzero; 00486 } 00487 00488 if (verboseLevel > 1) { 00489 G4cout << "Particles produced: " ; 00490 G4cout << pv[0].getName() << " " ; 00491 G4cout << pv[1].getName() << " " ; 00492 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ; 00493 G4cout << G4endl; 00494 } 00495 return; 00496 }
G4int G4HEOmegaMinusInelastic::GetNumberOfSecondaries | ( | ) | [inline] |
Definition at line 75 of file G4HEOmegaMinusInelastic.hh.
References vecLength.
00075 {return vecLength;}
void G4HEOmegaMinusInelastic::ModelDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 42 of file G4HEOmegaMinusInelastic.cc.
00043 { 00044 outFile << "G4HEOmegaMinusInelastic is one of the High Energy\n" 00045 << "Parameterized (HEP) models used to implement inelastic\n" 00046 << "Omega- 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 Omega- with initial energies\n" 00052 << "above 20 GeV.\n"; 00053 }
Definition at line 70 of file G4HEOmegaMinusInelastic.hh.
Referenced by ApplyYourself(), G4HEOmegaMinusInelastic(), and GetNumberOfSecondaries().