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