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