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 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00037 #include "G4INCLCrossSections.hh"
00038 #include "G4INCLKinematicsUtils.hh"
00039 #include "G4INCLParticleTable.hh"
00040 #include "G4INCLLogger.hh"
00041
00042
00043 namespace G4INCL {
00044
00045 G4double CrossSections::total(Particle const * const p1, Particle const * const p2) {
00046 G4double inelastic = 0.0;
00047 if(p1->isNucleon() && p2->isNucleon()) {
00048 inelastic = CrossSections::deltaProduction(p1, p2);
00049 } else if((p1->isNucleon() && p2->isDelta()) ||
00050 (p1->isDelta() && p2->isNucleon())) {
00051 inelastic = CrossSections::recombination(p1, p2);
00052 } else if((p1->isNucleon() && p2->isPion()) ||
00053 (p1->isPion() && p2->isNucleon())) {
00054 inelastic = CrossSections::pionNucleon(p1, p2);
00055 } else {
00056 inelastic = 0.0;
00057 }
00058
00059 return inelastic + CrossSections::elastic(p1, p2);
00060 }
00061
00062 G4double CrossSections::pionNucleon(Particle const * const particle1, Particle const * const particle2) {
00063
00064
00065
00066
00067
00068
00069
00070
00071 G4double x = KinematicsUtils::totalEnergyInCM(particle1, particle2);
00072 if(x>10000.) return 0.0;
00073
00074 G4int ipit3 = 0;
00075 G4int ind2t3 = 0;
00076 G4double ramass = 0.0;
00077
00078 if(particle1->isPion()) {
00079 ipit3 = ParticleTable::getIsospin(particle1->getType());
00080 } else if(particle2->isPion()) {
00081 ipit3 = ParticleTable::getIsospin(particle2->getType());
00082 }
00083
00084 if(particle1->isNucleon()) {
00085 ind2t3 = ParticleTable::getIsospin(particle1->getType());
00086 } else if(particle2->isNucleon()) {
00087 ind2t3 = ParticleTable::getIsospin(particle2->getType());
00088 }
00089
00090 G4double y=x*x;
00091 G4double q2=(y-1076.0*1076.0)*(y-800.0*800.0)/y/4.0;
00092 if (q2 <= 0.) {
00093 return 0.0;
00094 }
00095 G4double q3 = std::pow(std::sqrt(q2),3);
00096 G4double f3 = q3/(q3 + 5832000.);
00097 G4double spnResult = 326.5/(std::pow((x-1215.0-ramass)*2.0/(110.0-ramass), 2)+1.0);
00098 spnResult = spnResult*(1.0-5.0*ramass/1215.0);
00099 G4double cg = 4.0 + G4double(ind2t3)*G4double(ipit3);
00100 spnResult = spnResult*f3*cg/6.0;
00101
00102 if(x < 1200.0 && spnResult < 5.0) {
00103 spnResult = 5.0;
00104 }
00105
00106
00107 if(x > 1290.0) {
00108 if((ind2t3 == 1 && ipit3 == 2) || (ind2t3 == -1 && ipit3 == -2))
00109 spnResult=CrossSections::spnPiPlusPHE(x);
00110 else if((ind2t3 == 1 && ipit3 == -2) || (ind2t3 == -1 && ipit3 == 2))
00111 spnResult=CrossSections::spnPiMinusPHE(x);
00112 else if(ipit3 == 0) spnResult = (CrossSections::spnPiPlusPHE(x) + CrossSections::spnPiMinusPHE(x))/2.0;
00113 else {
00114 ERROR("Unknown configuration!" << std::endl);
00115 }
00116 }
00117
00118 return spnResult;
00119 }
00120
00121 G4double CrossSections::spnPiPlusPHE(const G4double x) {
00122
00123 if(x <= 1750.0) {
00124 return -2.33730e-06*std::pow(x, 3)+1.13819e-02*std::pow(x,2)
00125 -1.83993e+01*x+9893.4;
00126 } else if(x > 1750.0 && x <= 2175.0) {
00127 return 1.13531e-06*std::pow(x, 3)-6.91694e-03*std::pow(x, 2)
00128 +1.39907e+01*x-9360.76;
00129 } else {
00130 return -3.18087*std::log(x)+52.9784;
00131 }
00132 }
00133
00134 G4double CrossSections::spnPiMinusPHE(const G4double x) {
00135
00136 if(x <= 1475.0) {
00137 return 0.00120683*(x-1372.52)*(x-1372.52)+26.2058;
00138 } else if(x > 1475.0 && x <= 1565.0) {
00139 return 1.15873e-05*x*x+49965.6/((x-1519.59)*(x-1519.59)+2372.55);
00140 } else if(x > 1565.0 && x <= 2400.0) {
00141 return 34.0248+43262.2/((x-1681.65)*(x-1681.65)+1689.35);
00142 } else if(x > 2400.0 && x <= 7500.0) {
00143 return 3.3e-7*(x-7500.0)*(x-7500.0)+24.5;
00144 } else {
00145 return 24.5;
00146 }
00147 }
00148
00149 G4double CrossSections::recombination(Particle const * const p1, Particle const * const p2) {
00150 const G4int isospin = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
00151 if(isospin==4 || isospin==-4) return 0.0;
00152
00153 G4double s = KinematicsUtils::squareTotalEnergyInCM(p1, p2);
00154 G4double Ecm = std::sqrt(s);
00155 G4int deltaIsospin;
00156 G4double deltaMass;
00157 if(p1->isDelta()) {
00158 deltaIsospin = ParticleTable::getIsospin(p1->getType());
00159 deltaMass = p1->getMass();
00160 } else {
00161 deltaIsospin = ParticleTable::getIsospin(p2->getType());
00162 deltaMass = p2->getMass();
00163 }
00164
00165 if(Ecm <= 938.3 + deltaMass) {
00166 return 0.0;
00167 }
00168
00169 if(Ecm < 938.3 + deltaMass + 2.0) {
00170 Ecm = 938.3 + deltaMass + 2.0;
00171 s = Ecm*Ecm;
00172 }
00173
00174 const G4double x = (s - 4.*ParticleTable::effectiveNucleonMass2) /
00175 (s - std::pow(ParticleTable::effectiveNucleonMass + deltaMass, 2));
00176 const G4double y = s/(s - std::pow(deltaMass - ParticleTable::effectiveNucleonMass, 2));
00177
00178
00179
00180 const G4double pLab = KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
00181 G4double result = 0.5 * x * y * deltaProduction(isospin, pLab);
00182 result *= 3.*(32.0 + isospin * isospin * (deltaIsospin * deltaIsospin - 5))/64.0;
00183 result /= 1.0 + 0.25 * isospin * isospin;
00184 return result;
00185 }
00186
00187 G4double CrossSections::deltaProduction(Particle const * const p1, Particle const * const p2) {
00188
00189 const G4double sqrts = KinematicsUtils::totalEnergyInCM(p1,p2);
00190 if(sqrts < ParticleTable::effectivePionMass + 2*ParticleTable::effectiveNucleonMass + 50.) {
00191 return 0.0;
00192 } else {
00193 const G4double pLab = KinematicsUtils::momentumInLab(p1,p2);
00194 const G4int isospin = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
00195 return deltaProduction(isospin, pLab);
00196 }
00197 }
00198
00199 G4double CrossSections::deltaProduction(const G4int isospin, const G4double pLab) {
00200 G4double xs = 0.0;
00201
00202
00203 const G4double momentumGeV = 0.001 * pLab;
00204 if(pLab < 800.0) {
00205 return 0.0;
00206 }
00207
00208 if(isospin==2 || isospin==-2) {
00209 if(pLab >= 2000.0) {
00210 xs = (41.0 + (60.0*momentumGeV - 54.0)*std::exp(-1.2*momentumGeV) - 77.0/(momentumGeV + 1.5));
00211 } else if(pLab >= 1500.0 && pLab < 2000.0) {
00212 xs = (41.0 + 60.0*(momentumGeV - 0.9)*std::exp(-1.2*momentumGeV) - 1250.0/(momentumGeV+50.0)+ 4.0*std::pow(momentumGeV - 1.3, 2));
00213 } else if(pLab < 1500.0) {
00214 xs = (23.5 + 24.6/(1.0 + std::exp(-10.0*momentumGeV + 12.0))
00215 -1250.0/(momentumGeV +50.0)+4.0*std::pow(momentumGeV - 1.3,2));
00216 }
00217 } else if(isospin==0) {
00218 if(pLab >= 2000.0) {
00219 xs = (42.0 - 77.0/(momentumGeV + 1.5));
00220 } else if(pLab >= 1000.0 && pLab < 2000.0) {
00221 xs = (24.2 + 8.9*momentumGeV - 31.1/std::sqrt(momentumGeV));
00222 } else if(pLab < 1000.0) {
00223 xs = (33.0 + 196.0*std::sqrt(std::pow(std::abs(momentumGeV - 0.95),5))
00224 -31.1/std::sqrt(momentumGeV));
00225 }
00226 }
00227
00228 if(xs < 0.0) return 0.0;
00229 else return xs;
00230 }
00231
00232 G4double CrossSections::elasticNNHighEnergy(const G4double momentum) {
00233 return 77.0/(momentum + 1.5);
00234 }
00235
00236 G4double CrossSections::elasticProtonNeutron(const G4double momentum) {
00237 if(momentum < 0.450) {
00238 const G4double alp = std::log(momentum);
00239 return 6.3555*std::exp(-3.2481*alp-0.377*alp*alp);
00240 } else if(momentum >= 0.450 && momentum < 0.8) {
00241 return (33.0 + 196.0 * std::sqrt(std::pow(std::abs(momentum - 0.95), 5)));
00242 } else if(momentum > 2.0) {
00243 return CrossSections::elasticNNHighEnergy(momentum);
00244 } else {
00245 return 31.0/std::sqrt(momentum);
00246 }
00247 }
00248
00249 G4double CrossSections::elasticProtonProtonOrNeutronNeutron(const G4double momentum)
00250 {
00251 if(momentum < 0.440) {
00252 return 34.0*std::pow(momentum/0.4, -2.104);
00253 } else if(momentum < 0.8 && momentum >= 0.440) {
00254 return (23.5 + 1000.0*std::pow(momentum-0.7, 4));
00255 } else if(momentum < 2.0) {
00256 return (1250.0/(50.0 + momentum) - 4.0*std::pow(momentum-1.3, 2));
00257 } else {
00258 return CrossSections::elasticNNHighEnergy(momentum);
00259 }
00260 }
00261
00262 G4double CrossSections::elasticNN(Particle const * const p1, Particle const * const p2) {
00263 G4double momentum = 0.0;
00264 momentum = 0.001 * KinematicsUtils::momentumInLab(p1, p2);
00265 if((p1->getType() == Proton && p2->getType() == Proton) ||
00266 (p1->getType() == Neutron && p2->getType() == Neutron)) {
00267 return CrossSections::elasticProtonProtonOrNeutronNeutron(momentum);
00268 } else if((p1->getType() == Proton && p2->getType() == Neutron) ||
00269 (p1->getType() == Neutron && p2->getType() == Proton)) {
00270 return CrossSections::elasticProtonNeutron(momentum);
00271 } else {
00272 ERROR("G4INCL::CrossSections::elasticNN: Bad input!" << std::endl
00273 << p1->print() << p2->print() << std::endl);
00274 }
00275 return 0.0;
00276 }
00277
00278 G4double CrossSections::elasticNNLegacy(Particle const * const part1, Particle const * const part2) {
00279 G4double scale = 1.0;
00280
00281
00282 G4int i = ParticleTable::getIsospin(part1->getType())
00283 + ParticleTable::getIsospin(part2->getType());
00284 G4double sel = 0.0;
00285
00286
00287
00288
00289
00290
00291
00292
00293 const G4double s = KinematicsUtils::squareTotalEnergyInCM(part1, part2);
00294 G4double plab = KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
00295
00296 G4double p1=0.001*plab;
00297 if(plab > 2000.) goto sel13;
00298 if(part1->isNucleon() && part2->isNucleon())
00299 goto sel1;
00300 else
00301 goto sel3;
00302 sel1: if (i == 0) goto sel2;
00303 sel3: if (plab < 800.) goto sel4;
00304 if (plab > 2000.) goto sel13;
00305 sel=(1250./(50.+p1)-4.*std::pow(p1-1.3, 2))*scale;
00306 goto sel100;
00307 return sel;
00308 sel4: if (plab < 440.) {
00309 sel=34.*std::pow(p1/0.4, (-2.104))*scale;
00310 } else {
00311 sel=(23.5+1000.*std::pow(p1-0.7, 4))*scale;
00312 }
00313 goto sel100;
00314 return sel;
00315 sel13: sel=77./(p1+1.5)*scale;
00316 goto sel100;
00317 return sel;
00318 sel2: if (plab < 800.) goto sel11;
00319 if (plab > 2000.) goto sel13;
00320 sel=31./std::sqrt(p1)*scale;
00321 goto sel100;
00322 return sel;
00323 sel11: if (plab < 450.) {
00324 G4double alp=std::log(p1);
00325 sel=6.3555*std::exp(-3.2481*alp-0.377*alp*alp)*scale;
00326 } else {
00327 sel=(33.+196.*std::sqrt(std::pow(std::abs(p1-0.95),5)))*scale;
00328 }
00329
00330 sel100: return sel;
00331 }
00332
00333 G4double CrossSections::elastic(Particle const * const p1, Particle const * const p2) {
00334 if(!p1->isPion() && !p2->isPion())
00335
00336 return elasticNNLegacy(p1, p2);
00337 else
00338 return 0.0;
00339 }
00340
00341 G4double CrossSections::calculateNNDiffCrossSection(G4double pl, G4int iso) {
00342 G4double x = 0.001 * pl;
00343 if(iso != 0) {
00344 if(pl <= 2000.0) {
00345 x = std::pow(x, 8);
00346 return 5.5e-6 * x/(7.7 + x);
00347 } else {
00348 return (5.34 + 0.67*(x - 2.0)) * 1.0e-6;
00349 }
00350 } else {
00351 if(pl < 800.0) {
00352 G4double b = (7.16 - 1.63*x) * 1.0e-6;
00353 return b/(1.0 + std::exp(-(x - 0.45)/0.05));
00354 } else if(pl < 1100.0) {
00355 return (9.87 - 4.88 * x) * 1.0e-6;
00356 } else {
00357 return (3.68 + 0.76*x) * 1.0e-6;
00358 }
00359 }
00360 return 0.0;
00361 }
00362
00363 G4double CrossSections::interactionDistanceNN(const G4double projectileKineticEnergy) {
00364 ThreeVector nullVector;
00365 ThreeVector unitVector(0., 0., 1.);
00366
00367 Particle protonProjectile(Proton, unitVector, nullVector);
00368 protonProjectile.setEnergy(protonProjectile.getMass()+projectileKineticEnergy);
00369 protonProjectile.adjustMomentumFromEnergy();
00370 Particle neutronProjectile(Neutron, unitVector, nullVector);
00371 neutronProjectile.setEnergy(neutronProjectile.getMass()+projectileKineticEnergy);
00372 neutronProjectile.adjustMomentumFromEnergy();
00373
00374 Particle protonTarget(Proton, nullVector, nullVector);
00375 Particle neutronTarget(Neutron, nullVector, nullVector);
00376 const G4double sigmapp = total(&protonProjectile, &protonTarget);
00377 const G4double sigmapn = total(&protonProjectile, &neutronTarget);
00378 const G4double sigmanp = total(&neutronProjectile, &protonTarget);
00379 const G4double sigmann = total(&neutronProjectile, &neutronTarget);
00380
00381
00382
00383
00384
00385 const G4double largestSigma = std::max(sigmapp, std::max(sigmapn, std::max(sigmanp,sigmann)));
00386 const G4double interactionDistance = std::sqrt(largestSigma/Math::tenPi);
00387
00388 return interactionDistance;
00389 }
00390
00391 G4double CrossSections::interactionDistancePiN(const G4double projectileKineticEnergy) {
00392 ThreeVector nullVector;
00393 ThreeVector unitVector(0., 0., 1.);
00394
00395 Particle piPlusProjectile(PiPlus, unitVector, nullVector);
00396 piPlusProjectile.setEnergy(piPlusProjectile.getMass()+projectileKineticEnergy);
00397 piPlusProjectile.adjustMomentumFromEnergy();
00398 Particle piZeroProjectile(PiZero, unitVector, nullVector);
00399 piZeroProjectile.setEnergy(piZeroProjectile.getMass()+projectileKineticEnergy);
00400 piZeroProjectile.adjustMomentumFromEnergy();
00401 Particle piMinusProjectile(PiMinus, unitVector, nullVector);
00402 piMinusProjectile.setEnergy(piMinusProjectile.getMass()+projectileKineticEnergy);
00403 piMinusProjectile.adjustMomentumFromEnergy();
00404
00405 Particle protonTarget(Proton, nullVector, nullVector);
00406 Particle neutronTarget(Neutron, nullVector, nullVector);
00407 const G4double sigmapipp = total(&piPlusProjectile, &protonTarget);
00408 const G4double sigmapipn = total(&piPlusProjectile, &neutronTarget);
00409 const G4double sigmapi0p = total(&piZeroProjectile, &protonTarget);
00410 const G4double sigmapi0n = total(&piZeroProjectile, &neutronTarget);
00411 const G4double sigmapimp = total(&piMinusProjectile, &protonTarget);
00412 const G4double sigmapimn = total(&piMinusProjectile, &neutronTarget);
00413
00414
00415
00416
00417
00418 const G4double largestSigma = std::max(sigmapipp, std::max(sigmapipn, std::max(sigmapi0p, std::max(sigmapi0n, std::max(sigmapimp,sigmapimn)))));
00419 const G4double interactionDistance = std::sqrt(largestSigma/Math::tenPi);
00420
00421 return interactionDistance;
00422 }
00423
00424 }
00425