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