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 #include <iostream>
00036
00037 #include "G4LEAntiNeutronInelastic.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "Randomize.hh"
00041
00042 G4LEAntiNeutronInelastic::G4LEAntiNeutronInelastic(const G4String& name)
00043 :G4InelasticInteraction(name)
00044 {
00045 SetMinEnergy(0.0);
00046 SetMaxEnergy(25.*GeV);
00047 G4cout << "WARNING: model G4LEAntiNeutronInelastic is being deprecated and will\n"
00048 << "disappear in Geant4 version 10.0" << G4endl;
00049 }
00050
00051
00052 void G4LEAntiNeutronInelastic::ModelDescription(std::ostream& outFile) const
00053 {
00054 outFile << "G4LEAntiNeutronInelastic is one of the Low Energy\n"
00055 << "Parameterized (LEP) models used to implement inelastic\n"
00056 << "anti-neutron scattering from nuclei. It is a re-engineered\n"
00057 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00058 << "initial collision products into backward- and forward-going\n"
00059 << "clusters which are then decayed into final state hadrons.\n"
00060 << "The model does not conserve energy on an event-by-event\n"
00061 << "basis. It may be applied to anti-neutrons with initial\n"
00062 << "energies between 0 and 25 GeV.\n";
00063 }
00064
00065
00066 G4HadFinalState*
00067 G4LEAntiNeutronInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00068 G4Nucleus& targetNucleus)
00069 {
00070 const G4HadProjectile* originalIncident = &aTrack;
00071
00072
00073 G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00074
00075 if (verboseLevel > 1) {
00076 const G4Material *targetMaterial = aTrack.GetMaterial();
00077 G4cout << "G4LEAntiNeutronInelastic::ApplyYourself called" << G4endl;
00078 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00079 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00080 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00081 << G4endl;
00082 }
00083
00084
00085
00086 G4double ek = originalIncident->GetKineticEnergy()/MeV;
00087 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00088 G4ReactionProduct modifiedOriginal;
00089 modifiedOriginal = *originalIncident;
00090
00091 G4double tkin = targetNucleus.Cinema(ek);
00092 ek += tkin;
00093 modifiedOriginal.SetKineticEnergy(ek*MeV);
00094 G4double et = ek + amas;
00095 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00096 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00097 if (pp > 0.0) {
00098 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00099 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00100 }
00101
00102
00103 tkin = targetNucleus.EvaporationEffects(ek);
00104 ek -= tkin;
00105 modifiedOriginal.SetKineticEnergy(ek*MeV);
00106 et = ek + amas;
00107 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00108 pp = modifiedOriginal.GetMomentum().mag()/MeV;
00109 if (pp > 0.0) {
00110 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00111 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00112 }
00113
00114 G4ReactionProduct currentParticle = modifiedOriginal;
00115 G4ReactionProduct targetParticle;
00116 targetParticle = *originalTarget;
00117 currentParticle.SetSide(1);
00118 targetParticle.SetSide(-1);
00119 G4bool incidentHasChanged = false;
00120 G4bool targetHasChanged = false;
00121 G4bool quasiElastic = false;
00122 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;
00123 G4int vecLen = 0;
00124 vec.Initialize(0);
00125
00126 const G4double cutOff = 0.1*MeV;
00127 const G4double anni = std::min(1.3*currentParticle.GetTotalMomentum()/GeV, 0.4);
00128
00129 if ((currentParticle.GetKineticEnergy()/MeV > cutOff) ||
00130 (G4UniformRand() > anni) )
00131 Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00132 incidentHasChanged, targetHasChanged, quasiElastic);
00133 else
00134 quasiElastic = true;
00135
00136 CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00137 modifiedOriginal, targetNucleus, currentParticle,
00138 targetParticle, incidentHasChanged, targetHasChanged,
00139 quasiElastic);
00140
00141 SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00142
00143 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00144
00145 delete originalTarget;
00146 return &theParticleChange;
00147 }
00148
00149
00150 void G4LEAntiNeutronInelastic::Cascade(
00151 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00152 G4int& vecLen,
00153 const G4HadProjectile *originalIncident,
00154 G4ReactionProduct ¤tParticle,
00155 G4ReactionProduct &targetParticle,
00156 G4bool &incidentHasChanged,
00157 G4bool &targetHasChanged,
00158 G4bool &quasiElastic)
00159 {
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00172 const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00173 const G4double pOriginal = originalIncident->GetTotalMomentum()/MeV;
00174 const G4double targetMass = targetParticle.GetMass()/MeV;
00175 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00176 targetMass*targetMass +
00177 2.0*targetMass*etOriginal );
00178 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00179
00180 static G4bool first = true;
00181 const G4int numMul = 1200;
00182 const G4int numMulA = 400;
00183 const G4int numSec = 60;
00184 static G4double protmul[numMul], protnorm[numSec];
00185 static G4double neutmul[numMul], neutnorm[numSec];
00186 static G4double protmulA[numMulA], protnormA[numSec];
00187 static G4double neutmulA[numMulA], neutnormA[numSec];
00188
00189
00190 G4int counter;
00191 G4int nt=0;
00192 G4int npos = 0, nneg = 0, nzero = 0;
00193 G4double test;
00194 const G4double c = 1.25;
00195 const G4double b[] = { 0.70, 0.70 };
00196
00197 if (first) {
00198 first = false;
00199 G4int i;
00200 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00201 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00202 counter = -1;
00203 for (npos = 0; npos < (numSec/3); ++npos) {
00204 for (nneg = std::max(0,npos-2); nneg <= npos; ++nneg) {
00205 for (nzero = 0; nzero < numSec/3; ++nzero) {
00206 if (++counter < numMul) {
00207 nt = npos+nneg+nzero;
00208 if (nt > 0 && nt <= numSec) {
00209 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00210 protnorm[nt-1] += protmul[counter];
00211 }
00212 }
00213 }
00214 }
00215 }
00216 for (i = 0; i < numMul; ++i) neutmul[i] = 0.0;
00217 for (i = 0; i < numSec; ++i) neutnorm[i] = 0.0;
00218 counter = -1;
00219 for (npos = 0; npos < numSec/3; ++npos) {
00220 for (nneg = std::max(0,npos-1); nneg <= (npos+1); ++nneg) {
00221 for (nzero = 0; nzero < numSec/3; ++nzero) {
00222 if (++counter < numMul) {
00223 nt = npos+nneg+nzero;
00224 if ( (nt > 0) && (nt <= numSec) ) {
00225 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00226 neutnorm[nt-1] += neutmul[counter];
00227 }
00228 }
00229 }
00230 }
00231 }
00232 for (i = 0; i < numSec; ++i) {
00233 if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i];
00234 if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i];
00235 }
00236
00237
00238 for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
00239 for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
00240 counter = -1;
00241 for (npos = 1; npos < (numSec/3); ++npos) {
00242 nneg = npos-1;
00243 for (nzero = 0; nzero < numSec/3; ++nzero) {
00244 if (++counter < numMulA) {
00245 nt = npos+nneg+nzero;
00246 if (nt > 1 && nt <= numSec) {
00247 protmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00248 protnormA[nt-1] += protmulA[counter];
00249 }
00250 }
00251 }
00252 }
00253 for (i = 0; i < numMulA; ++i) neutmulA[i] = 0.0;
00254 for (i = 0; i < numSec; ++i) neutnormA[i] = 0.0;
00255 counter = -1;
00256 for (npos = 0; npos < numSec/3; ++npos) {
00257 nneg = npos;
00258 for (nzero = 0; nzero < numSec/3; ++nzero) {
00259 if (++counter < numMulA) {
00260 nt = npos+nneg+nzero;
00261 if (nt > 1 && nt <= numSec) {
00262 neutmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00263 neutnormA[nt-1] += neutmulA[counter];
00264 }
00265 }
00266 }
00267 }
00268 for (i = 0; i < numSec; ++i) {
00269 if (protnormA[i] > 0.0) protnormA[i] = 1.0/protnormA[i];
00270 if (neutnormA[i] > 0.0) neutnormA[i] = 1.0/neutnormA[i];
00271 }
00272 }
00273
00274 const G4double expxu = 82.;
00275 const G4double expxl = -expxu;
00276 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00277 G4ParticleDefinition *aProton = G4Proton::Proton();
00278 G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
00279 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00280
00281
00282
00283
00284 const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
00285 0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
00286 0.39,0.36,0.33,0.10,0.01};
00287 G4int iplab = G4int( pOriginal/GeV*10.0 );
00288 if( iplab > 9 )iplab = G4int( (pOriginal/GeV- 1.0)*5.0 ) + 10;
00289 if( iplab > 14 )iplab = G4int( pOriginal/GeV- 2.0 ) + 15;
00290 if( iplab > 22 )iplab = G4int( (pOriginal/GeV-10.0)/10.0 ) + 23;
00291 if( iplab > 24 )iplab = 24;
00292 if (G4UniformRand() > anhl[iplab]) {
00293 if (availableEnergy <= aPiPlus->GetPDGMass()/MeV) {
00294 quasiElastic = true;
00295 return;
00296 }
00297 G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
00298 const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
00299 G4double w0, wp, wt, wm;
00300 if ((availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) ) {
00301
00302
00303 npos = nneg = nzero = 0;
00304 if (targetParticle.GetDefinition() == aProton) {
00305 test = std::exp(std::min(expxu,
00306 std::max( expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
00307 w0 = test;
00308 wp = test;
00309 if (G4UniformRand() < w0/(w0+wp) ) nzero = 1;
00310 else npos = 1;
00311 } else {
00312
00313 test = std::exp(std::min(expxu,
00314 std::max( expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
00315 w0 = test;
00316 wp = test;
00317 test = std::exp(std::min(expxu,
00318 std::max( expxl, -(-1.0+b[1])*(-1.0+b[1])/(2.0*c*c) ) ) );
00319 wm = test;
00320 wt = w0+wp+wm;
00321 wp += w0;
00322 G4double ran = G4UniformRand();
00323 if( ran < w0/wt )
00324 nzero = 1;
00325 else if( ran < wp/wt )
00326 npos = 1;
00327 else
00328 nneg = 1;
00329 }
00330 } else {
00331
00332 G4double n, anpn;
00333 GetNormalizationConstant( availableEnergy, n, anpn );
00334 G4double ran = G4UniformRand();
00335 G4double dum, excs = 0.0;
00336 if (targetParticle.GetDefinition() == aProton) {
00337 counter = -1;
00338 for (npos = 0; npos < numSec/3 && ran>=excs; ++npos) {
00339 for (nneg = std::max(0,npos-2); nneg <= npos && ran >= excs; ++nneg) {
00340 for (nzero = 0; nzero < numSec/3 && ran >= excs; ++nzero) {
00341 if (++counter < numMul) {
00342 nt = npos + nneg + nzero;
00343 if (nt > 0) {
00344 test = std::exp(std::min(expxu,
00345 std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00346 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00347 if( std::fabs(dum) < 1.0 ) {
00348 if( test >= 1.0e-10 )excs += dum*test;
00349 }
00350 else
00351 excs += dum*test;
00352 }
00353 }
00354 }
00355 }
00356 }
00357 if (ran >= excs) {
00358
00359 quasiElastic = true;
00360 return;
00361 }
00362 npos--; nneg--; nzero--;
00363
00364 } else {
00365
00366 counter = -1;
00367 for (npos = 0; npos < numSec/3 && ran >= excs; ++npos) {
00368 for (nneg = std::max(0,npos-1); nneg <= (npos+1) && ran >= excs; ++nneg) {
00369 for (nzero = 0; nzero < numSec/3 && ran>=excs; ++nzero) {
00370 if (++counter < numMul) {
00371 nt = npos+nneg+nzero;
00372 if ((nt >= 1) && (nt <= numSec) ) {
00373 test = std::exp(std::min(expxu,
00374 std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00375 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00376 if (std::fabs(dum) < 1.0) {
00377 if (test >= 1.0e-10) excs += dum*test;
00378 }
00379 else
00380 excs += dum*test;
00381 }
00382 }
00383 }
00384 }
00385 }
00386 if (ran >= excs) {
00387
00388 quasiElastic = true;
00389 return;
00390 }
00391 npos--; nneg--; nzero--;
00392 }
00393 }
00394
00395 if (targetParticle.GetDefinition() == aProton) {
00396 switch (npos-nneg)
00397 {
00398 case 1:
00399 if( G4UniformRand() < 0.5 )
00400 {
00401 currentParticle.SetDefinitionAndUpdateE( anAntiProton );
00402 incidentHasChanged = true;
00403 }
00404 else
00405 {
00406 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00407 targetHasChanged = true;
00408 }
00409 break;
00410 case 2:
00411 currentParticle.SetDefinitionAndUpdateE( anAntiProton );
00412 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00413 incidentHasChanged = true;
00414 targetHasChanged = true;
00415 break;
00416 default:
00417 break;
00418 }
00419 } else {
00420
00421 switch (npos-nneg)
00422 {
00423 case 0:
00424 if (G4UniformRand() < 0.33) {
00425 currentParticle.SetDefinitionAndUpdateE( anAntiProton );
00426 targetParticle.SetDefinitionAndUpdateE( aProton );
00427 incidentHasChanged = true;
00428 targetHasChanged = true;
00429 }
00430 break;
00431 case 1:
00432 currentParticle.SetDefinitionAndUpdateE( anAntiProton );
00433 incidentHasChanged = true;
00434 break;
00435 default:
00436 targetParticle.SetDefinitionAndUpdateE( aProton );
00437 targetHasChanged = true;
00438 break;
00439 }
00440 }
00441 } else {
00442
00443 if (centerofmassEnergy <= 2*aPiPlus->GetPDGMass()/MeV) {
00444 quasiElastic = true;
00445 return;
00446 }
00447
00448
00449 G4double n, anpn;
00450 GetNormalizationConstant(-centerofmassEnergy, n, anpn);
00451 G4double ran = G4UniformRand();
00452 G4double dum, excs = 0.0;
00453 if (targetParticle.GetDefinition() == aProton) {
00454 counter = -1;
00455 for (npos = 1; (npos < numSec/3) && (ran>=excs); ++npos) {
00456 nneg = npos-1;
00457 for (nzero = 0; (nzero < numSec/3) && (ran>=excs); ++nzero) {
00458 if (++counter < numMulA) {
00459 nt = npos+nneg+nzero;
00460 if (nt > 1 && nt <= numSec) {
00461 test = std::exp(std::min(expxu,
00462 std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00463 dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
00464 if (std::fabs(dum) < 1.0)
00465 {
00466 if( test >= 1.0e-10 )excs += dum*test;
00467 }
00468 else
00469 excs += dum*test;
00470 }
00471 }
00472 }
00473 }
00474 } else {
00475
00476 counter = -1;
00477 for (npos = 0; (npos < numSec/3) && (ran >= excs); ++npos) {
00478 nneg = npos;
00479 for (nzero = 0; (nzero < numSec/3) && (ran >= excs); ++nzero) {
00480 if (++counter < numMulA) {
00481 nt = npos+nneg+nzero;
00482 if ((nt > 1) && (nt <= numSec) ) {
00483 test = std::exp(std::min(expxu,
00484 std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00485 dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
00486 if (std::fabs(dum) < 1.0) {
00487 if (test >= 1.0e-10 )excs += dum*test;
00488 }
00489 else
00490 excs += dum*test;
00491 }
00492 }
00493 }
00494 }
00495 }
00496 if (ran >= excs) {
00497
00498 quasiElastic = true;
00499 return;
00500 }
00501 npos--; nzero--;
00502 currentParticle.SetMass( 0.0 );
00503 targetParticle.SetMass( 0.0 );
00504 }
00505
00506 while(npos+nneg+nzero < 3) nzero++;
00507 SetUpPions(npos, nneg, nzero, vec, vecLen);
00508 return;
00509 }