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 "G4HEOmegaMinusInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041
00042 void G4HEOmegaMinusInelastic::ModelDescription(std::ostream& outFile) const
00043 {
00044 outFile << "G4HEOmegaMinusInelastic is one of the High Energy\n"
00045 << "Parameterized (HEP) models used to implement inelastic\n"
00046 << "Omega- 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 Omega- with initial energies\n"
00052 << "above 20 GeV.\n";
00053 }
00054
00055
00056 G4HadFinalState*
00057 G4HEOmegaMinusInelastic::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 << "GHEOmegaMinusInelastic: incident energy < 1 GeV" << G4endl;
00080
00081 if (verboseLevel > 1) {
00082 G4cout << "G4HEOmegaMinusInelastic::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 if (G4UniformRand() < atomicNumber/atomicWeight) {
00117 targetParticle.setDefinition("Proton");
00118 } else {
00119 targetParticle.setDefinition("Neutron");
00120 }
00121
00122 G4double targetMass = targetParticle.getMass();
00123 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00124 + targetMass*targetMass
00125 + 2.0*targetMass*incidentTotalEnergy);
00126 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00127
00128 G4bool inElastic = true;
00129 vecLength = 0;
00130
00131 if (verboseLevel > 1)
00132 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00133 << incidentCode << G4endl;
00134
00135 G4bool successful = false;
00136
00137 FirstIntInCasOmegaMinus(inElastic, availableEnergy, pv, vecLength,
00138 incidentParticle, targetParticle, atomicWeight);
00139
00140 if (verboseLevel > 1)
00141 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00142
00143 if ((vecLength > 0) && (availableEnergy > 1.))
00144 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00145 pv, vecLength,
00146 incidentParticle, targetParticle);
00147
00148 HighEnergyCascading(successful, pv, vecLength,
00149 excitationEnergyGNP, excitationEnergyDTA,
00150 incidentParticle, targetParticle,
00151 atomicWeight, atomicNumber);
00152 if (!successful)
00153 HighEnergyClusterProduction(successful, pv, vecLength,
00154 excitationEnergyGNP, excitationEnergyDTA,
00155 incidentParticle, targetParticle,
00156 atomicWeight, atomicNumber);
00157 if (!successful)
00158 MediumEnergyCascading(successful, pv, vecLength,
00159 excitationEnergyGNP, excitationEnergyDTA,
00160 incidentParticle, targetParticle,
00161 atomicWeight, atomicNumber);
00162
00163 if (!successful)
00164 MediumEnergyClusterProduction(successful, pv, vecLength,
00165 excitationEnergyGNP, excitationEnergyDTA,
00166 incidentParticle, targetParticle,
00167 atomicWeight, atomicNumber);
00168 if (!successful)
00169 QuasiElasticScattering(successful, pv, vecLength,
00170 excitationEnergyGNP, excitationEnergyDTA,
00171 incidentParticle, targetParticle,
00172 atomicWeight, atomicNumber);
00173 if (!successful)
00174 ElasticScattering(successful, pv, vecLength,
00175 incidentParticle,
00176 atomicWeight, atomicNumber);
00177
00178 if (!successful)
00179 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00180 << G4endl;
00181
00182 FillParticleChange(pv, vecLength);
00183 delete [] pv;
00184 theParticleChange.SetStatusChange(stopAndKill);
00185 return &theParticleChange;
00186 }
00187
00188
00189 void
00190 G4HEOmegaMinusInelastic::FirstIntInCasOmegaMinus(G4bool& inElastic,
00191 const G4double availableEnergy,
00192 G4HEVector pv[],
00193 G4int& vecLen,
00194 const G4HEVector& incidentParticle,
00195 const G4HEVector& targetParticle,
00196 const G4double atomicWeight)
00197
00198
00199
00200
00201
00202
00203
00204
00205 {
00206 static const G4double expxu = 82.;
00207 static const G4double expxl = -expxu;
00208
00209 static const G4double protb = 0.7;
00210 static const G4double neutb = 0.7;
00211 static const G4double c = 1.25;
00212
00213 static const G4int numMul = 1200;
00214 static const G4int numSec = 60;
00215
00216 G4int protonCode = Proton.getCode();
00217
00218 G4int targetCode = targetParticle.getCode();
00219 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00220
00221 static G4bool first = true;
00222 static G4double protmul[numMul], protnorm[numSec];
00223 static G4double neutmul[numMul], neutnorm[numSec];
00224
00225
00226
00227
00228 G4int i, counter, nt, npos, nneg, nzero;
00229
00230 if (first) {
00231
00232 first = false;
00233 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
00234 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00235 counter = -1;
00236 for( npos=0; npos<(numSec/3); npos++ )
00237 {
00238 for( nneg=std::max(0,npos-1); nneg<=npos; nneg++ )
00239 {
00240 for( nzero=0; nzero<numSec/3; nzero++ )
00241 {
00242 if( ++counter < numMul )
00243 {
00244 nt = npos+nneg+nzero;
00245 if( (nt>0) && (nt<=numSec) )
00246 {
00247 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00248 protnorm[nt-1] += protmul[counter];
00249 }
00250 }
00251 }
00252 }
00253 }
00254 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00255 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00256 counter = -1;
00257 for( npos=0; npos<numSec/3; npos++ )
00258 {
00259 for( nneg=npos; nneg<=(npos+1); nneg++ )
00260 {
00261 for( nzero=0; nzero<numSec/3; nzero++ )
00262 {
00263 if( ++counter < numMul )
00264 {
00265 nt = npos+nneg+nzero;
00266 if( (nt>0) && (nt<=numSec) )
00267 {
00268 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00269 neutnorm[nt-1] += neutmul[counter];
00270 }
00271 }
00272 }
00273 }
00274 }
00275 for( i=0; i<numSec; i++ )
00276 {
00277 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00278 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00279 }
00280 }
00281
00282
00283
00284
00285 pv[0] = incidentParticle;
00286 pv[1] = targetParticle;
00287 vecLen = 2;
00288
00289 if( !inElastic )
00290 {
00291 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00292 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
00293 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
00294 {
00295 G4double ran = G4UniformRand();
00296 if( targetCode == protonCode)
00297 {
00298 if (ran < 0.2)
00299 {
00300 pv[0] = XiZero;
00301 pv[1] = SigmaZero;
00302 }
00303 else if (ran < 0.4)
00304 {
00305 pv[0] = SigmaZero;
00306 pv[1] = XiZero;
00307 }
00308 else if (ran < 0.6)
00309 {
00310 pv[0] = XiZero;
00311 pv[1] = Lambda;
00312 }
00313 else if (ran < 0.8)
00314 {
00315 pv[0] = Lambda;
00316 pv[1] = XiZero;
00317 }
00318 else
00319 {
00320 pv[0] = Proton;
00321 pv[1] = OmegaMinus;
00322 }
00323 }
00324 else
00325 {
00326 if (ran < 0.2)
00327 {
00328 pv[0] = Neutron;
00329 pv[1] = OmegaMinus;
00330 }
00331 else if (ran < 0.4)
00332 {
00333 pv[0] = XiZero;
00334 pv[1] = SigmaMinus;
00335 }
00336 else if (ran < 0.6)
00337 {
00338 pv[0] = SigmaMinus;
00339 pv[1] = XiZero;
00340 }
00341 else if (ran < 0.8)
00342 {
00343 pv[0] = XiMinus;
00344 pv[1] = Lambda;
00345 }
00346 else
00347 {
00348 pv[0] = Lambda;
00349 pv[1] = XiMinus;
00350 }
00351 }
00352 }
00353 return;
00354 }
00355 else if (availableEnergy <= PionPlus.getMass())
00356 return;
00357
00358
00359 npos = 0; nneg = 0; nzero = 0;
00360
00361
00362 G4double aleab = std::log(availableEnergy);
00363 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00364 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00365
00366
00367
00368 G4double test, dum, anpn = 0.0;
00369
00370 for (nt=1; nt<=numSec; nt++) {
00371 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00372 dum = pi*nt/(2.0*n*n);
00373 if (std::fabs(dum) < 1.0) {
00374 if( test >= 1.0e-10 )anpn += dum*test;
00375 } else {
00376 anpn += dum*test;
00377 }
00378 }
00379
00380 G4double ran = G4UniformRand();
00381 G4double excs = 0.0;
00382 if (targetCode == protonCode) {
00383 counter = -1;
00384 for (npos = 0; npos < numSec/3; npos++) {
00385 for( nneg=std::max(0,npos-1); nneg<=npos; nneg++ )
00386 {
00387 for( nzero=0; nzero<numSec/3; nzero++ )
00388 {
00389 if( ++counter < numMul )
00390 {
00391 nt = npos+nneg+nzero;
00392 if ( (nt>0) && (nt<=numSec) ) {
00393 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00394 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00395 if (std::fabs(dum) < 1.0) {
00396 if( test >= 1.0e-10 )excs += dum*test;
00397 } else {
00398 excs += dum*test;
00399 }
00400 if (ran < excs) goto outOfLoop;
00401 }
00402 }
00403 }
00404 }
00405 }
00406
00407
00408 inElastic = false;
00409 return;
00410 }
00411 else
00412 {
00413 counter = -1;
00414 for( npos=0; npos<numSec/3; npos++ )
00415 {
00416 for( nneg=npos; nneg<=(npos+1); nneg++ )
00417 {
00418 for( nzero=0; nzero<numSec/3; nzero++ )
00419 {
00420 if( ++counter < numMul )
00421 {
00422 nt = npos+nneg+nzero;
00423 if ( (nt>=1) && (nt<=numSec) ) {
00424 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00425 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00426 if (std::fabs(dum) < 1.0) {
00427 if( test >= 1.0e-10 )excs += dum*test;
00428 } else {
00429 excs += dum*test;
00430 }
00431 if (ran < excs) goto outOfLoop;
00432 }
00433 }
00434 }
00435 }
00436 }
00437
00438 inElastic = false;
00439 return;
00440 }
00441
00442 outOfLoop:
00443
00444
00445
00446
00447 ran = G4UniformRand();
00448 if (targetCode == protonCode) {
00449 if( npos == nneg)
00450 {
00451 }
00452 else
00453 {
00454 pv[1] = Neutron;
00455 }
00456 } else {
00457 if (npos == nneg)
00458 {
00459 }
00460 else
00461 {
00462 pv[1] = Proton;
00463 }
00464 }
00465
00466 nt = npos + nneg + nzero;
00467 while (nt > 0) {
00468 G4double rnd = G4UniformRand();
00469 if (rnd < (G4double)npos/nt) {
00470 if (npos > 0) {
00471 pv[vecLen++] = PionPlus;
00472 npos--;
00473 }
00474 } else if (rnd < (G4double)(npos+nneg)/nt) {
00475 if (nneg > 0) {
00476 pv[vecLen++] = PionMinus;
00477 nneg--;
00478 }
00479 } else {
00480 if (nzero > 0) {
00481 pv[vecLen++] = PionZero;
00482 nzero--;
00483 }
00484 }
00485 nt = npos + nneg + nzero;
00486 }
00487
00488 if (verboseLevel > 1) {
00489 G4cout << "Particles produced: " ;
00490 G4cout << pv[0].getName() << " " ;
00491 G4cout << pv[1].getName() << " " ;
00492 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
00493 G4cout << G4endl;
00494 }
00495 return;
00496 }
00497