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