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 "G4HEKaonPlusInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041
00042 void G4HEKaonPlusInelastic::ModelDescription(std::ostream& outFile) const
00043 {
00044 outFile << "G4HEKaonPlusInelastic is one of the High Energy\n"
00045 << "Parameterized (HEP) models used to implement inelastic\n"
00046 << "K+ scattering from nuclei. It is a re-engineered\n"
00047 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00048 << "initial collision products into backward- and forward-going\n"
00049 << "clusters which are then decayed into final state hadrons.\n"
00050 << "The model does not conserve energy on an event-by-event\n"
00051 << "basis. It may be applied to K+ with initial energies\n"
00052 << "above 20 GeV.\n";
00053 }
00054
00055
00056 G4HadFinalState*
00057 G4HEKaonPlusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00058 G4Nucleus& targetNucleus)
00059 {
00060 G4HEVector* pv = new G4HEVector[MAXPART];
00061 const G4HadProjectile* aParticle = &aTrack;
00062 const G4double A = targetNucleus.GetA_asInt();
00063 const G4double Z = targetNucleus.GetZ_asInt();
00064 G4HEVector incidentParticle(aParticle);
00065
00066 G4double atomicNumber = Z;
00067 G4double atomicWeight = A;
00068
00069 G4int incidentCode = incidentParticle.getCode();
00070 G4double incidentMass = incidentParticle.getMass();
00071 G4double incidentTotalEnergy = incidentParticle.getEnergy();
00072
00073
00074
00075
00076 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00077
00078 if (incidentKineticEnergy < 1.)
00079 G4cout << "GHEKaonPlusInelastic: incident energy < 1 GeV" << G4endl;
00080
00081 if (verboseLevel > 1) {
00082 G4cout << "G4HEKaonPlusInelastic::ApplyYourself" << G4endl;
00083 G4cout << "incident particle " << incidentParticle.getName()
00084 << "mass " << incidentMass
00085 << "kinetic energy " << incidentKineticEnergy
00086 << G4endl;
00087 G4cout << "target material with (A,Z) = ("
00088 << atomicWeight << "," << atomicNumber << ")" << G4endl;
00089 }
00090
00091 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
00092 atomicWeight, atomicNumber);
00093 if (verboseLevel > 1)
00094 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00095
00096 incidentKineticEnergy -= inelasticity;
00097
00098 G4double excitationEnergyGNP = 0.;
00099 G4double excitationEnergyDTA = 0.;
00100
00101 G4double excitation = NuclearExcitation(incidentKineticEnergy,
00102 atomicWeight, atomicNumber,
00103 excitationEnergyGNP,
00104 excitationEnergyDTA);
00105 if (verboseLevel > 1)
00106 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
00107 << excitationEnergyDTA << G4endl;
00108
00109 incidentKineticEnergy -= excitation;
00110 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00111
00112
00113
00114
00115 G4HEVector targetParticle;
00116
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 vecLength = 0;
00131
00132 if (verboseLevel > 1)
00133 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00134 << incidentCode << G4endl;
00135
00136 G4bool successful = false;
00137
00138 FirstIntInCasKaonPlus(inElastic, availableEnergy, pv, vecLength,
00139 incidentParticle, targetParticle, atomicWeight);
00140
00141 if (verboseLevel > 1)
00142 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
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 if (!successful)
00175 ElasticScattering(successful, pv, vecLength,
00176 incidentParticle,
00177 atomicWeight, atomicNumber);
00178
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 G4HEKaonPlusInelastic::FirstIntInCasKaonPlus(G4bool& inElastic,
00192 const G4double availableEnergy,
00193 G4HEVector pv[],
00194 G4int& vecLen,
00195 const G4HEVector& incidentParticle,
00196 const G4HEVector& targetParticle,
00197 const G4double atomicWeight)
00198
00199
00200
00201
00202
00203
00204
00205
00206 {
00207 static const G4double expxu = 82.;
00208 static const G4double expxl = -expxu;
00209
00210 static const G4double protb = 0.7;
00211 static const G4double neutb = 0.7;
00212 static const G4double c = 1.25;
00213
00214 static const G4int numMul = 1200;
00215 static const G4int numSec = 60;
00216
00217 G4int neutronCode = Neutron.getCode();
00218 G4int protonCode = Proton.getCode();
00219
00220 G4int targetCode = targetParticle.getCode();
00221 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00222
00223 static G4bool first = true;
00224 static G4double protmul[numMul], protnorm[numSec];
00225 static G4double neutmul[numMul], neutnorm[numSec];
00226
00227
00228
00229
00230 G4int i, counter, nt, npos, nneg, nzero;
00231
00232 if( first )
00233 {
00234 first = false;
00235 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
00236 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00237 counter = -1;
00238 for( npos=0; npos<(numSec/3); npos++ )
00239 {
00240 for( nneg=Imax(0,npos-2); nneg<=npos; nneg++ )
00241 {
00242 for( nzero=0; nzero<numSec/3; nzero++ )
00243 {
00244 if( ++counter < numMul )
00245 {
00246 nt = npos+nneg+nzero;
00247 if( (nt>0) && (nt<=numSec) )
00248 {
00249 protmul[counter] =
00250 pmltpc(npos,nneg,nzero,nt,protb,c) ;
00251 protnorm[nt-1] += protmul[counter];
00252 }
00253 }
00254 }
00255 }
00256 }
00257 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00258 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00259 counter = -1;
00260 for( npos=0; npos<numSec/3; npos++ )
00261 {
00262 for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ )
00263 {
00264 for( nzero=0; nzero<numSec/3; nzero++ )
00265 {
00266 if( ++counter < numMul )
00267 {
00268 nt = npos+nneg+nzero;
00269 if( (nt>0) && (nt<=numSec) )
00270 {
00271 neutmul[counter] =
00272 pmltpc(npos,nneg,nzero,nt,neutb,c);
00273 neutnorm[nt-1] += neutmul[counter];
00274 }
00275 }
00276 }
00277 }
00278 }
00279 for( i=0; i<numSec; i++ )
00280 {
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
00287
00288
00289 pv[0] = incidentParticle;
00290 pv[1] = targetParticle;
00291 vecLen = 2;
00292
00293 if( !inElastic )
00294 {
00295 if( targetCode == neutronCode )
00296 {
00297 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07};
00298 G4int iplab = G4int( Amin( 9.0, incidentTotalMomentum*5. ) );
00299 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
00300 {
00301 pv[0] = KaonZero;
00302 pv[1] = Proton;
00303 }
00304 }
00305 return;
00306 }
00307 else if (availableEnergy <= PionPlus.getMass())
00308 return;
00309
00310
00311
00312 npos = 0, nneg = 0, nzero = 0;
00313 G4double eab = availableEnergy;
00314 G4int ieab = G4int( eab*5.0 );
00315
00316 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
00317 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
00318 {
00319
00320
00321 G4double w0, wp, wm, wt, ran;
00322 if( targetCode == protonCode )
00323 {
00324 w0 = - sqr(1.+protb)/(2.*c*c);
00325 wp = w0 = std::exp(w0);
00326 wp *= 2.;
00327 if( G4UniformRand() < w0/(w0+wp) )
00328 { npos = 0; nneg = 0; nzero = 1; }
00329 else
00330 { npos = 1; nneg = 0; nzero = 0; }
00331 }
00332 else
00333 {
00334 w0 = -sqr(1.+neutb)/(2.*c*c);
00335 wp = w0 = std::exp(w0);
00336 wm = -sqr(-1.+neutb)/(2.*c*c);
00337 wm = std::exp(wm);
00338 wt = w0+wp+wm;
00339 wp = w0+wp;
00340 ran = G4UniformRand();
00341 if( ran < w0/wt)
00342 { npos = 0; nneg = 0; nzero = 1; }
00343 else if( ran < wp/wt)
00344 { npos = 1; nneg = 0; nzero = 0; }
00345 else
00346 { npos = 0; nneg = 1; nzero = 0; }
00347 }
00348 }
00349 else
00350 {
00351
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-2); nneg<=npos; 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 if (ran < excs) goto outOfLoop;
00392 }
00393 }
00394 }
00395 }
00396 }
00397
00398
00399 inElastic = false;
00400 return;
00401 }
00402 else
00403 {
00404 counter = -1;
00405 for( npos=0; npos<numSec/3; npos++ )
00406 {
00407 for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ )
00408 {
00409 for (nzero=0; nzero<numSec/3; nzero++) {
00410 if (++counter < numMul) {
00411 nt = npos+nneg+nzero;
00412 if ( (nt>=1) && (nt<=numSec) ) {
00413 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00414 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00415 if (std::fabs(dum) < 1.0) {
00416 if( test >= 1.0e-10 )excs += dum*test;
00417 } else {
00418 excs += dum*test;
00419 }
00420 if (ran < excs) goto outOfLoop;
00421 }
00422 }
00423 }
00424 }
00425 }
00426
00427 inElastic = false;
00428 return;
00429 }
00430 }
00431 outOfLoop:
00432
00433 if( targetCode == protonCode)
00434 {
00435 if( npos == nneg)
00436 {
00437 }
00438 else if (npos == (1+nneg))
00439 {
00440 if( G4UniformRand() < 0.5)
00441 {
00442 pv[1] = Neutron;
00443 }
00444 else
00445 {
00446 pv[0] = KaonZero;
00447 }
00448 }
00449 else
00450 {
00451 pv[0] = KaonZero;
00452 pv[1] = Neutron;
00453 }
00454 }
00455 else
00456 {
00457 if( npos == nneg)
00458 {
00459 if( G4UniformRand() < 0.25)
00460 {
00461 pv[0] = KaonZero;
00462 pv[1] = Proton;
00463 }
00464 else
00465 {
00466 }
00467 }
00468 else if ( npos == (1+nneg))
00469 {
00470 pv[0] = KaonZero;
00471 }
00472 else
00473 {
00474 pv[1] = Proton;
00475 }
00476 }
00477
00478
00479 nt = npos + nneg + nzero;
00480 while ( nt > 0)
00481 {
00482 G4double ran = G4UniformRand();
00483 if ( ran < (G4double)npos/nt)
00484 {
00485 if( npos > 0 )
00486 { pv[vecLen++] = PionPlus;
00487 npos--;
00488 }
00489 }
00490 else if ( ran < (G4double)(npos+nneg)/nt)
00491 {
00492 if( nneg > 0 )
00493 {
00494 pv[vecLen++] = PionMinus;
00495 nneg--;
00496 }
00497 }
00498 else
00499 {
00500 if( nzero > 0 )
00501 {
00502 pv[vecLen++] = PionZero;
00503 nzero--;
00504 }
00505 }
00506 nt = npos + nneg + nzero;
00507 }
00508 if (verboseLevel > 1)
00509 {
00510 G4cout << "Particles produced: " ;
00511 G4cout << pv[0].getName() << " " ;
00512 G4cout << pv[1].getName() << " " ;
00513 for (i=2; i < vecLen; i++)
00514 {
00515 G4cout << pv[i].getName() << " " ;
00516 }
00517 G4cout << G4endl;
00518 }
00519 return;
00520 }
00521
00522
00523
00524
00525
00526
00527
00528
00529