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