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 "G4HEXiMinusInelastic.hh"
00037 #include "globals.hh"
00038 #include "G4ios.hh"
00039 #include "G4PhysicalConstants.hh"
00040
00041 void G4HEXiMinusInelastic::ModelDescription(std::ostream& outFile) const
00042 {
00043 outFile << "G4HEXiMinusInelastic is one of the High Energy\n"
00044 << "Parameterized (HEP) models used to implement inelastic\n"
00045 << "Xi- scattering from nuclei. It is a re-engineered\n"
00046 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00047 << "initial collision products into backward- and forward-going\n"
00048 << "clusters which are then decayed into final state hadrons.\n"
00049 << "The model does not conserve energy on an event-by-event\n"
00050 << "basis. It may be applied to Xi- with initial energies\n"
00051 << "above 20 GeV.\n";
00052 }
00053
00054
00055 G4HadFinalState*
00056 G4HEXiMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00057 G4Nucleus& targetNucleus)
00058 {
00059 G4HEVector* pv = new G4HEVector[MAXPART];
00060 const G4HadProjectile* aParticle = &aTrack;
00061 const G4double A = targetNucleus.GetA_asInt();
00062 const G4double Z = targetNucleus.GetZ_asInt();
00063 G4HEVector incidentParticle(aParticle);
00064
00065 G4double atomicNumber = Z;
00066 G4double atomicWeight = A;
00067
00068 G4int incidentCode = incidentParticle.getCode();
00069 G4double incidentMass = incidentParticle.getMass();
00070 G4double incidentTotalEnergy = incidentParticle.getEnergy();
00071
00072
00073
00074
00075 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00076
00077 if (incidentKineticEnergy < 1.)
00078 G4cout << "GHEXiMinusInelastic: incident energy < 1 GeV" << G4endl;
00079
00080 if (verboseLevel > 1) {
00081 G4cout << "G4HEXiMinusInelastic::ApplyYourself" << G4endl;
00082 G4cout << "incident particle " << incidentParticle.getName()
00083 << "mass " << incidentMass
00084 << "kinetic energy " << incidentKineticEnergy
00085 << G4endl;
00086 G4cout << "target material with (A,Z) = ("
00087 << atomicWeight << "," << atomicNumber << ")" << G4endl;
00088 }
00089
00090 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
00091 atomicWeight, atomicNumber);
00092 if (verboseLevel > 1)
00093 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00094
00095 incidentKineticEnergy -= inelasticity;
00096
00097 G4double excitationEnergyGNP = 0.;
00098 G4double excitationEnergyDTA = 0.;
00099
00100 G4double excitation = NuclearExcitation(incidentKineticEnergy,
00101 atomicWeight, atomicNumber,
00102 excitationEnergyGNP,
00103 excitationEnergyDTA);
00104 if (verboseLevel > 1)
00105 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
00106 << excitationEnergyDTA << G4endl;
00107
00108 incidentKineticEnergy -= excitation;
00109 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00110
00111
00112
00113
00114 G4HEVector targetParticle;
00115 if (G4UniformRand() < atomicNumber/atomicWeight) {
00116 targetParticle.setDefinition("Proton");
00117 } else {
00118 targetParticle.setDefinition("Neutron");
00119 }
00120
00121 G4double targetMass = targetParticle.getMass();
00122 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00123 + targetMass*targetMass
00124 + 2.0*targetMass*incidentTotalEnergy);
00125 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00126
00127 G4bool inElastic = true;
00128 vecLength = 0;
00129
00130 if (verboseLevel > 1)
00131 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00132 << incidentCode << G4endl;
00133
00134 G4bool successful = false;
00135
00136 FirstIntInCasXiMinus(inElastic, availableEnergy, pv, vecLength,
00137 incidentParticle, targetParticle, atomicWeight);
00138
00139 if (verboseLevel > 1)
00140 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00141
00142 if ((vecLength > 0) && (availableEnergy > 1.))
00143 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00144 pv, vecLength,
00145 incidentParticle, targetParticle);
00146
00147 HighEnergyCascading(successful, pv, vecLength,
00148 excitationEnergyGNP, excitationEnergyDTA,
00149 incidentParticle, targetParticle,
00150 atomicWeight, atomicNumber);
00151 if (!successful)
00152 HighEnergyClusterProduction(successful, pv, vecLength,
00153 excitationEnergyGNP, excitationEnergyDTA,
00154 incidentParticle, targetParticle,
00155 atomicWeight, atomicNumber);
00156 if (!successful)
00157 MediumEnergyCascading(successful, pv, vecLength,
00158 excitationEnergyGNP, excitationEnergyDTA,
00159 incidentParticle, targetParticle,
00160 atomicWeight, atomicNumber);
00161
00162 if (!successful)
00163 MediumEnergyClusterProduction(successful, pv, vecLength,
00164 excitationEnergyGNP, excitationEnergyDTA,
00165 incidentParticle, targetParticle,
00166 atomicWeight, atomicNumber);
00167 if (!successful)
00168 QuasiElasticScattering(successful, pv, vecLength,
00169 excitationEnergyGNP, excitationEnergyDTA,
00170 incidentParticle, targetParticle,
00171 atomicWeight, atomicNumber);
00172 if (!successful)
00173 ElasticScattering(successful, pv, vecLength,
00174 incidentParticle,
00175 atomicWeight, atomicNumber);
00176
00177 if (!successful)
00178 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00179 << G4endl;
00180
00181 FillParticleChange(pv, vecLength);
00182 delete [] pv;
00183 theParticleChange.SetStatusChange(stopAndKill);
00184 return &theParticleChange;
00185 }
00186
00187
00188 void
00189 G4HEXiMinusInelastic::FirstIntInCasXiMinus(G4bool& inElastic,
00190 const G4double availableEnergy,
00191 G4HEVector pv[],
00192 G4int& vecLen,
00193 const G4HEVector& incidentParticle,
00194 const G4HEVector& targetParticle,
00195 const G4double atomicWeight)
00196
00197
00198
00199
00200
00201
00202
00203
00204 {
00205 static const G4double expxu = 82.;
00206 static const G4double expxl = -expxu;
00207
00208 static const G4double protb = 0.7;
00209 static const G4double neutb = 0.7;
00210 static const G4double c = 1.25;
00211
00212 static const G4int numMul = 1200;
00213 static const G4int numSec = 60;
00214
00215 G4int neutronCode = Neutron.getCode();
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+1); 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+2); 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 == neutronCode)
00297 {
00298 if (ran < 0.2)
00299 {
00300 pv[0] = SigmaMinus;
00301 pv[1] = SigmaZero;
00302 }
00303 else if (ran < 0.4)
00304 {
00305 pv[0] = SigmaZero;
00306 pv[1] = SigmaMinus;
00307 }
00308 else if (ran < 0.6)
00309 {
00310 pv[0] = SigmaMinus;
00311 pv[1] = Lambda;
00312 }
00313 else if (ran < 0.8)
00314 {
00315 pv[0] = Lambda;
00316 pv[1] = SigmaMinus;
00317 }
00318 else
00319 {
00320 pv[0] = Neutron;
00321 pv[1] = XiMinus;
00322 }
00323 }
00324 else
00325 {
00326 if (ran < 0.2)
00327 {
00328 pv[0] = SigmaZero;
00329 pv[1] = SigmaZero;
00330 }
00331 else if (ran < 0.3)
00332 {
00333 pv[0] = Lambda;
00334 pv[1] = Lambda;
00335 }
00336 else if (ran < 0.4)
00337 {
00338 pv[0] = SigmaZero;
00339 pv[1] = Lambda;
00340 }
00341 else if (ran < 0.5)
00342 {
00343 pv[0] = Lambda;
00344 pv[1] = SigmaZero;
00345 }
00346 else if (ran < 0.7)
00347 {
00348 pv[0] = SigmaZero;
00349 pv[1] = Neutron;
00350 }
00351 else if (ran < 0.8)
00352 {
00353 pv[0] = Neutron;
00354 pv[1] = SigmaZero;
00355 }
00356 else
00357 {
00358 pv[0] = Proton;
00359 pv[1] = XiMinus;
00360 }
00361 }
00362 }
00363 return;
00364 }
00365 else if (availableEnergy <= PionPlus.getMass())
00366 return;
00367
00368
00369
00370 npos = 0; nneg = 0; nzero = 0;
00371
00372
00373
00374 G4double aleab = std::log(availableEnergy);
00375 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00376 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00377
00378
00379
00380 G4double test, dum, anpn = 0.0;
00381
00382 for (nt=1; nt<=numSec; nt++) {
00383 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00384 dum = pi*nt/(2.0*n*n);
00385 if (std::fabs(dum) < 1.0) {
00386 if (test >= 1.0e-10) anpn += dum*test;
00387 } else {
00388 anpn += dum*test;
00389 }
00390 }
00391
00392 G4double ran = G4UniformRand();
00393 G4double excs = 0.0;
00394 if( targetCode == protonCode )
00395 {
00396 counter = -1;
00397 for (npos=0; npos<numSec/3; npos++) {
00398 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00399 for (nzero=0; nzero<numSec/3; nzero++) {
00400 if (++counter < numMul) {
00401 nt = npos+nneg+nzero;
00402 if ( (nt>0) && (nt<=numSec) ) {
00403 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00404 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00405 if (std::fabs(dum) < 1.0) {
00406 if (test >= 1.0e-10) excs += dum*test;
00407 } else {
00408 excs += dum*test;
00409 }
00410 if (ran < excs) goto outOfLoop;
00411 }
00412 }
00413 }
00414 }
00415 }
00416
00417
00418 inElastic = false;
00419 return;
00420 }
00421 else
00422 {
00423 counter = -1;
00424 for (npos=0; npos<numSec/3; npos++) {
00425 for (nneg=npos; nneg<=(npos+2); nneg++) {
00426 for( nzero=0; nzero<numSec/3; nzero++) {
00427 if (++counter < numMul) {
00428 nt = npos+nneg+nzero;
00429 if ( (nt>=1) && (nt<=numSec) ) {
00430 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00431 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00432 if (std::fabs(dum) < 1.0) {
00433 if( test >= 1.0e-10 )excs += dum*test;
00434 } else {
00435 excs += dum*test;
00436 }
00437 if (ran < excs) goto outOfLoop;
00438 }
00439 }
00440 }
00441 }
00442 }
00443
00444 inElastic = false;
00445 return;
00446 }
00447
00448 outOfLoop:
00449
00450
00451
00452
00453
00454 ran = G4UniformRand();
00455 if (targetCode == neutronCode) {
00456 if( npos == nneg)
00457 {
00458 }
00459 else if (npos == (nneg-1))
00460 {
00461 if( ran < 0.50)
00462 {
00463 pv[0] = XiZero;
00464 }
00465 else
00466 {
00467 pv[1] = Proton;
00468 }
00469 }
00470 else
00471 {
00472 pv[0] = XiZero;
00473 pv[1] = Proton;
00474 }
00475 } else {
00476 if (npos == nneg)
00477 {
00478 if (ran < 0.5)
00479 {
00480 }
00481 else
00482 {
00483 pv[0] = XiZero;
00484 pv[1] = Neutron;
00485 }
00486 }
00487 else if (npos == (nneg+1))
00488 {
00489 pv[1] = Neutron;
00490 }
00491 else
00492 {
00493 pv[0] = XiZero;
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