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
00037 #include "G4HENeutronInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041
00042 void G4HENeutronInelastic::ModelDescription(std::ostream& outFile) const
00043 {
00044 outFile << "G4HENeutronInelastic is one of the High Energy\n"
00045 << "Parameterized (HEP) models used to implement inelastic\n"
00046 << "neutron 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 neutrons with initial energies\n"
00052 << "above 20 GeV.\n";
00053 }
00054
00055
00056 G4HadFinalState*
00057 G4HENeutronInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00058 G4Nucleus& targetNucleus)
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
00074
00075
00076 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00077
00078 if (incidentKineticEnergy < 1.)
00079 G4cout << "GHENeutronInelastic: incident energy < 1 GeV" << G4endl;
00080
00081 if (verboseLevel > 1) {
00082 G4cout << "G4HENeutronInelastic::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
00112
00113
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
00129
00130
00131
00132
00133 G4bool inElastic = true;
00134
00135
00136 vecLength = 0;
00137
00138 if (verboseLevel > 1)
00139 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00140 << incidentCode << G4endl;
00141
00142 G4bool successful = false;
00143
00144 FirstIntInCasNeutron(inElastic, availableEnergy, pv, vecLength,
00145 incidentParticle, targetParticle, atomicWeight);
00146
00147 if (verboseLevel > 1)
00148 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00149
00150 if ((vecLength > 0) && (availableEnergy > 1.))
00151 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00152 pv, vecLength,
00153 incidentParticle, targetParticle);
00154
00155 HighEnergyCascading(successful, pv, vecLength,
00156 excitationEnergyGNP, excitationEnergyDTA,
00157 incidentParticle, targetParticle,
00158 atomicWeight, atomicNumber);
00159 if (!successful)
00160 HighEnergyClusterProduction(successful, pv, vecLength,
00161 excitationEnergyGNP, excitationEnergyDTA,
00162 incidentParticle, targetParticle,
00163 atomicWeight, atomicNumber);
00164 if (!successful)
00165 MediumEnergyCascading(successful, pv, vecLength,
00166 excitationEnergyGNP, excitationEnergyDTA,
00167 incidentParticle, targetParticle,
00168 atomicWeight, atomicNumber);
00169
00170 if (!successful)
00171 MediumEnergyClusterProduction(successful, pv, vecLength,
00172 excitationEnergyGNP, excitationEnergyDTA,
00173 incidentParticle, targetParticle,
00174 atomicWeight, atomicNumber);
00175 if (!successful)
00176 QuasiElasticScattering(successful, pv, vecLength,
00177 excitationEnergyGNP, excitationEnergyDTA,
00178 incidentParticle, targetParticle,
00179 atomicWeight, atomicNumber);
00180 if (!successful)
00181 ElasticScattering(successful, pv, vecLength,
00182 incidentParticle,
00183 atomicWeight, atomicNumber);
00184
00185 if (!successful)
00186 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00187 << G4endl;
00188
00189 FillParticleChange(pv, vecLength);
00190 delete [] pv;
00191 theParticleChange.SetStatusChange(stopAndKill);
00192 return &theParticleChange;
00193 }
00194
00195
00196 void
00197 G4HENeutronInelastic::FirstIntInCasNeutron(G4bool& inElastic,
00198 const G4double availableEnergy,
00199 G4HEVector pv[],
00200 G4int& vecLen,
00201 const G4HEVector& incidentParticle,
00202 const G4HEVector& targetParticle,
00203 const G4double atomicWeight)
00204
00205
00206
00207
00208
00209
00210
00211
00212 {
00213 static const G4double expxu = 82.;
00214 static const G4double expxl = -expxu;
00215
00216 static const G4double protb = 0.35;
00217 static const G4double neutb = 0.35;
00218 static const G4double c = 1.25;
00219
00220 static const G4int numMul = 1200;
00221 static const G4int numSec = 60;
00222
00223 G4int neutronCode = Neutron.getCode();
00224 G4int protonCode = Proton.getCode();
00225
00226 G4int targetCode = targetParticle.getCode();
00227 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00228
00229 static G4bool first = true;
00230 static G4double protmul[numMul], protnorm[numSec];
00231 static G4double neutmul[numMul], neutnorm[numSec];
00232
00233
00234
00235
00236 G4int i, counter, nt, npos, nneg, nzero;
00237
00238 if (first) {
00239
00240 first = false;
00241 for (i = 0; i < numMul; i++) protmul[i] = 0.0;
00242 for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
00243 counter = -1;
00244 for (npos = 0; npos < (numSec/3); npos++) {
00245 for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
00246 for( nzero=0; nzero<numSec/3; nzero++ )
00247 {
00248 if( ++counter < numMul )
00249 {
00250 nt = npos+nneg+nzero;
00251 if( (nt>0) && (nt<=numSec) )
00252 {
00253 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c)
00254 /(Factorial(1-npos+nneg)*Factorial(1+npos-nneg)) ;
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 {
00267 for( nneg=npos; nneg<=(npos+2); nneg++ )
00268 {
00269 for( nzero=0; nzero<numSec/3; nzero++ )
00270 {
00271 if( ++counter < numMul )
00272 {
00273 nt = npos+nneg+nzero;
00274 if( (nt>0) && (nt<=numSec) )
00275 {
00276 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c)
00277 /(Factorial(-npos+nneg)*Factorial(2+npos-nneg));
00278 neutnorm[nt-1] += neutmul[counter];
00279 }
00280 }
00281 }
00282 }
00283 }
00284 for( i=0; i<numSec; i++ )
00285 {
00286 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00287 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00288 }
00289 }
00290
00291
00292
00293
00294 pv[0] = incidentParticle;
00295 pv[1] = targetParticle;
00296 vecLen = 2;
00297
00298 if( !inElastic )
00299 {
00300 if( targetCode == protonCode )
00301 {
00302 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00303 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5) );
00304 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
00305 {
00306 pv[0] = Proton;
00307 pv[1] = Neutron;
00308 }
00309 }
00310 return;
00311 }
00312 else if (availableEnergy <= PionPlus.getMass())
00313 return;
00314
00315
00316
00317 npos = 0, nneg = 0, nzero = 0;
00318 G4double eab = availableEnergy;
00319 G4int ieab = G4int( eab*5.0 );
00320
00321 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
00322 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
00323 {
00324
00325
00326 G4double w0, wp, wm, wt, ran;
00327 if( targetCode == neutronCode )
00328 {
00329 w0 = - sqr(1.+neutb)/(2.*c*c);
00330 wm = w0 = std::exp(w0);
00331 w0 = w0/2.;
00332 if( G4UniformRand() < w0/(w0+wm) )
00333 { npos = 0; nneg = 0; nzero = 1; }
00334 else
00335 { npos = 0; nneg = 1; nzero = 0; }
00336 }
00337 else
00338 {
00339 w0 = -sqr(1.+protb)/(2.*c*c);
00340 w0 = std::exp(w0);
00341 wp = w0/2.;
00342 wm = -sqr(-1.+protb)/(2.*c*c);
00343 wm = std::exp(wm)/2.;
00344 wt = w0+wp+wm;
00345 wp = w0+wp;
00346 ran = G4UniformRand();
00347 if( ran < w0/wt)
00348 { npos = 0; nneg = 0; nzero = 1; }
00349 else if( ran < wp/wt)
00350 { npos = 1; nneg = 0; nzero = 0; }
00351 else
00352 { npos = 0; nneg = 1; nzero = 0; }
00353 }
00354 }
00355 else
00356 {
00357
00358
00359 G4double aleab = std::log(availableEnergy);
00360 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00361 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00362
00363
00364
00365 G4double test, dum, anpn = 0.0;
00366
00367 for (nt=1; nt<=numSec; nt++) {
00368 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00369 dum = pi*nt/(2.0*n*n);
00370 if (std::fabs(dum) < 1.0) {
00371 if( test >= 1.0e-10 )anpn += dum*test;
00372 } else {
00373 anpn += dum*test;
00374 }
00375 }
00376
00377 G4double ran = G4UniformRand();
00378 G4double excs = 0.0;
00379 if( targetCode == protonCode )
00380 {
00381 counter = -1;
00382 for (npos=0; npos<numSec/3; npos++) {
00383 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00384 for (nzero=0; nzero<numSec/3; nzero++) {
00385 if (++counter < numMul) {
00386 nt = npos+nneg+nzero;
00387 if ( (nt>0) && (nt<=numSec) ) {
00388 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00389 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00390 if (std::fabs(dum) < 1.0) {
00391 if( test >= 1.0e-10 )excs += dum*test;
00392 } else {
00393 excs += dum*test;
00394 }
00395 if (ran < excs) goto outOfLoop;
00396 }
00397 }
00398 }
00399 }
00400 }
00401
00402
00403 inElastic = false;
00404 return;
00405 }
00406 else
00407 {
00408 counter = -1;
00409 for (npos=0; npos<numSec/3; npos++) {
00410 for (nneg=npos; nneg<=(npos+2); nneg++) {
00411 for (nzero=0; nzero<numSec/3; nzero++) {
00412 if (++counter < numMul) {
00413 nt = npos+nneg+nzero;
00414 if ( (nt>=1) && (nt<=numSec) ) {
00415 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00416 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00417 if (std::fabs(dum) < 1.0) {
00418 if( test >= 1.0e-10 )excs += dum*test;
00419 } else {
00420 excs += dum*test;
00421 }
00422 if (ran < excs) goto outOfLoop;
00423 }
00424 }
00425 }
00426 }
00427 }
00428
00429 inElastic = false;
00430 return;
00431 }
00432 }
00433 outOfLoop:
00434
00435 if( targetCode == neutronCode)
00436 {
00437 if( npos == nneg)
00438 {
00439 }
00440 else if (npos == (nneg-1))
00441 {
00442 if( G4UniformRand() < 0.5)
00443 {
00444 pv[1] = Proton;
00445 }
00446 else
00447 {
00448 pv[0] = Proton;
00449 }
00450 }
00451 else
00452 {
00453 pv[0] = Proton;
00454 pv[1] = Proton;
00455 }
00456 }
00457 else
00458 {
00459 if( npos == nneg)
00460 {
00461 if( G4UniformRand() < 0.25)
00462 {
00463 pv[0] = Proton;
00464 pv[1] = Neutron;
00465 }
00466 else
00467 {
00468 }
00469 }
00470 else if ( npos == (1+nneg))
00471 {
00472 pv[1] = Neutron;
00473 }
00474 else
00475 {
00476 pv[0] = Proton;
00477 }
00478 }
00479
00480
00481 nt = npos + nneg + nzero;
00482 while ( nt > 0)
00483 {
00484 G4double ran = G4UniformRand();
00485 if ( ran < (G4double)npos/nt)
00486 {
00487 if( npos > 0 )
00488 { pv[vecLen++] = PionPlus;
00489 npos--;
00490 }
00491 }
00492 else if ( ran < (G4double)(npos+nneg)/nt)
00493 {
00494 if( nneg > 0 )
00495 {
00496 pv[vecLen++] = PionMinus;
00497 nneg--;
00498 }
00499 }
00500 else
00501 {
00502 if( nzero > 0 )
00503 {
00504 pv[vecLen++] = PionZero;
00505 nzero--;
00506 }
00507 }
00508 nt = npos + nneg + nzero;
00509 }
00510 if (verboseLevel > 1)
00511 {
00512 G4cout << "Particles produced: " ;
00513 G4cout << pv[0].getName() << " " ;
00514 G4cout << pv[1].getName() << " " ;
00515 for (i=2; i < vecLen; i++)
00516 {
00517 G4cout << pv[i].getName() << " " ;
00518 }
00519 G4cout << G4endl;
00520 }
00521 return;
00522 }
00523
00524
00525
00526
00527
00528
00529
00530
00531