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 "G4HEKaonZeroShortInelastic.hh"
00039 #include "globals.hh"
00040 #include "G4ios.hh"
00041 #include "G4PhysicalConstants.hh"
00042
00043 void G4HEKaonZeroShortInelastic::ModelDescription(std::ostream& outFile) const
00044 {
00045 outFile << "G4HEKaonZeroShortInelastic is one of the High Energy\n"
00046 << "Parameterized (HEP) models used to implement inelastic\n"
00047 << "K0S 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 K0S with initial energies\n"
00053 << "above 20 GeV.\n";
00054 }
00055
00056
00057 G4HadFinalState*
00058 G4HEKaonZeroShortInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00059 G4Nucleus& targetNucleus)
00060 {
00061 G4HEVector* pv = new G4HEVector[MAXPART];
00062 const G4HadProjectile* aParticle = &aTrack;
00063 const G4double atomicWeight = targetNucleus.GetA_asInt();
00064 const G4double atomicNumber = targetNucleus.GetZ_asInt();
00065 G4HEVector incidentParticle(aParticle);
00066
00067 G4int incidentCode = incidentParticle.getCode();
00068 G4double incidentMass = incidentParticle.getMass();
00069 G4double incidentTotalEnergy = incidentParticle.getEnergy();
00070
00071
00072
00073
00074 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00075
00076 if(incidentKineticEnergy < 1)
00077 G4cout << "GHEKaonZeroShortInelastic: incident energy < 1 GeV" << G4endl;
00078
00079 if(verboseLevel > 1) {
00080 G4cout << "G4HEKaonZeroShortInelastic::ApplyYourself" << G4endl;
00081 G4cout << "incident particle " << incidentParticle.getName()
00082 << "mass " << incidentMass
00083 << "kinetic energy " << incidentKineticEnergy
00084 << G4endl;
00085 G4cout << "target material with (A,Z) = ("
00086 << atomicWeight << "," << atomicNumber << ")" << G4endl;
00087 }
00088
00089 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
00090 atomicWeight, atomicNumber);
00091 if(verboseLevel > 1)
00092 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00093
00094 incidentKineticEnergy -= inelasticity;
00095
00096 G4double excitationEnergyGNP = 0.;
00097 G4double excitationEnergyDTA = 0.;
00098
00099 G4double excitation = NuclearExcitation(incidentKineticEnergy,
00100 atomicWeight, atomicNumber,
00101 excitationEnergyGNP,
00102 excitationEnergyDTA);
00103 if(verboseLevel > 1)
00104 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
00105 << excitationEnergyDTA << G4endl;
00106
00107 incidentKineticEnergy -= excitation;
00108 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00109
00110
00111
00112
00113 G4HEVector targetParticle;
00114 if(G4UniformRand() < atomicNumber/atomicWeight) {
00115 targetParticle.setDefinition("Proton");
00116 } else {
00117 targetParticle.setDefinition("Neutron");
00118 }
00119
00120 G4double targetMass = targetParticle.getMass();
00121 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00122 + targetMass*targetMass
00123 + 2.0*targetMass*incidentTotalEnergy);
00124 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00125
00126 G4bool inElastic = true;
00127 vecLength = 0;
00128
00129 if(verboseLevel > 1)
00130 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00131 << incidentCode << G4endl;
00132
00133 G4bool successful = false;
00134
00135
00136 if (G4UniformRand() < 0.5)
00137 FirstIntInCasAntiKaonZero(inElastic, availableEnergy, pv, vecLength,
00138 incidentParticle, targetParticle );
00139 else
00140 FirstIntInCasKaonZero(inElastic, availableEnergy, pv, vecLength,
00141 incidentParticle, targetParticle, atomicWeight );
00142
00143
00144 if ((vecLength > 0) && (availableEnergy > 1.))
00145 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00146 pv, vecLength,
00147 incidentParticle, targetParticle);
00148
00149 HighEnergyCascading(successful, pv, vecLength,
00150 excitationEnergyGNP, excitationEnergyDTA,
00151 incidentParticle, targetParticle,
00152 atomicWeight, atomicNumber);
00153 if (!successful)
00154 HighEnergyClusterProduction(successful, pv, vecLength,
00155 excitationEnergyGNP, excitationEnergyDTA,
00156 incidentParticle, targetParticle,
00157 atomicWeight, atomicNumber);
00158 if (!successful)
00159 MediumEnergyCascading(successful, pv, vecLength,
00160 excitationEnergyGNP, excitationEnergyDTA,
00161 incidentParticle, targetParticle,
00162 atomicWeight, atomicNumber);
00163
00164 if (!successful)
00165 MediumEnergyClusterProduction(successful, pv, vecLength,
00166 excitationEnergyGNP, excitationEnergyDTA,
00167 incidentParticle, targetParticle,
00168 atomicWeight, atomicNumber);
00169 if (!successful)
00170 QuasiElasticScattering(successful, pv, vecLength,
00171 excitationEnergyGNP, excitationEnergyDTA,
00172 incidentParticle, targetParticle,
00173 atomicWeight, atomicNumber);
00174
00175 if (!successful)
00176 ElasticScattering(successful, pv, vecLength,
00177 incidentParticle,
00178 atomicWeight, atomicNumber);
00179
00180 if (!successful)
00181 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00182 << G4endl;
00183
00184
00185 G4int kcode;
00186 for (G4int i = 0; i < vecLength; i++) {
00187 kcode = pv[i].getCode();
00188 if (kcode == KaonZero.getCode() || kcode == AntiKaonZero.getCode()) {
00189 if (G4UniformRand() < 0.5)
00190 pv[i] = KaonZeroShort;
00191 else
00192 pv[i] = KaonZeroLong;
00193 }
00194 }
00195
00196
00197
00198 FillParticleChange(pv, vecLength);
00199 delete [] pv;
00200 theParticleChange.SetStatusChange(stopAndKill);
00201 return &theParticleChange;
00202 }
00203
00204
00205 void
00206 G4HEKaonZeroShortInelastic::FirstIntInCasKaonZero(G4bool& inElastic,
00207 const G4double availableEnergy,
00208 G4HEVector pv[],
00209 G4int& vecLen,
00210 const G4HEVector& incidentParticle,
00211 const G4HEVector& targetParticle,
00212 const G4double atomicWeight)
00213
00214
00215
00216
00217
00218
00219
00220
00221 {
00222 static const G4double expxu = 82.;
00223 static const G4double expxl = -expxu;
00224
00225 static const G4double protb = 0.7;
00226 static const G4double neutb = 0.7;
00227 static const G4double c = 1.25;
00228
00229 static const G4int numMul = 1200;
00230 static const G4int numSec = 60;
00231
00232 G4int neutronCode = Neutron.getCode();
00233 G4int protonCode = Proton.getCode();
00234
00235 G4int targetCode = targetParticle.getCode();
00236 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00237
00238 static G4bool first = true;
00239 static G4double protmul[numMul], protnorm[numSec];
00240 static G4double neutmul[numMul], neutnorm[numSec];
00241
00242
00243
00244
00245 G4int i, counter, nt, npos, nneg, nzero;
00246
00247 if (first) {
00248
00249 first = false;
00250 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
00251 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00252 counter = -1;
00253 for (npos=0; npos<(numSec/3); npos++) {
00254 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00255 for (nzero=0; nzero<numSec/3; nzero++) {
00256 if (++counter < numMul) {
00257 nt = npos+nneg+nzero;
00258 if( (nt>0) && (nt<=numSec) ) {
00259 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c) ;
00260 protnorm[nt-1] += protmul[counter];
00261 }
00262 }
00263 }
00264 }
00265 }
00266
00267 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00268 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00269 counter = -1;
00270 for (npos=0; npos<numSec/3; npos++) {
00271 for (nneg=npos; nneg<=(npos+2); nneg++) {
00272 for (nzero=0; nzero<numSec/3; nzero++) {
00273 if (++counter < numMul) {
00274 nt = npos+nneg+nzero;
00275 if( (nt>0) && (nt<=numSec) ) {
00276 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00277 neutnorm[nt-1] += neutmul[counter];
00278 }
00279 }
00280 }
00281 }
00282 }
00283
00284 for (i=0; i<numSec; i++) {
00285 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00286 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00287 }
00288 }
00289
00290
00291
00292
00293 pv[0] = incidentParticle;
00294 pv[1] = targetParticle;
00295 vecLen = 2;
00296
00297 if( !inElastic ) {
00298
00299 if( targetCode == protonCode) {
00300 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07};
00301 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*5. ) );
00302 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42)) {
00303
00304 pv[0] = KaonPlus;
00305 pv[1] = Neutron;
00306 }
00307 }
00308 return;
00309 } else if (availableEnergy <= PionPlus.getMass()) {
00310 return;
00311 }
00312
00313
00314
00315 npos = 0, nneg = 0, nzero = 0;
00316 G4double eab = availableEnergy;
00317 G4int ieab = G4int( eab*5.0 );
00318
00319 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
00320 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab])) {
00321
00322
00323 G4double w0, wp, wm, wt, ran;
00324 if (targetCode == neutronCode) {
00325
00326 w0 = - sqr(1.+protb)/(2.*c*c);
00327 w0 = std::exp(w0);
00328 wm = - sqr(-1.+protb)/(2.*c*c);
00329 wm = std::exp(wm);
00330 w0 = w0/2.;
00331 wm = wm*1.5;
00332 if (G4UniformRand() < w0/(w0+wm) ) {
00333 npos = 0;
00334 nneg = 0;
00335 nzero = 1;
00336 } else {
00337 npos = 0;
00338 nneg = 1;
00339 nzero = 0;
00340 }
00341
00342 } else {
00343
00344 w0 = -sqr(1.+neutb)/(2.*c*c);
00345 wp = w0 = std::exp(w0);
00346 wm = -sqr(-1.+neutb)/(2.*c*c);
00347 wm = std::exp(wm);
00348 wt = w0+wp+wm;
00349 wp = w0+wp;
00350 ran = G4UniformRand();
00351 if ( ran < w0/wt) {
00352 npos = 0;
00353 nneg = 0;
00354 nzero = 1;
00355 } else if (ran < wp/wt) {
00356 npos = 1;
00357 nneg = 0;
00358 nzero = 0;
00359 } else {
00360 npos = 0;
00361 nneg = 1;
00362 nzero = 0;
00363 }
00364 }
00365 } else {
00366
00367
00368 G4double aleab = std::log(availableEnergy);
00369 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00370 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00371
00372
00373
00374 G4double test, dum, anpn = 0.0;
00375
00376 for (nt=1; nt<=numSec; nt++) {
00377 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00378 dum = pi*nt/(2.0*n*n);
00379 if (std::fabs(dum) < 1.0) {
00380 if( test >= 1.0e-10 )anpn += dum*test;
00381 } else {
00382 anpn += dum*test;
00383 }
00384 }
00385
00386 G4double ran = G4UniformRand();
00387 G4double excs = 0.0;
00388 if( targetCode == protonCode )
00389 {
00390 counter = -1;
00391 for( npos=0; npos<numSec/3; npos++ )
00392 {
00393 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
00394 {
00395 for (nzero=0; nzero<numSec/3; nzero++) {
00396 if (++counter < numMul) {
00397 nt = npos+nneg+nzero;
00398 if ( (nt>0) && (nt<=numSec) ) {
00399 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00400 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00401 if (std::fabs(dum) < 1.0) {
00402 if( test >= 1.0e-10 )excs += dum*test;
00403 } else {
00404 excs += dum*test;
00405 }
00406 if (ran < excs) goto outOfLoop;
00407 }
00408 }
00409 }
00410 }
00411 }
00412
00413
00414 inElastic = false;
00415 return;
00416 }
00417 else
00418 {
00419 counter = -1;
00420 for( npos=0; npos<numSec/3; npos++ )
00421 {
00422 for( nneg=npos; nneg<=(npos+2); nneg++ )
00423 {
00424 for (nzero=0; nzero<numSec/3; nzero++) {
00425 if (++counter < numMul) {
00426 nt = npos+nneg+nzero;
00427 if ( (nt>=1) && (nt<=numSec) ) {
00428 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00429 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00430 if (std::fabs(dum) < 1.0) {
00431 if( test >= 1.0e-10 )excs += dum*test;
00432 } else {
00433 excs += dum*test;
00434 }
00435 if (ran < excs) goto outOfLoop;
00436 }
00437 }
00438 }
00439 }
00440 }
00441
00442 inElastic = false;
00443 return;
00444 }
00445 }
00446 outOfLoop:
00447
00448 if( targetCode == neutronCode)
00449 {
00450 if( npos == nneg)
00451 {
00452 }
00453 else if (npos == (nneg-1))
00454 {
00455 if( G4UniformRand() < 0.5)
00456 {
00457 pv[0] = KaonPlus;
00458 }
00459 else
00460 {
00461 pv[1] = Proton;
00462 }
00463 }
00464 else
00465 {
00466 pv[0] = KaonPlus;
00467 pv[1] = Proton;
00468 }
00469 }
00470 else
00471 {
00472 if( npos == nneg )
00473 {
00474 if( G4UniformRand() < 0.25)
00475 {
00476 pv[0] = KaonPlus;
00477 pv[1] = Neutron;
00478 }
00479 else
00480 {
00481 }
00482 }
00483 else if ( npos == (nneg+1))
00484 {
00485 pv[1] = Neutron;
00486 }
00487 else
00488 {
00489 pv[0] = KaonPlus;
00490 }
00491 }
00492
00493 nt = npos + nneg + nzero;
00494 while (nt > 0) {
00495 G4double ran = G4UniformRand();
00496 if (ran < (G4double)npos/nt) {
00497 if (npos > 0) {
00498 pv[vecLen++] = PionPlus;
00499 npos--;
00500 }
00501 } else if ( ran < (G4double)(npos+nneg)/nt) {
00502 if (nneg > 0) {
00503 pv[vecLen++] = PionMinus;
00504 nneg--;
00505 }
00506 } else {
00507 if (nzero > 0) {
00508 pv[vecLen++] = PionZero;
00509 nzero--;
00510 }
00511 }
00512 nt = npos + nneg + nzero;
00513 }
00514
00515 if (verboseLevel > 1) {
00516 G4cout << "Particles produced: " ;
00517 G4cout << pv[0].getName() << " " ;
00518 G4cout << pv[1].getName() << " " ;
00519 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
00520 G4cout << G4endl;
00521 }
00522
00523 return;
00524 }
00525
00526
00527 void
00528 G4HEKaonZeroShortInelastic::FirstIntInCasAntiKaonZero(G4bool& inElastic,
00529 const G4double availableEnergy,
00530 G4HEVector pv[],
00531 G4int& vecLen,
00532 const G4HEVector& incidentParticle,
00533 const G4HEVector& targetParticle)
00534
00535
00536
00537
00538
00539
00540
00541
00542 {
00543 static const G4double expxu = 82.;
00544 static const G4double expxl = -expxu;
00545
00546 static const G4double protb = 0.7;
00547 static const G4double neutb = 0.7;
00548 static const G4double c = 1.25;
00549
00550 static const G4int numMul = 1200;
00551 static const G4int numSec = 60;
00552
00553 G4int neutronCode = Neutron.getCode();
00554 G4int protonCode = Proton.getCode();
00555 G4int kaonMinusCode = KaonMinus.getCode();
00556 G4int kaonZeroCode = KaonZero.getCode();
00557 G4int antiKaonZeroCode = AntiKaonZero.getCode();
00558
00559 G4int targetCode = targetParticle.getCode();
00560 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00561
00562 static G4bool first = true;
00563 static G4double protmul[numMul], protnorm[numSec];
00564 static G4double neutmul[numMul], neutnorm[numSec];
00565
00566
00567
00568
00569 G4int i, counter, nt, npos, nneg, nzero;
00570
00571 if (first) {
00572
00573 first = false;
00574 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
00575 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00576 counter = -1;
00577 for(npos=0; npos<(numSec/3); npos++) {
00578 for(nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
00579 for(nzero=0; nzero<numSec/3; nzero++) {
00580 if(++counter < numMul) {
00581 nt = npos+nneg+nzero;
00582 if( (nt>0) && (nt<=numSec) ) {
00583 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c) ;
00584 protnorm[nt-1] += protmul[counter];
00585 }
00586 }
00587 }
00588 }
00589 }
00590
00591 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00592 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00593 counter = -1;
00594 for(npos=0; npos<numSec/3; npos++) {
00595 for(nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00596 for(nzero=0; nzero<numSec/3; nzero++) {
00597 if(++counter < numMul) {
00598 nt = npos+nneg+nzero;
00599 if( (nt>0) && (nt<=numSec) ) {
00600 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00601 neutnorm[nt-1] += neutmul[counter];
00602 }
00603 }
00604 }
00605 }
00606 }
00607
00608 for(i=0; i<numSec; i++) {
00609 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00610 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00611 }
00612 }
00613
00614
00615
00616 pv[0] = incidentParticle;
00617 pv[1] = targetParticle;
00618 vecLen = 2;
00619
00620 if (!inElastic || (availableEnergy <= PionPlus.getMass()))
00621 return;
00622
00623
00624
00625 npos = 0, nneg = 0, nzero = 0;
00626 G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
00627 G4int iplab = G4int( incidentTotalMomentum*5.);
00628 if( (iplab < 10) && (G4UniformRand() < cech[iplab]) ) {
00629 G4int ipl = std::min(19, G4int( incidentTotalMomentum*5.));
00630 G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
00631 0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
00632 if(G4UniformRand() < cnk0[ipl]) {
00633 if(targetCode == protonCode) {
00634 return;
00635 } else {
00636 pv[0] = KaonMinus;
00637 pv[1] = Proton;
00638 return;
00639 }
00640 }
00641
00642 G4double ran = G4UniformRand();
00643 if(targetCode == protonCode) {
00644
00645
00646 if( ran < 0.25 ) {
00647 ;
00648 } else if (ran < 0.50) {
00649 pv[0] = PionPlus;
00650 pv[1] = SigmaZero;
00651 } else if (ran < 0.75) {
00652 ;
00653 } else {
00654 pv[0] = PionPlus;
00655 pv[1] = Lambda;
00656 }
00657 } else {
00658
00659
00660 if( ran < 0.25 ) {
00661 pv[0] = PionMinus;
00662 pv[1] = SigmaPlus;
00663 } else if (ran < 0.50) {
00664 pv[0] = PionZero;
00665 pv[1] = SigmaZero;
00666 } else if (ran < 0.75) {
00667 pv[0] = PionPlus;
00668 pv[1] = SigmaMinus;
00669 } else {
00670 pv[0] = PionZero;
00671 pv[1] = Lambda;
00672 }
00673 }
00674 return;
00675
00676 } else {
00677
00678
00679 G4double aleab = std::log(availableEnergy);
00680 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00681 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00682
00683
00684
00685 G4double test, dum, anpn = 0.0;
00686
00687 for (nt=1; nt<=numSec; nt++) {
00688 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00689 dum = pi*nt/(2.0*n*n);
00690 if (std::fabs(dum) < 1.0) {
00691 if( test >= 1.0e-10 )anpn += dum*test;
00692 } else {
00693 anpn += dum*test;
00694 }
00695 }
00696
00697 G4double ran = G4UniformRand();
00698 G4double excs = 0.0;
00699 if (targetCode == protonCode) {
00700 counter = -1;
00701 for (npos=0; npos<numSec/3; npos++) {
00702 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
00703 for (nzero=0; nzero<numSec/3; nzero++) {
00704 if (++counter < numMul) {
00705 nt = npos+nneg+nzero;
00706 if( (nt>0) && (nt<=numSec) ) {
00707 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00708 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00709
00710 if (std::fabs(dum) < 1.0) {
00711 if( test >= 1.0e-10 )excs += dum*test;
00712 } else {
00713 excs += dum*test;
00714 }
00715
00716 if (ran < excs) goto outOfLoop;
00717 }
00718 }
00719 }
00720 }
00721 }
00722
00723 inElastic = false;
00724 return;
00725
00726 } else {
00727 counter = -1;
00728 for (npos=0; npos<numSec/3; npos++) {
00729 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00730 for (nzero=0; nzero<numSec/3; nzero++) {
00731 if (++counter < numMul) {
00732 nt = npos+nneg+nzero;
00733 if( (nt>=1) && (nt<=numSec) ) {
00734 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00735 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00736
00737 if (std::fabs(dum) < 1.0) {
00738 if( test >= 1.0e-10 )excs += dum*test;
00739 } else {
00740 excs += dum*test;
00741 }
00742
00743 if (ran < excs) goto outOfLoop;
00744 }
00745 }
00746 }
00747 }
00748 }
00749
00750 inElastic = false;
00751 return;
00752 }
00753 }
00754 outOfLoop:
00755
00756 if( targetCode == protonCode)
00757 {
00758 if( npos == nneg)
00759 {
00760 }
00761 else if (npos == (1+nneg))
00762 {
00763 if( G4UniformRand() < 0.5)
00764 {
00765 pv[0] = KaonMinus;
00766 }
00767 else
00768 {
00769 pv[1] = Neutron;
00770 }
00771 }
00772 else
00773 {
00774 pv[0] = KaonMinus;
00775 pv[1] = Neutron;
00776 }
00777 }
00778 else
00779 {
00780 if( npos == nneg)
00781 {
00782 if( G4UniformRand() < 0.75)
00783 {
00784 }
00785 else
00786 {
00787 pv[0] = KaonMinus;
00788 pv[1] = Proton;
00789 }
00790 }
00791 else if ( npos == (1+nneg))
00792 {
00793 pv[0] = KaonMinus;
00794 }
00795 else
00796 {
00797 pv[1] = Proton;
00798 }
00799 }
00800
00801
00802 if( G4UniformRand() < 0.5 )
00803 {
00804 if( ( (pv[0].getCode() == kaonMinusCode)
00805 && (pv[1].getCode() == neutronCode) )
00806 || ( (pv[0].getCode() == kaonZeroCode)
00807 && (pv[1].getCode() == protonCode) )
00808 || ( (pv[0].getCode() == antiKaonZeroCode)
00809 && (pv[1].getCode() == protonCode) ) )
00810 {
00811 G4double ran = G4UniformRand();
00812 if( pv[1].getCode() == protonCode)
00813 {
00814 if(ran < 0.68)
00815 {
00816 pv[0] = PionPlus;
00817 pv[1] = Lambda;
00818 }
00819 else if (ran < 0.84)
00820 {
00821 pv[0] = PionZero;
00822 pv[1] = SigmaPlus;
00823 }
00824 else
00825 {
00826 pv[0] = PionPlus;
00827 pv[1] = SigmaZero;
00828 }
00829 }
00830 else
00831 {
00832 if(ran < 0.68)
00833 {
00834 pv[0] = PionMinus;
00835 pv[1] = Lambda;
00836 }
00837 else if (ran < 0.84)
00838 {
00839 pv[0] = PionMinus;
00840 pv[1] = SigmaZero;
00841 }
00842 else
00843 {
00844 pv[0] = PionZero;
00845 pv[1] = SigmaMinus;
00846 }
00847 }
00848 }
00849 else
00850 {
00851 G4double ran = G4UniformRand();
00852 if (ran < 0.67)
00853 {
00854 pv[0] = PionZero;
00855 pv[1] = Lambda;
00856 }
00857 else if (ran < 0.78)
00858 {
00859 pv[0] = PionMinus;
00860 pv[1] = SigmaPlus;
00861 }
00862 else if (ran < 0.89)
00863 {
00864 pv[0] = PionZero;
00865 pv[1] = SigmaZero;
00866 }
00867 else
00868 {
00869 pv[0] = PionPlus;
00870 pv[1] = SigmaMinus;
00871 }
00872 }
00873 }
00874
00875 nt = npos + nneg + nzero;
00876 while ( nt > 0) {
00877 G4double ran = G4UniformRand();
00878 if ( ran < (G4double)npos/nt) {
00879 if( npos > 0 ) {
00880 pv[vecLen++] = PionPlus;
00881 npos--;
00882 }
00883 } else if (ran < (G4double)(npos+nneg)/nt) {
00884 if( nneg > 0 ) {
00885 pv[vecLen++] = PionMinus;
00886 nneg--;
00887 }
00888 } else {
00889 if( nzero > 0 ) {
00890 pv[vecLen++] = PionZero;
00891 nzero--;
00892 }
00893 }
00894 nt = npos + nneg + nzero;
00895 }
00896
00897 if (verboseLevel > 1) {
00898 G4cout << "Particles produced: " ;
00899 G4cout << pv[0].getName() << " " ;
00900 G4cout << pv[1].getName() << " " ;
00901 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
00902 G4cout << G4endl;
00903 }
00904
00905 return;
00906 }