G4INCLCrossSections.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // INCL++ intra-nuclear cascade model
00027 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
00028 // Davide Mancusi, CEA
00029 // Alain Boudard, CEA
00030 // Sylvie Leray, CEA
00031 // Joseph Cugnon, University of Liege
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 // #include <cassert>
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     //      FUNCTION SPN(X,IND2T3,IPIT3,f17)
00064     // SIGMA(PI+ + P) IN THE (3,3) REGION
00065     // NEW FIT BY J.VANDERMEULEN  + FIT BY Th AOUST ABOVE (3,3) RES
00066     //                              CONST AT LOW AND VERY HIGH ENERGY
00067     //      COMMON/BL8/RATHR,RAMASS                                           REL21800
00068     //      integer f17
00069     // RATHR and RAMASS are always 0.0!!!
00070 
00071     G4double x = KinematicsUtils::totalEnergyInCM(particle1, particle2);
00072     if(x>10000.) return 0.0; // no cross section above this value
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.); // 5832000 = 180^3
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     // HE pi+ p and pi- n
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; // (spnpipphe(x)+spnpimphe(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     // HE and LE pi- p and pi+ n
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     // HE pi- p and pi+ n
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     /* Concerning the way we calculate the lab momentum, see the considerations
00178      * in CrossSections::elasticNNLegacy().
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 // assert(p1->isNucleon() && p2->isNucleon());
00189     const G4double sqrts = KinematicsUtils::totalEnergyInCM(p1,p2);
00190     if(sqrts < ParticleTable::effectivePionMass + 2*ParticleTable::effectiveNucleonMass + 50.) { // approximately yields INCL4.6's hard-coded threshold in collis, 2065 MeV
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 // assert(isospin==-2 || isospin==0 || isospin==2);
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) { // pp, nn
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) { // pn
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     /* The NN cross section is parametrised as a function of the lab momentum
00287      * of one of the nucleons. For NDelta or DeltaDelta, the physical
00288      * assumption is that the cross section is the same as NN *for the same
00289      * total CM energy*. Thus, we calculate s from the particles involved, and
00290      * we convert this value to the lab momentum of a nucleon *as if this were
00291      * an NN collision*.
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       //    return elasticNN(p1, p2); // New implementation
00336       return elasticNNLegacy(p1, p2); // Translated from INCL4.6 FORTRAN
00337     else
00338       return 0.0; // No pion-nucleon elastic scattering
00339   }
00340 
00341   G4double CrossSections::calculateNNDiffCrossSection(G4double pl, G4int iso) {
00342     G4double x = 0.001 * pl; // Change to GeV
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; // Should never reach this point
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     /* We compute the interaction distance from the largest of the NN cross
00381      * sections. Note that this is different from INCL4.6, which just takes the
00382      * average of the four, and will in general lead to a different geometrical
00383      * cross section.
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     /* We compute the interaction distance from the largest of the pi-N cross
00414      * sections. Note that this is different from INCL4.6, which just takes the
00415      * average of the six, and will in general lead to a different geometrical
00416      * cross section.
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 

Generated on Mon May 27 17:48:34 2013 for Geant4 by  doxygen 1.4.7