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