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
00035
00036
00037
00038 #include "G4HEAntiKaonZeroInelastic.hh"
00039 #include "globals.hh"
00040 #include "G4ios.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043
00044 G4HEAntiKaonZeroInelastic::G4HEAntiKaonZeroInelastic(const G4String& name)
00045 : G4HEInelastic(name)
00046 {
00047 vecLength = 0;
00048 theMinEnergy = 20*GeV;
00049 theMaxEnergy = 10*TeV;
00050 MAXPART = 2048;
00051 verboseLevel = 0;
00052 G4cout << "WARNING: model G4HEAntiKaonZeroInelastic is being deprecated and will\n"
00053 << "disappear in Geant4 version 10.0" << G4endl;
00054 }
00055
00056
00057 void G4HEAntiKaonZeroInelastic::ModelDescription(std::ostream& outFile) const
00058 {
00059 outFile << "G4HEAntiKaonZeroInelastic is one of the High Energy\n"
00060 << "Parameterized (HEP) models used to implement inelastic\n"
00061 << "anti-K0 scattering from nuclei. It is a re-engineered\n"
00062 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00063 << "initial collision products into backward- and forward-going\n"
00064 << "clusters which are then decayed into final state hadrons.\n"
00065 << "The model does not conserve energy on an event-by-event\n"
00066 << "basis. It may be applied to anti-K0 with initial energies\n"
00067 << "above 20 GeV.\n";
00068 }
00069
00070
00071 G4HadFinalState*
00072 G4HEAntiKaonZeroInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00073 G4Nucleus& targetNucleus)
00074 {
00075 G4HEVector* pv = new G4HEVector[MAXPART];
00076 const G4HadProjectile* aParticle = &aTrack;
00077 const G4double atomicWeight = targetNucleus.GetA_asInt();
00078 const G4double atomicNumber = targetNucleus.GetZ_asInt();
00079 G4HEVector incidentParticle(aParticle);
00080
00081 G4int incidentCode = incidentParticle.getCode();
00082 G4double incidentMass = incidentParticle.getMass();
00083 G4double incidentTotalEnergy = incidentParticle.getEnergy();
00084
00085
00086
00087
00088 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00089
00090 if (incidentKineticEnergy < 1.)
00091 G4cout << "GHEAntiKaonZeroInelastic: incident energy < 1 GeV" << G4endl;
00092
00093 if (verboseLevel > 1) {
00094 G4cout << "G4HEAntiKaonZeroInelastic::ApplyYourself" << G4endl;
00095 G4cout << "incident particle " << incidentParticle.getName()
00096 << "mass " << incidentMass
00097 << "kinetic energy " << incidentKineticEnergy
00098 << G4endl;
00099 G4cout << "target material with (A,Z) = ("
00100 << atomicWeight << "," << atomicNumber << ")" << G4endl;
00101 }
00102
00103 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
00104 atomicWeight, atomicNumber);
00105 if (verboseLevel > 1)
00106 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00107
00108 incidentKineticEnergy -= inelasticity;
00109
00110 G4double excitationEnergyGNP = 0.;
00111 G4double excitationEnergyDTA = 0.;
00112
00113 G4double excitation = NuclearExcitation(incidentKineticEnergy,
00114 atomicWeight, atomicNumber,
00115 excitationEnergyGNP,
00116 excitationEnergyDTA);
00117 if (verboseLevel > 1)
00118 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
00119 << excitationEnergyDTA << G4endl;
00120
00121
00122 incidentKineticEnergy -= excitation;
00123 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00124
00125
00126
00127
00128 G4HEVector targetParticle;
00129 if (G4UniformRand() < atomicNumber/atomicWeight) {
00130 targetParticle.setDefinition("Proton");
00131 } else {
00132 targetParticle.setDefinition("Neutron");
00133 }
00134
00135 G4double targetMass = targetParticle.getMass();
00136 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00137 + targetMass*targetMass
00138 + 2.0*targetMass*incidentTotalEnergy);
00139 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00140
00141 vecLength = 0;
00142
00143 if (verboseLevel > 1)
00144 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00145 << incidentCode << G4endl;
00146
00147 G4bool successful = false;
00148
00149 G4bool inElastic = true;
00150 FirstIntInCasAntiKaonZero(inElastic, availableEnergy, pv, vecLength,
00151 incidentParticle, targetParticle );
00152
00153 if (verboseLevel > 1)
00154 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00155
00156 if ((vecLength > 0) && (availableEnergy > 1.))
00157 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00158 pv, vecLength,
00159 incidentParticle, targetParticle);
00160 HighEnergyCascading(successful, pv, vecLength,
00161 excitationEnergyGNP, excitationEnergyDTA,
00162 incidentParticle, targetParticle,
00163 atomicWeight, atomicNumber);
00164 if (!successful)
00165 HighEnergyClusterProduction(successful, pv, vecLength,
00166 excitationEnergyGNP, excitationEnergyDTA,
00167 incidentParticle, targetParticle,
00168 atomicWeight, atomicNumber);
00169 if (!successful)
00170 MediumEnergyCascading(successful, pv, vecLength,
00171 excitationEnergyGNP, excitationEnergyDTA,
00172 incidentParticle, targetParticle,
00173 atomicWeight, atomicNumber);
00174 if (!successful)
00175 MediumEnergyClusterProduction(successful, pv, vecLength,
00176 excitationEnergyGNP, excitationEnergyDTA,
00177 incidentParticle, targetParticle,
00178 atomicWeight, atomicNumber);
00179 if (!successful)
00180 QuasiElasticScattering(successful, pv, vecLength,
00181 excitationEnergyGNP, excitationEnergyDTA,
00182 incidentParticle, targetParticle,
00183 atomicWeight, atomicNumber);
00184 if (!successful)
00185 ElasticScattering(successful, pv, vecLength,
00186 incidentParticle,
00187 atomicWeight, atomicNumber);
00188
00189 if (!successful)
00190 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00191 << G4endl;
00192
00193 FillParticleChange(pv, vecLength);
00194 delete [] pv;
00195 theParticleChange.SetStatusChange(stopAndKill);
00196 return &theParticleChange;
00197 }
00198
00199
00200 void
00201 G4HEAntiKaonZeroInelastic::FirstIntInCasAntiKaonZero(G4bool& inElastic,
00202 const G4double availableEnergy,
00203 G4HEVector pv[],
00204 G4int &vecLen,
00205 const G4HEVector& incidentParticle,
00206 const G4HEVector& targetParticle)
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 {
00217 static const G4double expxu = 82.;
00218 static const G4double expxl = -expxu;
00219
00220 static const G4double protb = 0.7;
00221 static const G4double neutb = 0.7;
00222 static const G4double c = 1.25;
00223
00224 static const G4int numMul = 1200;
00225 static const G4int numSec = 60;
00226
00227 G4int neutronCode = Neutron.getCode();
00228 G4int protonCode = Proton.getCode();
00229 G4int kaonMinusCode = KaonMinus.getCode();
00230 G4int kaonZeroCode = KaonZero.getCode();
00231 G4int antiKaonZeroCode = AntiKaonZero.getCode();
00232
00233 G4int targetCode = targetParticle.getCode();
00234 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00235
00236 static G4bool first = true;
00237 static G4double protmul[numMul], protnorm[numSec];
00238 static G4double neutmul[numMul], neutnorm[numSec];
00239
00240
00241
00242
00243 G4int i, counter, nt, npos, nneg, nzero;
00244
00245 if (first) {
00246 first = false;
00247 for (i = 0; i < numMul; i++) protmul[i] = 0.0;
00248 for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
00249 counter = -1;
00250 for (npos = 0; npos < (numSec/3); npos++) {
00251 for (nneg = std::max(0,npos-2); nneg <= npos; nneg++) {
00252 for (nzero = 0; nzero < numSec/3; nzero++) {
00253 if (++counter < numMul) {
00254 nt = npos+nneg+nzero;
00255 if ((nt > 0) && (nt <= numSec) ) {
00256 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00257 protnorm[nt-1] += protmul[counter];
00258 }
00259 }
00260 }
00261 }
00262 }
00263
00264 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00265 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00266 counter = -1;
00267 for (npos = 0; npos < numSec/3; npos++) {
00268 for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
00269 for (nzero = 0; nzero < numSec/3; nzero++) {
00270 if (++counter < numMul) {
00271 nt = npos+nneg+nzero;
00272 if ((nt > 0) && (nt <= numSec) ) {
00273 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00274 neutnorm[nt-1] += neutmul[counter];
00275 }
00276 }
00277 }
00278 }
00279 }
00280 for (i = 0; i < numSec; i++) {
00281 if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i];
00282 if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i];
00283 }
00284 }
00285
00286 pv[0] = incidentParticle;
00287 pv[1] = targetParticle;
00288 vecLen = 2;
00289
00290 if (!inElastic || (availableEnergy <= PionPlus.getMass())) return;
00291
00292
00293 npos = 0, nneg = 0, nzero = 0;
00294 G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
00295 G4int iplab = G4int( incidentTotalMomentum*5.);
00296 if ( (iplab < 10) && (G4UniformRand() < cech[iplab]) ) {
00297 G4int ipl = std::min(19, G4int( incidentTotalMomentum*5.));
00298 G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
00299 0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
00300 if (G4UniformRand() < cnk0[ipl] ) {
00301 if (targetCode == protonCode) {
00302 return;
00303 } else {
00304 pv[0] = KaonMinus;
00305 pv[1] = Proton;
00306 return;
00307 }
00308 }
00309
00310 G4double ran = G4UniformRand();
00311 if (targetCode == protonCode) {
00312 if (ran < 0.25) {
00313 } else if (ran < 0.50) {
00314 pv[0] = PionPlus;
00315 pv[1] = SigmaZero;
00316 } else if (ran < 0.75) {
00317 } else {
00318 pv[0] = PionPlus;
00319 pv[1] = Lambda;
00320 }
00321 } else {
00322 if( ran < 0.25 )
00323 {
00324 pv[0] = PionMinus;
00325 pv[1] = SigmaPlus;
00326 }
00327 else if (ran < 0.50)
00328 {
00329 pv[0] = PionZero;
00330 pv[1] = SigmaZero;
00331 }
00332 else if (ran < 0.75)
00333 {
00334 pv[0] = PionPlus;
00335 pv[1] = SigmaMinus;
00336 }
00337 else
00338 {
00339 pv[0] = PionZero;
00340 pv[1] = Lambda;
00341 }
00342 }
00343 return;
00344
00345 } else {
00346
00347 G4double aleab = std::log(availableEnergy);
00348 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00349 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00350
00351
00352
00353 G4double test, dum, anpn = 0.0;
00354
00355 for (nt=1; nt<=numSec; nt++) {
00356 test = std::exp(std::min(expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00357 dum = pi*nt/(2.0*n*n);
00358 if (std::fabs(dum) < 1.0) {
00359 if( test >= 1.0e-10 )anpn += dum*test;
00360 } else {
00361 anpn += dum*test;
00362 }
00363 }
00364
00365 G4double ran = G4UniformRand();
00366 G4double excs = 0.0;
00367 if (targetCode == protonCode) {
00368 counter = -1;
00369 for (npos=0; npos<numSec/3; npos++) {
00370 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
00371 for (nzero=0; nzero<numSec/3; nzero++) {
00372 if (++counter < numMul) {
00373 nt = npos+nneg+nzero;
00374 if ((nt>0) && (nt<=numSec) ) {
00375 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00376 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00377
00378 if (std::fabs(dum) < 1.0) {
00379 if( test >= 1.0e-10 )excs += dum*test;
00380 } else {
00381 excs += dum*test;
00382 }
00383
00384 if (ran < excs) goto outOfLoop;
00385 }
00386 }
00387 }
00388 }
00389 }
00390
00391 inElastic = false;
00392 return;
00393
00394 } else {
00395 counter = -1;
00396 for (npos=0; npos<numSec/3; npos++) {
00397 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00398 for (nzero=0; nzero<numSec/3; nzero++) {
00399 if (++counter < numMul) {
00400 nt = npos+nneg+nzero;
00401 if( (nt>=1) && (nt<=numSec) ) {
00402 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00403 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00404
00405 if (std::fabs(dum) < 1.0) {
00406 if( test >= 1.0e-10 )excs += dum*test;
00407 } else {
00408 excs += dum*test;
00409 }
00410
00411 if (ran < excs) goto outOfLoop;
00412 }
00413 }
00414 }
00415 }
00416 }
00417
00418 inElastic = false;
00419 return;
00420 }
00421 }
00422
00423 outOfLoop:
00424
00425 if (targetCode == protonCode) {
00426 if (npos == nneg) {
00427
00428 } else if (npos == (1+nneg)) {
00429 if (G4UniformRand() < 0.5) {
00430 pv[0] = KaonMinus;
00431 } else {
00432 pv[1] = Neutron;
00433 }
00434 } else {
00435 pv[0] = KaonMinus;
00436 pv[1] = Neutron;
00437 }
00438
00439 } else {
00440 if( npos == nneg)
00441 {
00442 if( G4UniformRand() < 0.75)
00443 {
00444 }
00445 else
00446 {
00447 pv[0] = KaonMinus;
00448 pv[1] = Proton;
00449 }
00450 }
00451 else if ( npos == (1+nneg))
00452 {
00453 pv[0] = KaonMinus;
00454 }
00455 else
00456 {
00457 pv[1] = Proton;
00458 }
00459 }
00460
00461 if (G4UniformRand() < 0.5) {
00462 if (((pv[0].getCode() == kaonMinusCode)
00463 && (pv[1].getCode() == neutronCode) )
00464 || ((pv[0].getCode() == kaonZeroCode)
00465 && (pv[1].getCode() == protonCode) )
00466 || ((pv[0].getCode() == antiKaonZeroCode)
00467 && (pv[1].getCode() == protonCode) ) ) {
00468
00469 G4double ran = G4UniformRand();
00470 if (pv[1].getCode() == protonCode) {
00471 if (ran < 0.68) {
00472 pv[0] = PionPlus;
00473 pv[1] = Lambda;
00474 } else if (ran < 0.84) {
00475 pv[0] = PionZero;
00476 pv[1] = SigmaPlus;
00477 } else {
00478 pv[0] = PionPlus;
00479 pv[1] = SigmaZero;
00480 }
00481 } else {
00482 if(ran < 0.68)
00483 {
00484 pv[0] = PionMinus;
00485 pv[1] = Lambda;
00486 }
00487 else if (ran < 0.84)
00488 {
00489 pv[0] = PionMinus;
00490 pv[1] = SigmaZero;
00491 }
00492 else
00493 {
00494 pv[0] = PionZero;
00495 pv[1] = SigmaMinus;
00496 }
00497 }
00498 } else {
00499 G4double ran = G4UniformRand();
00500 if (ran < 0.67)
00501 {
00502 pv[0] = PionZero;
00503 pv[1] = Lambda;
00504 }
00505 else if (ran < 0.78)
00506 {
00507 pv[0] = PionMinus;
00508 pv[1] = SigmaPlus;
00509 }
00510 else if (ran < 0.89)
00511 {
00512 pv[0] = PionZero;
00513 pv[1] = SigmaZero;
00514 }
00515 else
00516 {
00517 pv[0] = PionPlus;
00518 pv[1] = SigmaMinus;
00519 }
00520 }
00521 }
00522
00523 nt = npos + nneg + nzero;
00524 while (nt > 0) {
00525 G4double ran = G4UniformRand();
00526 if (ran < (G4double)npos/nt) {
00527 if (npos > 0) {
00528 pv[vecLen++] = PionPlus;
00529 npos--;
00530 }
00531 } else if (ran < (G4double)(npos+nneg)/nt) {
00532 if (nneg > 0) {
00533 pv[vecLen++] = PionMinus;
00534 nneg--;
00535 }
00536 } else {
00537 if (nzero > 0) {
00538 pv[vecLen++] = PionZero;
00539 nzero--;
00540 }
00541 }
00542 nt = npos + nneg + nzero;
00543 }
00544
00545 if (verboseLevel > 1) {
00546 G4cout << "Particles produced: " ;
00547 G4cout << pv[0].getName() << " " ;
00548 G4cout << pv[1].getName() << " " ;
00549 for (i=2; i < vecLen; i++) {
00550 G4cout << pv[i].getName() << " " ;
00551 }
00552 G4cout << G4endl;
00553 }
00554
00555 return;
00556 }
00557