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 "G4HEAntiNeutronInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042
00043 G4HEAntiNeutronInelastic::G4HEAntiNeutronInelastic(const G4String& name)
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 G4HEAntiNeutronInelastic is being deprecated and will\n"
00052 << "disappear in Geant4 version 10.0" << G4endl;
00053 }
00054
00055
00056 void G4HEAntiNeutronInelastic::ModelDescription(std::ostream& outFile) const
00057 {
00058 outFile << "G4HEAntiNeutronInelastic is one of the High Energy\n"
00059 << "Parameterized (HEP) models used to implement inelastic\n"
00060 << "anti-neutron scattering from nuclei. It is a re-engineered\n"
00061 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00062 << "initial collision products into backward- and forward-going\n"
00063 << "clusters which are then decayed into final state hadrons.\n"
00064 << "The model does not conserve energy on an event-by-event\n"
00065 << "basis. It may be applied to anti-neutrons with initial energies\n"
00066 << "above 20 GeV.\n";
00067 }
00068
00069
00070 G4HadFinalState*
00071 G4HEAntiNeutronInelastic::ApplyYourself(const G4HadProjectile &aTrack,
00072 G4Nucleus &targetNucleus)
00073 {
00074 G4HEVector* pv = new G4HEVector[MAXPART];
00075 const G4HadProjectile* aParticle = &aTrack;
00076 const G4double atomicWeight = targetNucleus.GetA_asInt();
00077 const G4double atomicNumber = targetNucleus.GetZ_asInt();
00078 G4HEVector incidentParticle(aParticle);
00079
00080 G4int incidentCode = incidentParticle.getCode();
00081 G4double incidentMass = incidentParticle.getMass();
00082 G4double incidentTotalEnergy = incidentParticle.getEnergy();
00083
00084
00085
00086
00087 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00088
00089 if (incidentKineticEnergy < 1.)
00090 G4cout << "GHEAntiNeutronInelastic: incident energy < 1 GeV" << G4endl;;
00091
00092 if (verboseLevel > 1) {
00093 G4cout << "G4HEAntiNeutronInelastic::ApplyYourself" << G4endl;
00094 G4cout << "incident particle " << incidentParticle.getName()
00095 << "mass " << incidentMass
00096 << "kinetic energy " << incidentKineticEnergy
00097 << G4endl;
00098 G4cout << "target material with (A,Z) = ("
00099 << atomicWeight << "," << atomicNumber << ")" << G4endl;
00100 }
00101
00102 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
00103 atomicWeight, atomicNumber);
00104 if (verboseLevel > 1)
00105 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00106
00107 incidentKineticEnergy -= inelasticity;
00108
00109 G4double excitationEnergyGNP = 0.;
00110 G4double excitationEnergyDTA = 0.;
00111
00112 G4double excitation = NuclearExcitation(incidentKineticEnergy,
00113 atomicWeight, atomicNumber,
00114 excitationEnergyGNP,
00115 excitationEnergyDTA);
00116 if (verboseLevel > 1)
00117 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
00118 << excitationEnergyDTA << G4endl;
00119
00120 incidentKineticEnergy -= excitation;
00121 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00122
00123
00124
00125
00126 G4HEVector targetParticle;
00127 if (G4UniformRand() < atomicNumber/atomicWeight) {
00128 targetParticle.setDefinition("Proton");
00129 } else {
00130 targetParticle.setDefinition("Neutron");
00131 }
00132
00133 G4double targetMass = targetParticle.getMass();
00134 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00135 + targetMass*targetMass
00136 + 2.0*targetMass*incidentTotalEnergy);
00137 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00138
00139 G4bool inElastic = true;
00140 vecLength = 0;
00141
00142 if (verboseLevel > 1)
00143 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00144 << incidentCode << G4endl;
00145
00146 G4bool successful = false;
00147
00148 FirstIntInCasAntiNeutron(inElastic, availableEnergy, pv, vecLength,
00149 incidentParticle, targetParticle, atomicWeight);
00150
00151 if (verboseLevel > 1)
00152 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00153
00154 if ((vecLength > 0) && (availableEnergy > 1.))
00155 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00156 pv, vecLength,
00157 incidentParticle, targetParticle);
00158 HighEnergyCascading(successful, pv, vecLength,
00159 excitationEnergyGNP, excitationEnergyDTA,
00160 incidentParticle, targetParticle,
00161 atomicWeight, atomicNumber);
00162 if (!successful)
00163 HighEnergyClusterProduction(successful, pv, vecLength,
00164 excitationEnergyGNP, excitationEnergyDTA,
00165 incidentParticle, targetParticle,
00166 atomicWeight, atomicNumber);
00167 if (!successful)
00168 MediumEnergyCascading(successful, pv, vecLength,
00169 excitationEnergyGNP, excitationEnergyDTA,
00170 incidentParticle, targetParticle,
00171 atomicWeight, atomicNumber);
00172
00173 if (!successful)
00174 MediumEnergyClusterProduction(successful, pv, vecLength,
00175 excitationEnergyGNP, excitationEnergyDTA,
00176 incidentParticle, targetParticle,
00177 atomicWeight, atomicNumber);
00178 if (!successful)
00179 QuasiElasticScattering(successful, pv, vecLength,
00180 excitationEnergyGNP, excitationEnergyDTA,
00181 incidentParticle, targetParticle,
00182 atomicWeight, atomicNumber);
00183 if (!successful)
00184 ElasticScattering(successful, pv, vecLength,
00185 incidentParticle,
00186 atomicWeight, atomicNumber);
00187
00188 if (!successful)
00189 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00190 << G4endl;
00191
00192 FillParticleChange(pv, vecLength);
00193 delete [] pv;
00194 theParticleChange.SetStatusChange(stopAndKill);
00195 return &theParticleChange;
00196 }
00197
00198
00199 void
00200 G4HEAntiNeutronInelastic::FirstIntInCasAntiNeutron(G4bool& inElastic,
00201 const G4double availableEnergy,
00202 G4HEVector pv[],
00203 G4int& vecLen,
00204 const G4HEVector& incidentParticle,
00205 const G4HEVector& targetParticle,
00206 const G4double atomicWeight)
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 {
00217 static const G4double expxu = 82.;
00218 static const G4double expxl = -expxu;
00219
00220 static const G4double protb = 0.7;
00221 static const G4double neutb = 0.7;
00222 static const G4double c = 1.25;
00223
00224 static const G4int numMul = 1200;
00225 static const G4int numMulAn = 400;
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 protmulAn[numMulAn],protnormAn[numSec];
00237 static G4double neutmul[numMul], neutnorm[numSec];
00238 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
00239
00240
00241
00242
00243 G4int i, counter, nt, npos, nneg, nzero;
00244
00245 if( first )
00246 {
00247 first = false;
00248 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
00249 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
00250 counter = -1;
00251 for( npos=0; npos<(numSec/3); npos++ )
00252 {
00253 for( nneg=std::max(0,npos-2); nneg<=npos; 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 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00263 protnorm[nt-1] += protmul[counter];
00264 }
00265 }
00266 }
00267 }
00268 }
00269 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00270 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00271 counter = -1;
00272 for( npos=0; npos<numSec/3; npos++ )
00273 {
00274 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
00275 {
00276 for( nzero=0; nzero<numSec/3; nzero++ )
00277 {
00278 if( ++counter < numMul )
00279 {
00280 nt = npos+nneg+nzero;
00281 if( (nt>0) && (nt<=numSec) )
00282 {
00283 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00284 neutnorm[nt-1] += neutmul[counter];
00285 }
00286 }
00287 }
00288 }
00289 }
00290 for( i=0; i<numSec; i++ )
00291 {
00292 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00293 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00294 }
00295
00296 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
00297 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
00298 counter = -1;
00299 for( npos=1; npos<(numSec/3); npos++ )
00300 {
00301 nneg = std::max(0,npos-1);
00302 for( nzero=0; nzero<numSec/3; nzero++ )
00303 {
00304 if( ++counter < numMulAn )
00305 {
00306 nt = npos+nneg+nzero;
00307 if( (nt>1) && (nt<=numSec) )
00308 {
00309 protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00310 protnormAn[nt-1] += protmulAn[counter];
00311 }
00312 }
00313 }
00314 }
00315 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
00316 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
00317 counter = -1;
00318 for( npos=0; npos<numSec/3; npos++ )
00319 {
00320 nneg = npos;
00321 for( nzero=0; nzero<numSec/3; nzero++ )
00322 {
00323 if( ++counter < numMulAn )
00324 {
00325 nt = npos+nneg+nzero;
00326 if( (nt>1) && (nt<=numSec) )
00327 {
00328 neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00329 neutnormAn[nt-1] += neutmulAn[counter];
00330 }
00331 }
00332 }
00333 }
00334 for( i=0; i<numSec; i++ )
00335 {
00336 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
00337 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
00338 }
00339 }
00340
00341
00342
00343
00344 pv[0] = incidentParticle;
00345 pv[1] = targetParticle;
00346 vecLen = 2;
00347
00348 if( !inElastic )
00349 {
00350 if( targetCode == neutronCode )
00351 {
00352 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00353
00354 G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5));
00355 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
00356 {
00357 pv[0] = AntiProton;
00358 pv[1] = Proton;
00359 }
00360 }
00361 return;
00362 }
00363 else if (availableEnergy <= PionPlus.getMass())
00364 return;
00365
00366
00367
00368 npos = 0, nneg = 0, nzero = 0;
00369 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
00370 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
00371 0.39, 0.36, 0.33, 0.10, 0.01};
00372 G4int iplab = G4int( incidentTotalMomentum*10.);
00373 if ( iplab > 9) iplab = 10 + G4int( (incidentTotalMomentum -1.)*5. );
00374 if ( iplab > 14) iplab = 15 + G4int( incidentTotalMomentum -2. );
00375 if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.);
00376 iplab = std::min(24, iplab);
00377
00378 if ( G4UniformRand() > anhl[iplab] )
00379 {
00380
00381 G4double eab = availableEnergy;
00382 G4int ieab = G4int( eab*5.0 );
00383
00384 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
00385 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
00386 {
00387
00388
00389 G4double w0, wp, wm, wt, ran;
00390 if( targetCode == protonCode )
00391 {
00392 w0 = - sqr(1.+protb)/(2.*c*c);
00393 w0 = wp = std::exp(w0);
00394 if( G4UniformRand() < w0/(w0+wp) )
00395 { npos = 0; nneg = 0; nzero = 1; }
00396 else
00397 { npos = 1; nneg = 0; nzero = 0; }
00398 }
00399 else
00400 {
00401 w0 = -sqr(1.+neutb)/(2.*c*c);
00402 w0 = wp = std::exp(w0);
00403 wm = -sqr(-1.+neutb)/(2.*c*c);
00404 wm = std::exp(wm);
00405 wt = w0+wp+wm;
00406 wp = w0+wp;
00407 ran = G4UniformRand();
00408 if( ran < w0/wt)
00409 { npos = 0; nneg = 0; nzero = 1; }
00410 else if( ran < wp/wt)
00411 { npos = 1; nneg = 0; nzero = 0; }
00412 else
00413 { npos = 0; nneg = 1; nzero = 0; }
00414 }
00415 }
00416 else
00417 {
00418
00419
00420 G4double aleab = std::log(availableEnergy);
00421 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00422 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00423
00424
00425
00426 G4double test, dum, anpn = 0.0;
00427
00428 for (nt=1; nt<=numSec; nt++) {
00429 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00430 dum = pi*nt/(2.0*n*n);
00431
00432 if (std::fabs(dum) < 1.0) {
00433 if( test >= 1.0e-10 )anpn += dum*test;
00434 } else {
00435 anpn += dum*test;
00436 }
00437 }
00438
00439 G4double ran = G4UniformRand();
00440 G4double excs = 0.0;
00441 if (targetCode == protonCode) {
00442 counter = -1;
00443 for (npos=0; npos<numSec/3; npos++) {
00444 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
00445 for (nzero=0; nzero<numSec/3; nzero++) {
00446 if (++counter < numMul) {
00447 nt = npos+nneg+nzero;
00448 if ( (nt>0) && (nt<=numSec) ) {
00449 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00450 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00451
00452 if (std::fabs(dum) < 1.0) {
00453 if( test >= 1.0e-10 )excs += dum*test;
00454 } else {
00455 excs += dum*test;
00456 }
00457
00458 if (ran < excs) goto outOfLoop;
00459 }
00460 }
00461 }
00462 }
00463 }
00464
00465
00466 inElastic = false;
00467 return;
00468
00469 } else {
00470 counter = -1;
00471 for (npos=0; npos<numSec/3; npos++) {
00472 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00473 for (nzero=0; nzero<numSec/3; nzero++) {
00474 if (++counter < numMul) {
00475 nt = npos+nneg+nzero;
00476 if ((nt>0) && (nt<=numSec) ) {
00477 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00478 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00479 if (std::fabs(dum) < 1.0) {
00480 if( test >= 1.0e-10 )excs += dum*test;
00481 } else {
00482 excs += dum*test;
00483 }
00484
00485 if (ran < excs) goto outOfLoop;
00486 }
00487 }
00488 }
00489 }
00490 }
00491
00492 inElastic = false;
00493 return;
00494 }
00495 }
00496 outOfLoop:
00497
00498 if( targetCode == protonCode)
00499 {
00500 if( npos == nneg)
00501 {
00502 }
00503 else if (npos == (nneg+1))
00504 {
00505 if( G4UniformRand() < 0.5)
00506 {
00507 pv[1] = Neutron;
00508 }
00509 else
00510 {
00511 pv[0] = AntiProton;
00512 }
00513 }
00514 else
00515 {
00516 pv[0] = AntiProton;
00517 pv[1] = Neutron;
00518 }
00519 }
00520 else
00521 {
00522 if( npos == nneg)
00523 {
00524 if( G4UniformRand() < 0.25)
00525 {
00526 pv[0] = AntiProton;
00527 pv[1] = Proton;
00528 }
00529 else
00530 {
00531 }
00532 }
00533 else if ( npos == (nneg-1))
00534 {
00535 pv[1] = Proton;
00536 }
00537 else
00538 {
00539 pv[0] = AntiProton;
00540 }
00541 }
00542
00543 }
00544 else
00545 {
00546 if ( availableEnergy > 2. * PionPlus.getMass() )
00547 {
00548
00549 G4double aleab = std::log(availableEnergy);
00550 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00551 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00552
00553
00554
00555 G4double test, dum, anpn = 0.0;
00556
00557 for (nt=2; nt<=numSec; nt++) {
00558 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00559 dum = pi*nt/(2.0*n*n);
00560 if (std::fabs(dum) < 1.0) {
00561 if( test >= 1.0e-10 )anpn += dum*test;
00562 } else {
00563 anpn += dum*test;
00564 }
00565 }
00566
00567 G4double ran = G4UniformRand();
00568 G4double excs = 0.0;
00569 if (targetCode == protonCode) {
00570 counter = -1;
00571 for (npos=1; npos<numSec/3; npos++) {
00572 nneg = npos-1;
00573 for (nzero=0; nzero<numSec/3; nzero++) {
00574 if (++counter < numMulAn) {
00575 nt = npos+nneg+nzero;
00576 if ( (nt>1) && (nt<=numSec) ) {
00577 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00578 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
00579
00580 if (std::fabs(dum) < 1.0) {
00581 if( test >= 1.0e-10 )excs += dum*test;
00582 } else {
00583 excs += dum*test;
00584 }
00585
00586 if (ran < excs) goto outOfLoopAn;
00587 }
00588 }
00589 }
00590 }
00591
00592 inElastic = false;
00593 return;
00594
00595 } else {
00596 counter = -1;
00597 for (npos=0; npos<numSec/3; npos++) {
00598 nneg = npos;
00599 for (nzero=0; nzero<numSec/3; nzero++) {
00600 if (++counter < numMulAn) {
00601 nt = npos+nneg+nzero;
00602 if ( (nt>1) && (nt<=numSec) ) {
00603 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00604 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
00605
00606 if (std::fabs(dum) < 1.0) {
00607 if( test >= 1.0e-10 )excs += dum*test;
00608 } else {
00609 excs += dum*test;
00610 }
00611
00612 if (ran < excs) goto outOfLoopAn;
00613 }
00614 }
00615 }
00616 }
00617 inElastic = false;
00618 return;
00619 }
00620 outOfLoopAn:
00621 vecLen = 0;
00622 }
00623 }
00624
00625 nt = npos + nneg + nzero;
00626 while ( nt > 0)
00627 {
00628 G4double ran = G4UniformRand();
00629 if ( ran < (G4double)npos/nt)
00630 {
00631 if( npos > 0 )
00632 { pv[vecLen++] = PionPlus;
00633 npos--;
00634 }
00635 }
00636 else if ( ran < (G4double)(npos+nneg)/nt)
00637 {
00638 if( nneg > 0 )
00639 {
00640 pv[vecLen++] = PionMinus;
00641 nneg--;
00642 }
00643 }
00644 else
00645 {
00646 if( nzero > 0 )
00647 {
00648 pv[vecLen++] = PionZero;
00649 nzero--;
00650 }
00651 }
00652 nt = npos + nneg + nzero;
00653 }
00654 if (verboseLevel > 1)
00655 {
00656 G4cout << "Particles produced: " ;
00657 G4cout << pv[0].getName() << " " ;
00658 G4cout << pv[1].getName() << " " ;
00659 for (i=2; i < vecLen; i++)
00660 {
00661 G4cout << pv[i].getName() << " " ;
00662 }
00663 G4cout << G4endl;
00664 }
00665 return;
00666 }
00667
00668
00669
00670
00671
00672
00673
00674
00675