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