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