00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "G4HESigmaMinusInelastic.hh"
00037 #include "globals.hh"
00038 #include "G4ios.hh"
00039 #include "G4PhysicalConstants.hh"
00040 #include "G4SystemOfUnits.hh"
00041
00042 G4HESigmaMinusInelastic::G4HESigmaMinusInelastic(const G4String& name)
00043 : G4HEInelastic(name)
00044 {
00045 vecLength = 0;
00046 theMinEnergy = 20*GeV;
00047 theMaxEnergy = 10*TeV;
00048 MAXPART = 2048;
00049 verboseLevel = 0;
00050 G4cout << "WARNING: model G4HESigmaMinusInelastic is being deprecated and will\n"
00051 << "disappear in Geant4 version 10.0" << G4endl;
00052 }
00053
00054
00055 void G4HESigmaMinusInelastic::ModelDescription(std::ostream& outFile) const
00056 {
00057 outFile << "G4HESigmaMinusInelastic is one of the High Energy Parameterized\n"
00058 << "(HEP) models used to implement inelastic Sigma- scattering\n"
00059 << "from nuclei. It is a re-engineered version of the GHEISHA\n"
00060 << "code of H. Fesefeldt. It divides the initial collision\n"
00061 << "products into backward- and forward-going clusters which are\n"
00062 << "then decayed into final state hadrons. The model does not\n"
00063 << "conserve energy on an event-by-event basis. It may be\n"
00064 << "applied to sigmas with initial energies above 20 GeV.\n";
00065 }
00066
00067
00068 G4HadFinalState*
00069 G4HESigmaMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00070 G4Nucleus& targetNucleus)
00071 {
00072 G4HEVector* pv = new G4HEVector[MAXPART];
00073 const G4HadProjectile* aParticle = &aTrack;
00074 const G4double A = targetNucleus.GetA_asInt();
00075 const G4double Z = targetNucleus.GetZ_asInt();
00076 G4HEVector incidentParticle(aParticle);
00077
00078 G4double atomicNumber = Z;
00079 G4double atomicWeight = A;
00080
00081 G4int incidentCode = incidentParticle.getCode();
00082 G4double incidentMass = incidentParticle.getMass();
00083 G4double incidentTotalEnergy = incidentParticle.getEnergy();
00084
00085
00086
00087
00088 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00089
00090 if (incidentKineticEnergy < 1.)
00091 G4cout << "GHESigmaMinusInelastic: incident energy < 1 GeV" << G4endl;
00092
00093 if (verboseLevel > 1) {
00094 G4cout << "G4HESigmaMinusInelastic::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 incidentKineticEnergy -= excitation;
00122 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00123
00124
00125
00126
00127 G4HEVector targetParticle;
00128 if (G4UniformRand() < atomicNumber/atomicWeight) {
00129 targetParticle.setDefinition("Proton");
00130 } else {
00131 targetParticle.setDefinition("Neutron");
00132 }
00133
00134 G4double targetMass = targetParticle.getMass();
00135 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00136 + targetMass*targetMass
00137 + 2.0*targetMass*incidentTotalEnergy);
00138 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00139
00140 G4bool inElastic = true;
00141 vecLength = 0;
00142
00143 if (verboseLevel > 1)
00144 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00145 << incidentCode << G4endl;
00146
00147 G4bool successful = false;
00148
00149 FirstIntInCasSigmaMinus(inElastic, availableEnergy, pv, vecLength,
00150 incidentParticle, targetParticle, atomicWeight);
00151
00152 if (verboseLevel > 1)
00153 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00154
00155 if ((vecLength > 0) && (availableEnergy > 1.))
00156 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00157 pv, vecLength,
00158 incidentParticle, targetParticle);
00159
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
00175 if (!successful)
00176 MediumEnergyClusterProduction(successful, pv, vecLength,
00177 excitationEnergyGNP, excitationEnergyDTA,
00178 incidentParticle, targetParticle,
00179 atomicWeight, atomicNumber);
00180 if (!successful)
00181 QuasiElasticScattering(successful, pv, vecLength,
00182 excitationEnergyGNP, excitationEnergyDTA,
00183 incidentParticle, targetParticle,
00184 atomicWeight, atomicNumber);
00185 if (!successful)
00186 ElasticScattering(successful, pv, vecLength,
00187 incidentParticle,
00188 atomicWeight, atomicNumber);
00189
00190 if (!successful)
00191 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00192 << G4endl;
00193
00194 FillParticleChange(pv, vecLength);
00195 delete [] pv;
00196 theParticleChange.SetStatusChange(stopAndKill);
00197 return &theParticleChange;
00198 }
00199
00200
00201 void
00202 G4HESigmaMinusInelastic::FirstIntInCasSigmaMinus(G4bool& inElastic,
00203 const G4double availableEnergy,
00204 G4HEVector pv[],
00205 G4int& vecLen,
00206 const G4HEVector& incidentParticle,
00207 const G4HEVector& targetParticle,
00208 const G4double atomicWeight)
00209
00210
00211
00212
00213
00214
00215
00216
00217 {
00218 static const G4double expxu = 82.;
00219 static const G4double expxl = -expxu;
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 neutronCode = Neutron.getCode();
00229 G4int protonCode = Proton.getCode();
00230
00231 G4int targetCode = targetParticle.getCode();
00232 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00233
00234 static G4bool first = true;
00235 static G4double protmul[numMul], protnorm[numSec];
00236 static G4double neutmul[numMul], neutnorm[numSec];
00237
00238
00239
00240
00241 G4int i, counter, nt, npos, nneg, nzero;
00242
00243 if (first) {
00244 first = false;
00245 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
00246 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00247 counter = -1;
00248 for (npos = 0; npos < (numSec/3); npos++) {
00249 for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
00250 for (nzero = 0; nzero < numSec/3; nzero++) {
00251 if (++counter < numMul) {
00252 nt = npos+nneg+nzero;
00253 if ((nt > 0) && (nt<=numSec) ) {
00254 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00255 protnorm[nt-1] += protmul[counter];
00256 }
00257 }
00258 }
00259 }
00260 }
00261
00262 for (i = 0; i < numMul; i++) neutmul[i] = 0.0;
00263 for (i = 0; i < numSec; i++) neutnorm[i] = 0.0;
00264 counter = -1;
00265 for (npos = 0; npos < numSec/3; npos++) {
00266 for (nneg = npos; nneg <= (npos+2); nneg++) {
00267 for (nzero = 0; nzero < numSec/3; nzero++) {
00268 if (++counter < numMul) {
00269 nt = npos+nneg+nzero;
00270 if ((nt>0) && (nt<=numSec) ) {
00271 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00272 neutnorm[nt-1] += neutmul[counter];
00273 }
00274 }
00275 }
00276 }
00277 }
00278 for (i = 0; i < numSec; i++) {
00279 if (protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00280 if (neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00281 }
00282 }
00283
00284 pv[0] = incidentParticle;
00285 pv[1] = targetParticle;
00286 vecLen = 2;
00287
00288 if (!inElastic) {
00289 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00290 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
00291 if (G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) {
00292 G4double ran = G4UniformRand();
00293 if (targetCode == neutronCode) {
00294 pv[0] = Neutron;
00295 pv[1] = SigmaMinus;
00296 } else {
00297 if (ran < 0.2) {
00298 pv[0] = SigmaZero;
00299 pv[1] = Neutron;
00300 } else if (ran < 0.4) {
00301 pv[0] = Lambda;
00302 pv[1] = Neutron;
00303 } else if (ran < 0.6) {
00304 pv[0] = Proton;
00305 pv[1] = SigmaMinus;
00306 } else if (ran < 0.8) {
00307 pv[0] = Neutron;
00308 pv[1] = SigmaZero;
00309 } else {
00310 pv[0] = Neutron;
00311 pv[1] = Lambda;
00312 }
00313 }
00314 }
00315 return;
00316
00317 } else if (availableEnergy <= PionPlus.getMass()) return;
00318
00319
00320 npos = 0, nneg = 0, nzero = 0;
00321
00322
00323 G4double aleab = std::log(availableEnergy);
00324 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00325 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00326
00327
00328
00329 G4double test, dum, anpn = 0.0;
00330
00331 for (nt = 1; nt <= numSec; nt++) {
00332 test = std::exp(std::min(expxu, std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00333 dum = pi*nt/(2.0*n*n);
00334 if (std::fabs(dum) < 1.0) {
00335 if (test >= 1.0e-10) anpn += dum*test;
00336 } else {
00337 anpn += dum*test;
00338 }
00339 }
00340
00341 G4double ran = G4UniformRand();
00342 G4double excs = 0.0;
00343 if (targetCode == protonCode) {
00344 counter = -1;
00345 for (npos=0; npos<numSec/3; npos++) {
00346 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00347 for (nzero=0; nzero<numSec/3; nzero++) {
00348 if (++counter < numMul) {
00349 nt = npos+nneg+nzero;
00350 if ( (nt>0) && (nt<=numSec) ) {
00351 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00352 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00353 if (std::fabs(dum) < 1.0) {
00354 if( test >= 1.0e-10 )excs += dum*test;
00355 } else {
00356 excs += dum*test;
00357 }
00358 if (ran < excs) goto outOfLoop;
00359 }
00360 }
00361 }
00362 }
00363 }
00364
00365
00366 inElastic = false;
00367 return;
00368
00369 } else {
00370 counter = -1;
00371 for (npos=0; npos<numSec/3; npos++) {
00372 for (nneg=npos; nneg<=(npos+2); nneg++) {
00373 for (nzero=0; nzero<numSec/3; nzero++) {
00374 if (++counter < numMul) {
00375 nt = npos+nneg+nzero;
00376 if ( (nt>=1) && (nt<=numSec) ) {
00377 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00378 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00379 if (std::fabs(dum) < 1.0) {
00380 if( test >= 1.0e-10 )excs += dum*test;
00381 } else {
00382 excs += dum*test;
00383 }
00384 if (ran < excs) goto outOfLoop;
00385 }
00386 }
00387 }
00388 }
00389 }
00390
00391
00392 inElastic = false;
00393 return;
00394 }
00395
00396 outOfLoop:
00397
00398 ran = G4UniformRand();
00399 if (targetCode == neutronCode) {
00400 if( npos == nneg)
00401 {
00402 }
00403 else if (npos == (nneg-1))
00404 {
00405 if( ran < 0.25)
00406 {
00407 pv[0] = SigmaZero;
00408 }
00409 else if (ran < 0.5)
00410 {
00411 pv[0] = Lambda;
00412 }
00413 else
00414 {
00415 pv[1] = Proton;
00416 }
00417 }
00418 else
00419 {
00420 if(ran < 0.5)
00421 {
00422 pv[0] = SigmaZero;
00423 pv[1] = Proton;
00424 }
00425 else
00426 {
00427 pv[0] = Lambda;
00428 pv[1] = Proton;
00429 }
00430 }
00431 } else {
00432 if (npos == nneg)
00433 {
00434 if (ran < 0.5)
00435 {
00436 }
00437 else if (ran < 0.75)
00438 {
00439 pv[0] = SigmaZero;
00440 pv[1] = Neutron;
00441 }
00442 else
00443 {
00444 pv[0] = Lambda;
00445 pv[1] = Neutron;
00446 }
00447 }
00448 else if (npos == (nneg+1))
00449 {
00450 pv[1] = Neutron;
00451 }
00452 else
00453 {
00454 if (ran < 0.5)
00455 {
00456 pv[0] = SigmaZero;
00457 }
00458 else
00459 {
00460 pv[0] = Lambda;
00461 }
00462 }
00463 }
00464
00465 nt = npos + nneg + nzero;
00466 while (nt > 0) {
00467 G4double rnd = G4UniformRand();
00468 if (rnd < (G4double)npos/nt) {
00469 if (npos > 0) {
00470 pv[vecLen++] = PionPlus;
00471 npos--;
00472 }
00473 } else if (rnd < (G4double)(npos+nneg)/nt) {
00474 if (nneg > 0) {
00475 pv[vecLen++] = PionMinus;
00476 nneg--;
00477 }
00478 } else {
00479 if (nzero > 0) {
00480 pv[vecLen++] = PionZero;
00481 nzero--;
00482 }
00483 }
00484 nt = npos + nneg + nzero;
00485 }
00486
00487 if (verboseLevel > 1) {
00488 G4cout << "Particles produced: " ;
00489 G4cout << pv[0].getCode() << " " ;
00490 G4cout << pv[1].getCode() << " " ;
00491 for (i = 2; i < vecLen; i++) G4cout << pv[i].getCode() << " " ;
00492 G4cout << G4endl;
00493 }
00494
00495 return;
00496 }
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506