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