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 "G4HEXiZeroInelastic.hh"
00037 #include "globals.hh"
00038 #include "G4ios.hh"
00039 #include "G4PhysicalConstants.hh"
00040
00041 void G4HEXiZeroInelastic::ModelDescription(std::ostream& outFile) const
00042 {
00043 outFile << "G4HEXiZeroInelastic is one of the High Energy\n"
00044 << "Parameterized (HEP) models used to implement inelastic\n"
00045 << "Xi0 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 Xi0 with initial energies\n"
00051 << "above 20 GeV.\n";
00052 }
00053
00054
00055 G4HadFinalState*
00056 G4HEXiZeroInelastic::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 << "GHEXiZeroInelastic: incident energy < 1 GeV" << G4endl;
00079
00080 if (verboseLevel > 1) {
00081 G4cout << "G4HEXiZeroInelastic::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 FirstIntInCasXiZero(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 G4HEXiZeroInelastic::FirstIntInCasXiZero(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 protonCode = Proton.getCode();
00216 G4int targetCode = targetParticle.getCode();
00217 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00218
00219 static G4bool first = true;
00220 static G4double protmul[numMul], protnorm[numSec];
00221 static G4double neutmul[numMul], neutnorm[numSec];
00222
00223
00224
00225
00226 G4int i, counter, nt, npos, nneg, nzero;
00227
00228 if (first) {
00229
00230 first = false;
00231 for (i = 0; i < numMul; i++) protmul[i] = 0.0;
00232 for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
00233 counter = -1;
00234 for (npos = 0; npos < (numSec/3); npos++) {
00235 for (nneg = std::max(0,npos-2); nneg <= npos; nneg++) {
00236 for (nzero = 0; nzero < numSec/3; nzero++) {
00237 if (++counter < numMul) {
00238 nt = npos+nneg+nzero;
00239 if ((nt>0) && (nt<=numSec) ) {
00240 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00241 protnorm[nt-1] += protmul[counter];
00242 }
00243 }
00244 }
00245 }
00246 }
00247
00248 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00249 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00250 counter = -1;
00251 for( npos=0; npos<numSec/3; npos++ )
00252 {
00253 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
00254 {
00255 for( nzero=0; nzero<numSec/3; nzero++ )
00256 {
00257 if( ++counter < numMul )
00258 {
00259 nt = npos+nneg+nzero;
00260 if( (nt>0) && (nt<=numSec) )
00261 {
00262 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00263 neutnorm[nt-1] += neutmul[counter];
00264 }
00265 }
00266 }
00267 }
00268 }
00269 for (i = 0; i < numSec; i++) {
00270 if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i];
00271 if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i];
00272 }
00273 }
00274
00275
00276 pv[0] = incidentParticle;
00277 pv[1] = targetParticle;
00278 vecLen = 2;
00279
00280 if( !inElastic )
00281 {
00282 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00283 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
00284 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
00285 {
00286 G4double ran = G4UniformRand();
00287 if( targetCode == protonCode)
00288 {
00289 if (ran < 0.2)
00290 {
00291 pv[0] = SigmaPlus;
00292 pv[1] = SigmaZero;
00293 }
00294 else if (ran < 0.4)
00295 {
00296 pv[0] = SigmaZero;
00297 pv[1] = SigmaPlus;
00298 }
00299 else if (ran < 0.6)
00300 {
00301 pv[0] = SigmaPlus;
00302 pv[1] = Lambda;
00303 }
00304 else if (ran < 0.8)
00305 {
00306 pv[0] = Lambda;
00307 pv[1] = SigmaPlus;
00308 }
00309 else
00310 {
00311 pv[0] = Proton;
00312 pv[1] = XiZero;
00313 }
00314 }
00315 else
00316 {
00317 if (ran < 0.2)
00318 {
00319 pv[0] = Neutron;
00320 pv[1] = XiZero;
00321 }
00322 else if (ran < 0.3)
00323 {
00324 pv[0] = SigmaZero;
00325 pv[1] = SigmaZero;
00326 }
00327 else if (ran < 0.4)
00328 {
00329 pv[0] = Lambda;
00330 pv[1] = Lambda;
00331 }
00332 else if (ran < 0.5)
00333 {
00334 pv[0] = SigmaZero;
00335 pv[1] = Lambda;
00336 }
00337 else if (ran < 0.6)
00338 {
00339 pv[0] = Lambda;
00340 pv[1] = SigmaZero;
00341 }
00342 else if (ran < 0.7)
00343 {
00344 pv[0] = SigmaPlus;
00345 pv[1] = SigmaMinus;
00346 }
00347 else if (ran < 0.8)
00348 {
00349 pv[0] = SigmaMinus;
00350 pv[1] = SigmaPlus;
00351 }
00352 else if (ran < 0.9)
00353 {
00354 pv[0] = XiMinus;
00355 pv[1] = Proton;
00356 }
00357 else
00358 {
00359 pv[0] = Proton;
00360 pv[1] = XiMinus;
00361 }
00362 }
00363 }
00364 return;
00365 }
00366 else if (availableEnergy <= PionPlus.getMass())
00367 return;
00368
00369
00370
00371 npos = 0; nneg = 0; nzero = 0;
00372
00373
00374
00375 G4double aleab = std::log(availableEnergy);
00376 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00377 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00378
00379
00380
00381 G4double test, dum, anpn = 0.0;
00382
00383 for (nt=1; nt<=numSec; nt++) {
00384 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00385 dum = pi*nt/(2.0*n*n);
00386 if (std::fabs(dum) < 1.0) {
00387 if (test >= 1.0e-10) anpn += dum*test;
00388 } else {
00389 anpn += dum*test;
00390 }
00391 }
00392
00393 G4double ran = G4UniformRand();
00394 G4double excs = 0.0;
00395 if( targetCode == protonCode )
00396 {
00397 counter = -1;
00398 for (npos=0; npos<numSec/3; npos++) {
00399 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
00400 for (nzero=0; nzero<numSec/3; nzero++) {
00401 if (++counter < numMul) {
00402 nt = npos+nneg+nzero;
00403 if ( (nt>0) && (nt<=numSec) ) {
00404 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00405 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00406 if (std::fabs(dum) < 1.0) {
00407 if (test >= 1.0e-10) excs += dum*test;
00408 } else {
00409 excs += dum*test;
00410 }
00411 if (ran < excs) goto outOfLoop;
00412 }
00413 }
00414 }
00415 }
00416 }
00417
00418
00419
00420 inElastic = false;
00421 return;
00422 }
00423 else
00424 {
00425 counter = -1;
00426 for (npos=0; npos<numSec/3; npos++) {
00427 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00428 for (nzero=0; nzero<numSec/3; nzero++) {
00429 if (++counter < numMul) {
00430 nt = npos+nneg+nzero;
00431 if ( (nt>=1) && (nt<=numSec) ) {
00432 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00433 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00434 if (std::fabs(dum) < 1.0) {
00435 if (test >= 1.0e-10) excs += dum*test;
00436 } else {
00437 excs += dum*test;
00438 }
00439 if (ran < excs) goto outOfLoop;
00440 }
00441 }
00442 }
00443 }
00444 }
00445
00446
00447 inElastic = false;
00448 return;
00449 }
00450
00451 outOfLoop:
00452
00453
00454
00455
00456
00457 ran = G4UniformRand();
00458 if (targetCode == protonCode) {
00459 if( npos == nneg)
00460 {
00461 }
00462 else if (npos == (nneg+1))
00463 {
00464 if( ran < 0.50)
00465 {
00466 pv[0] = XiMinus;
00467 }
00468 else
00469 {
00470 pv[1] = Neutron;
00471 }
00472 }
00473 else
00474 {
00475 pv[0] = XiMinus;
00476 pv[1] = Neutron;
00477 }
00478 } else {
00479 if (npos == nneg)
00480 {
00481 if (ran < 0.5)
00482 {
00483 }
00484 else
00485 {
00486 pv[0] = XiMinus;
00487 pv[1] = Proton;
00488 }
00489 }
00490 else if (npos == (nneg-1))
00491 {
00492 pv[1] = Proton;
00493 }
00494 else
00495 {
00496 pv[0] = XiMinus;
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].getCode() << " " ;
00525 G4cout << pv[1].getCode() << " " ;
00526 for ( i = 2; i < vecLen; i++) G4cout << pv[i].getCode() << " " ;
00527 G4cout << G4endl;
00528 }
00529
00530 return;
00531 }
00532