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