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 #include "G4HESigmaPlusInelastic.hh"
00037 #include "globals.hh"
00038 #include "G4ios.hh"
00039 #include "G4PhysicalConstants.hh"
00040 #include "G4SystemOfUnits.hh"
00041
00042
00043 G4HESigmaPlusInelastic::G4HESigmaPlusInelastic(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 G4HESigmaPlusInelastic is being deprecated and will\n"
00052 << "disappear in Geant4 version 10.0" << G4endl;
00053 }
00054
00055 void G4HESigmaPlusInelastic::ModelDescription(std::ostream& outFile) const
00056 {
00057 outFile << "G4HESigmaPlusInelastic is one of the High Energy\n"
00058 << "Parameterized (HEP) models used to implement inelastic\n"
00059 << "Sigma+ scattering from nuclei. It is a re-engineered\n"
00060 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00061 << "initial collision products into backward- and forward-going\n"
00062 << "clusters which are then decayed into final state hadrons.\n"
00063 << "The model does not conserve energy on an event-by-event\n"
00064 << "basis. It may be applied to Sigma+ with initial energies\n"
00065 << "above 20 GeV.\n";
00066 }
00067
00068
00069 G4HadFinalState*
00070 G4HESigmaPlusInelastic::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 << "GHESigmaPlusInelastic: incident energy < 1 GeV" << G4endl;
00093
00094 if (verboseLevel > 1) {
00095 G4cout << "G4HESigmaPlusInelastic::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 FirstIntInCasSigmaPlus(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
00195 FillParticleChange(pv, vecLength);
00196 delete [] pv;
00197 theParticleChange.SetStatusChange(stopAndKill);
00198 return &theParticleChange;
00199 }
00200
00201
00202 void
00203 G4HESigmaPlusInelastic::FirstIntInCasSigmaPlus(G4bool& inElastic,
00204 const G4double availableEnergy,
00205 G4HEVector pv[],
00206 G4int& vecLen,
00207 const G4HEVector& incidentParticle,
00208 const G4HEVector& targetParticle,
00209 const G4double atomicWeight)
00210
00211
00212
00213
00214
00215
00216
00217
00218 {
00219 static const G4double expxu = 82.;
00220 static const G4double expxl = -expxu;
00221
00222 static const G4double protb = 0.7;
00223 static const G4double neutb = 0.7;
00224 static const G4double c = 1.25;
00225
00226 static const G4int numMul = 1200;
00227 static const G4int numSec = 60;
00228
00229 G4int protonCode = Proton.getCode();
00230
00231 G4int targetCode = targetParticle.getCode();
00232 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00233
00234 static G4bool first = true;
00235 static G4double protmul[numMul], protnorm[numSec];
00236 static G4double neutmul[numMul], neutnorm[numSec];
00237
00238
00239
00240
00241 G4int i, counter, nt, npos, nneg, nzero;
00242
00243 if (first) {
00244
00245 first = false;
00246 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
00247 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00248 counter = -1;
00249 for( npos=0; npos<(numSec/3); npos++ )
00250 {
00251 for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ )
00252 {
00253 for( nzero=0; nzero<numSec/3; nzero++ )
00254 {
00255 if( ++counter < numMul )
00256 {
00257 nt = npos+nneg+nzero;
00258 if( (nt>0) && (nt<=numSec) )
00259 {
00260 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00261 protnorm[nt-1] += protmul[counter];
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 {
00272 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
00273 {
00274 for( nzero=0; nzero<numSec/3; nzero++ )
00275 {
00276 if( ++counter < numMul )
00277 {
00278 nt = npos+nneg+nzero;
00279 if( (nt>0) && (nt<=numSec) )
00280 {
00281 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00282 neutnorm[nt-1] += neutmul[counter];
00283 }
00284 }
00285 }
00286 }
00287 }
00288 for( i=0; i<numSec; i++ )
00289 {
00290 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00291 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00292 }
00293 }
00294
00295
00296
00297
00298 pv[0] = incidentParticle;
00299 pv[1] = targetParticle;
00300 vecLen = 2;
00301
00302 if( !inElastic )
00303 {
00304 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00305 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
00306 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
00307 {
00308 G4double ran = G4UniformRand();
00309 if( targetCode == protonCode)
00310 {
00311 pv[0] = Proton;
00312 pv[1] = SigmaPlus;
00313 }
00314 else
00315 {
00316 if(ran < 0.2)
00317 {
00318 pv[0] = SigmaZero;
00319 pv[1] = Proton;
00320 }
00321 else if(ran < 0.4)
00322 {
00323 pv[0] = Lambda;
00324 pv[1] = Proton;
00325 }
00326 else if(ran < 0.6)
00327 {
00328 pv[0] = Neutron;
00329 pv[1] = SigmaPlus;
00330 }
00331 else if(ran < 0.8)
00332 {
00333 pv[0] = Proton;
00334 pv[1] = SigmaZero;
00335 }
00336 else
00337 {
00338 pv[0] = Proton;
00339 pv[1] = Lambda;
00340 }
00341 }
00342 }
00343 return;
00344 }
00345 else if (availableEnergy <= PionPlus.getMass())
00346 return;
00347
00348
00349
00350 npos = 0; nneg = 0; nzero = 0;
00351
00352
00353
00354 G4double aleab = std::log(availableEnergy);
00355 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00356 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00357
00358
00359
00360 G4double test, dum, anpn = 0.0;
00361
00362 for (nt=1; nt<=numSec; nt++) {
00363 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00364 dum = pi*nt/(2.0*n*n);
00365 if (std::fabs(dum) < 1.0) {
00366 if (test >= 1.0e-10) anpn += dum*test;
00367 } else {
00368 anpn += dum*test;
00369 }
00370 }
00371
00372 G4double ran = G4UniformRand();
00373 G4double excs = 0.0;
00374 if( targetCode == protonCode )
00375 {
00376 counter = -1;
00377 for (npos=0; npos<numSec/3; npos++) {
00378 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
00379 for (nzero=0; nzero<numSec/3; nzero++) {
00380 if (++counter < numMul) {
00381 nt = npos+nneg+nzero;
00382 if ( (nt>0) && (nt<=numSec) ) {
00383 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00384 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00385 if (std::fabs(dum) < 1.0) {
00386 if (test >= 1.0e-10) excs += dum*test;
00387 } else {
00388 excs += dum*test;
00389 }
00390 if (ran < excs) goto outOfLoop;
00391 }
00392 }
00393 }
00394 }
00395 }
00396
00397
00398 inElastic = false;
00399 return;
00400 }
00401 else
00402 {
00403 counter = -1;
00404 for (npos=0; npos<numSec/3; npos++) {
00405 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00406 for (nzero=0; nzero<numSec/3; nzero++) {
00407 if (++counter < numMul) {
00408 nt = npos+nneg+nzero;
00409 if ( (nt>=1) && (nt<=numSec) ) {
00410 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00411 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00412 if (std::fabs(dum) < 1.0) {
00413 if( test >= 1.0e-10 )excs += dum*test;
00414 } else {
00415 excs += dum*test;
00416 }
00417 if (ran < excs) goto outOfLoop;
00418 }
00419 }
00420 }
00421 }
00422 }
00423
00424 inElastic = false;
00425 return;
00426 }
00427
00428 outOfLoop:
00429
00430 ran = G4UniformRand();
00431 if (targetCode == protonCode) {
00432 if( npos == nneg)
00433 {
00434 }
00435 else if (npos == (nneg+1))
00436 {
00437 if( ran < 0.25)
00438 {
00439 pv[0] = SigmaZero;
00440 }
00441 else if(ran < 0.5)
00442 {
00443 pv[0] = Lambda;
00444 }
00445 else
00446 {
00447 pv[1] = Neutron;
00448 }
00449 }
00450 else
00451 {
00452 if(ran < 0.5)
00453 {
00454 pv[0] = SigmaZero;
00455 pv[1] = Neutron;
00456 }
00457 else
00458 {
00459 pv[0] = Lambda;
00460 pv[1] = Neutron;
00461 }
00462 }
00463 } else {
00464 if (npos == nneg)
00465 {
00466 if (ran < 0.5)
00467 {
00468 }
00469 else if (ran < 0.75)
00470 {
00471 pv[0] = SigmaZero;
00472 pv[1] = Proton;
00473 }
00474 else
00475 {
00476 pv[0] = Lambda;
00477 pv[1] = Proton;
00478 }
00479 }
00480 else if (npos == (nneg-1))
00481 {
00482 pv[1] = Proton;
00483 }
00484 else
00485 {
00486 if (ran < 0.5)
00487 {
00488 pv[0] = SigmaZero;
00489 }
00490 else
00491 {
00492 pv[0] = Lambda;
00493 }
00494 }
00495 }
00496
00497 nt = npos + nneg + nzero;
00498 while (nt > 0) {
00499 G4double rnd = G4UniformRand();
00500 if (rnd < (G4double)npos/nt) {
00501 if (npos > 0) {
00502 pv[vecLen++] = PionPlus;
00503 npos--;
00504 }
00505 } else if (rnd < (G4double)(npos+nneg)/nt) {
00506 if (nneg > 0) {
00507 pv[vecLen++] = PionMinus;
00508 nneg--;
00509 }
00510 } else {
00511 if (nzero > 0) {
00512 pv[vecLen++] = PionZero;
00513 nzero--;
00514 }
00515 }
00516 nt = npos + nneg + nzero;
00517 }
00518
00519 if (verboseLevel > 1) {
00520 G4cout << "Particles produced: " ;
00521 G4cout << pv[0].getName() << " " ;
00522 G4cout << pv[1].getName() << " " ;
00523 for (i = 2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
00524 G4cout << G4endl;
00525 }
00526
00527 return;
00528 }
00529