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 "G4INCLKinematicsUtils.hh"
00038 #include "G4INCLParticleTable.hh"
00039
00040 namespace G4INCL {
00041
00042 void KinematicsUtils::transformToLocalEnergyFrame(Nucleus const * const n, Particle * const p) {
00043 const G4double localEnergy = KinematicsUtils::getLocalEnergy(n, p);
00044 const G4double localTotalEnergy = p->getEnergy() - localEnergy;
00045 p->setEnergy(localTotalEnergy);
00046 p->adjustMomentumFromEnergy();
00047 }
00048
00049 G4double KinematicsUtils::getLocalEnergy(Nucleus const * const n, Particle * const p) {
00050
00051
00052 G4double vloc = 0.0;
00053 const G4double r = p->getPosition().mag();
00054 const G4double mass = p->getMass();
00055
00056
00057 if(r > n->getUniverseRadius()) {
00058 WARN("Tried to evaluate local energy for a particle outside the maximum radius."
00059 << std::endl << p->print() << std::endl
00060 << "Maximum radius = " << n->getDensity()->getMaximumRadius() << std::endl
00061 << "Universe radius = " << n->getUniverseRadius() << std::endl);
00062 return 0.0;
00063 }
00064
00065 G4double pfl0 = 0.0;
00066 const G4double kinE = p->getKineticEnergy();
00067 if(kinE <= n->getPotential()->getFermiEnergy(p->getType())) {
00068 pfl0 = n->getPotential()->getFermiMomentum(p);
00069 } else {
00070 const G4double tf0 = p->getPotentialEnergy() - n->getPotential()->getSeparationEnergy(p);
00071 if(tf0<0.0) return 0.0;
00072 pfl0 = std::sqrt(tf0*(tf0 + 2.0*mass));
00073 }
00074 const G4double pl = pfl0*n->getDensity()->getMaxTFromR(r);
00075 vloc = std::sqrt(pl*pl + mass*mass) - mass;
00076
00077 return vloc;
00078 }
00079
00080 ThreeVector KinematicsUtils::makeBoostVector(Particle const * const p1, Particle const * const p2){
00081 const G4double totalEnergy = p1->getEnergy() + p2->getEnergy();
00082 return ((p1->getMomentum() + p2->getMomentum())/totalEnergy);
00083 }
00084
00085 G4double KinematicsUtils::totalEnergyInCM(Particle const * const p1, Particle const * const p2){
00086 return std::sqrt(squareTotalEnergyInCM(p1,p2));
00087 }
00088
00089 G4double KinematicsUtils::squareTotalEnergyInCM(Particle const * const p1, Particle const * const p2) {
00090 G4double beta2 = KinematicsUtils::makeBoostVector(p1, p2).mag2();
00091 if(beta2 > 1.0) {
00092 ERROR("KinematicsUtils::squareTotalEnergyInCM: beta2 == " << beta2 << " > 1.0" << std::endl);
00093 beta2 = 0.0;
00094 }
00095 return (1.0 - beta2)*std::pow(p1->getEnergy() + p2->getEnergy(), 2);
00096 }
00097
00098 G4double KinematicsUtils::momentumInCM(Particle const * const p1, Particle const * const p2) {
00099 const G4double m1sq = std::pow(p1->getMass(),2);
00100 const G4double m2sq = std::pow(p2->getMass(),2);
00101 const G4double z = p1->getEnergy()*p2->getEnergy() - p1->getMomentum().dot(p2->getMomentum());
00102 G4double pcm2 = (z*z-m1sq*m2sq)/(2*z+m1sq+m2sq);
00103 if(pcm2 < 0.0) {
00104 ERROR("KinematicsUtils::momentumInCM: pcm2 == " << pcm2 << " < 0.0" << std::endl);
00105 pcm2 = 0.0;
00106 }
00107 return std::sqrt(pcm2);
00108 }
00109
00110 G4double KinematicsUtils::momentumInCM(const G4double E, const G4double M1, const G4double M2) {
00111 return 0.5*std::sqrt((E*E - std::pow(M1 + M2, 2))
00112 *(E*E - std::pow(M1 - M2, 2)))/E;
00113 }
00114
00115 G4double KinematicsUtils::momentumInLab(const G4double s, const G4double m1, const G4double m2) {
00116 const G4double m1sq = m1*m1;
00117 const G4double m2sq = m2*m2;
00118 G4double plab2 = (s*s-2*s*(m1sq+m2sq)+(m1sq-m2sq)*(m1sq-m2sq))/(4*m2sq);
00119 if(plab2 < 0.0) {
00120 ERROR("KinematicsUtils::momentumInLab: plab2 == " << plab2 << " < 0.0; m1sq == " << m1sq << "; m2sq == " << m2sq << "; s == " << s << std::endl);
00121 plab2 = 0.0;
00122 }
00123 return std::sqrt(plab2);
00124 }
00125
00126 G4double KinematicsUtils::momentumInLab(Particle const * const p1, Particle const * const p2) {
00127 const G4double m1 = p1->getMass();
00128 const G4double m2 = p2->getMass();
00129 const G4double s = squareTotalEnergyInCM(p1, p2);
00130 return KinematicsUtils::momentumInLab(s, m1, m2);
00131 }
00132
00133 G4double KinematicsUtils::sumTotalEnergies(const ParticleList &pl) {
00134 G4double E = 0.0;
00135 for(ParticleIter i = pl.begin(); i != pl.end(); ++i) {
00136 E += (*i)->getEnergy();
00137 }
00138 return E;
00139 }
00140
00141 ThreeVector KinematicsUtils::sumMomenta(const ParticleList &pl) {
00142 ThreeVector p(0.0, 0.0, 0.0);
00143 for(ParticleIter i = pl.begin(); i != pl.end(); ++i) {
00144 p += (*i)->getMomentum();
00145 }
00146 return p;
00147 }
00148
00149 G4double KinematicsUtils::energy(const ThreeVector &p, const G4double m) {
00150 return std::sqrt(p.mag2() + m*m);
00151 }
00152
00153 G4double KinematicsUtils::invariantMass(const G4double E, const ThreeVector & p) {
00154 return std::sqrt(E*E - p.mag2());
00155 }
00156
00157 G4double KinematicsUtils::gammaFromKineticEnergy(const ParticleSpecies &p, const G4double EKin) {
00158 G4double mass;
00159 if(p.theType==Composite)
00160 mass = ParticleTable::getTableMass(p.theA, p.theZ);
00161 else
00162 mass = ParticleTable::getTableParticleMass(p.theType);
00163 return (1.+EKin/mass);
00164 }
00165
00166 }