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 "G4HEKaonZeroLongInelastic.hh"
00039 #include "globals.hh"
00040 #include "G4ios.hh"
00041 #include "G4PhysicalConstants.hh"
00042
00043 void G4HEKaonZeroLongInelastic::ModelDescription(std::ostream& outFile) const
00044 {
00045 outFile << "G4HEKaonZeroLongInelastic is one of the High Energy\n"
00046 << "Parameterized (HEP) models used to implement inelastic\n"
00047 << "K0L 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 K0L with initial energies\n"
00053 << "above 20 GeV.\n";
00054 }
00055
00056
00057 G4HadFinalState*
00058 G4HEKaonZeroLongInelastic::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 << "GHEKaonZeroLongInelastic: incident energy < 1 GeV " << G4endl;
00078
00079 if(verboseLevel > 1) {
00080 G4cout << "G4HEKaonZeroLongInelastic::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 G4HEKaonZeroLongInelastic::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 G4HEKaonZeroLongInelastic::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 {
00544 static const G4double expxu = 82.;
00545 static const G4double expxl = -expxu;
00546
00547 static const G4double protb = 0.7;
00548 static const G4double neutb = 0.7;
00549 static const G4double c = 1.25;
00550
00551 static const G4int numMul = 1200;
00552 static const G4int numSec = 60;
00553
00554 G4int neutronCode = Neutron.getCode();
00555 G4int protonCode = Proton.getCode();
00556 G4int kaonMinusCode = KaonMinus.getCode();
00557 G4int kaonZeroCode = KaonZero.getCode();
00558 G4int antiKaonZeroCode = AntiKaonZero.getCode();
00559
00560 G4int targetCode = targetParticle.getCode();
00561 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00562
00563 static G4bool first = true;
00564 static G4double protmul[numMul], protnorm[numSec];
00565 static G4double neutmul[numMul], neutnorm[numSec];
00566
00567
00568
00569
00570 G4int i, counter, nt, npos, nneg, nzero;
00571
00572 if (first) {
00573
00574 first = false;
00575 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
00576 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00577 counter = -1;
00578 for(npos=0; npos<(numSec/3); npos++) {
00579 for(nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
00580 for(nzero=0; nzero<numSec/3; nzero++) {
00581 if(++counter < numMul) {
00582 nt = npos+nneg+nzero;
00583 if( (nt>0) && (nt<=numSec) ) {
00584 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c) ;
00585 protnorm[nt-1] += protmul[counter];
00586 }
00587 }
00588 }
00589 }
00590 }
00591
00592 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00593 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00594 counter = -1;
00595 for(npos=0; npos<numSec/3; npos++) {
00596 for(nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00597 for(nzero=0; nzero<numSec/3; nzero++) {
00598 if(++counter < numMul) {
00599 nt = npos+nneg+nzero;
00600 if( (nt>0) && (nt<=numSec) ) {
00601 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00602 neutnorm[nt-1] += neutmul[counter];
00603 }
00604 }
00605 }
00606 }
00607 }
00608
00609 for(i=0; i<numSec; i++) {
00610 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00611 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00612 }
00613 }
00614
00615
00616
00617 pv[0] = incidentParticle;
00618 pv[1] = targetParticle;
00619 vecLen = 2;
00620
00621 if (!inElastic || (availableEnergy <= PionPlus.getMass()))
00622 return;
00623
00624
00625
00626 npos = 0, nneg = 0, nzero = 0;
00627 G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
00628 G4int iplab = G4int( incidentTotalMomentum*5.);
00629 if ((iplab < 10) && (G4UniformRand() < cech[iplab]) ) {
00630 G4int ipl = std::min(19, G4int( incidentTotalMomentum*5.));
00631 G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
00632 0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
00633 if (G4UniformRand() < cnk0[ipl]) {
00634 if(targetCode == protonCode) {
00635 return;
00636 } else {
00637 pv[0] = KaonMinus;
00638 pv[1] = Proton;
00639 return;
00640 }
00641 }
00642
00643 G4double ran = G4UniformRand();
00644 if(targetCode == protonCode) {
00645
00646
00647 if( ran < 0.25 ) {
00648 ;
00649 } else if (ran < 0.50) {
00650 pv[0] = PionPlus;
00651 pv[1] = SigmaZero;
00652 } else if (ran < 0.75) {
00653 ;
00654 } else {
00655 pv[0] = PionPlus;
00656 pv[1] = Lambda;
00657 }
00658 } else {
00659
00660
00661 if( ran < 0.25 ) {
00662 pv[0] = PionMinus;
00663 pv[1] = SigmaPlus;
00664 } else if (ran < 0.50) {
00665 pv[0] = PionZero;
00666 pv[1] = SigmaZero;
00667 } else if (ran < 0.75) {
00668 pv[0] = PionPlus;
00669 pv[1] = SigmaMinus;
00670 } else {
00671 pv[0] = PionZero;
00672 pv[1] = Lambda;
00673 }
00674 }
00675 return;
00676
00677 } else {
00678
00679
00680 G4double aleab = std::log(availableEnergy);
00681 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00682 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00683
00684
00685
00686 G4double test, dum, anpn = 0.0;
00687
00688 for (nt=1; nt<=numSec; nt++) {
00689 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00690 dum = pi*nt/(2.0*n*n);
00691 if (std::fabs(dum) < 1.0) {
00692 if( test >= 1.0e-10 )anpn += dum*test;
00693 } else {
00694 anpn += dum*test;
00695 }
00696 }
00697
00698 G4double ran = G4UniformRand();
00699 G4double excs = 0.0;
00700 if (targetCode == protonCode) {
00701 counter = -1;
00702 for (npos=0; npos<numSec/3; npos++) {
00703 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
00704 for (nzero=0; nzero<numSec/3; nzero++) {
00705 if (++counter < numMul) {
00706 nt = npos+nneg+nzero;
00707 if( (nt>0) && (nt<=numSec) ) {
00708 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00709 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00710
00711 if (std::fabs(dum) < 1.0) {
00712 if( test >= 1.0e-10 )excs += dum*test;
00713 } else {
00714 excs += dum*test;
00715 }
00716
00717 if (ran < excs) goto outOfLoop;
00718 }
00719 }
00720 }
00721 }
00722 }
00723
00724 inElastic = false;
00725 return;
00726
00727 } else {
00728 counter = -1;
00729 for (npos=0; npos<numSec/3; npos++) {
00730 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00731 for (nzero=0; nzero<numSec/3; nzero++) {
00732 if (++counter < numMul) {
00733 nt = npos+nneg+nzero;
00734 if( (nt>=1) && (nt<=numSec) ) {
00735 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00736 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00737
00738 if (std::fabs(dum) < 1.0) {
00739 if( test >= 1.0e-10 )excs += dum*test;
00740 } else {
00741 excs += dum*test;
00742 }
00743
00744 if (ran < excs) goto outOfLoop;
00745 }
00746 }
00747 }
00748 }
00749 }
00750
00751 inElastic = false;
00752 return;
00753 }
00754 }
00755 outOfLoop:
00756
00757 if( targetCode == protonCode)
00758 {
00759 if( npos == nneg)
00760 {
00761 }
00762 else if (npos == (1+nneg))
00763 {
00764 if( G4UniformRand() < 0.5)
00765 {
00766 pv[0] = KaonMinus;
00767 }
00768 else
00769 {
00770 pv[1] = Neutron;
00771 }
00772 }
00773 else
00774 {
00775 pv[0] = KaonMinus;
00776 pv[1] = Neutron;
00777 }
00778 }
00779 else
00780 {
00781 if( npos == nneg)
00782 {
00783 if( G4UniformRand() < 0.75)
00784 {
00785 }
00786 else
00787 {
00788 pv[0] = KaonMinus;
00789 pv[1] = Proton;
00790 }
00791 }
00792 else if ( npos == (1+nneg))
00793 {
00794 pv[0] = KaonMinus;
00795 }
00796 else
00797 {
00798 pv[1] = Proton;
00799 }
00800 }
00801
00802
00803 if( G4UniformRand() < 0.5 )
00804 {
00805 if( ( (pv[0].getCode() == kaonMinusCode)
00806 && (pv[1].getCode() == neutronCode) )
00807 || ( (pv[0].getCode() == kaonZeroCode)
00808 && (pv[1].getCode() == protonCode) )
00809 || ( (pv[0].getCode() == antiKaonZeroCode)
00810 && (pv[1].getCode() == protonCode) ) )
00811 {
00812 G4double ran = G4UniformRand();
00813 if( pv[1].getCode() == protonCode)
00814 {
00815 if(ran < 0.68)
00816 {
00817 pv[0] = PionPlus;
00818 pv[1] = Lambda;
00819 }
00820 else if (ran < 0.84)
00821 {
00822 pv[0] = PionZero;
00823 pv[1] = SigmaPlus;
00824 }
00825 else
00826 {
00827 pv[0] = PionPlus;
00828 pv[1] = SigmaZero;
00829 }
00830 }
00831 else
00832 {
00833 if(ran < 0.68)
00834 {
00835 pv[0] = PionMinus;
00836 pv[1] = Lambda;
00837 }
00838 else if (ran < 0.84)
00839 {
00840 pv[0] = PionMinus;
00841 pv[1] = SigmaZero;
00842 }
00843 else
00844 {
00845 pv[0] = PionZero;
00846 pv[1] = SigmaMinus;
00847 }
00848 }
00849 }
00850 else
00851 {
00852 G4double ran = G4UniformRand();
00853 if (ran < 0.67)
00854 {
00855 pv[0] = PionZero;
00856 pv[1] = Lambda;
00857 }
00858 else if (ran < 0.78)
00859 {
00860 pv[0] = PionMinus;
00861 pv[1] = SigmaPlus;
00862 }
00863 else if (ran < 0.89)
00864 {
00865 pv[0] = PionZero;
00866 pv[1] = SigmaZero;
00867 }
00868 else
00869 {
00870 pv[0] = PionPlus;
00871 pv[1] = SigmaMinus;
00872 }
00873 }
00874 }
00875
00876 nt = npos + nneg + nzero;
00877 while ( nt > 0) {
00878 G4double ran = G4UniformRand();
00879 if ( ran < (G4double)npos/nt) {
00880 if( npos > 0 ) {
00881 pv[vecLen++] = PionPlus;
00882 npos--;
00883 }
00884 } else if (ran < (G4double)(npos+nneg)/nt) {
00885 if( nneg > 0 ) {
00886 pv[vecLen++] = PionMinus;
00887 nneg--;
00888 }
00889 } else {
00890 if( nzero > 0 ) {
00891 pv[vecLen++] = PionZero;
00892 nzero--;
00893 }
00894 }
00895 nt = npos + nneg + nzero;
00896 }
00897
00898 if (verboseLevel > 1) {
00899 G4cout << "Particles produced: " ;
00900 G4cout << pv[0].getName() << " " ;
00901 G4cout << pv[1].getName() << " " ;
00902 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
00903 G4cout << G4endl;
00904 }
00905
00906 return;
00907 }