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