Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLKinematicsUtils.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
38 #include "G4INCLKinematicsUtils.hh"
39 #include "G4INCLParticleTable.hh"
40 
41 namespace G4INCL {
42 
43  namespace KinematicsUtils {
44 
45  void transformToLocalEnergyFrame(Nucleus const * const n, Particle * const p) {
46  const G4double localEnergy = getLocalEnergy(n, p);
47  const G4double localTotalEnergy = p->getEnergy() - localEnergy;
48  p->setEnergy(localTotalEnergy);
50  }
51 
52  G4double getLocalEnergy(Nucleus const * const n, Particle * const p) {
53 // assert(!p->isPion()); // No local energy for pions
54 
55  G4double vloc = 0.0;
56  const G4double r = p->getPosition().mag();
57  const G4double mass = p->getMass();
58 
59  // Local energy is constant outside the surface
60  if(r > n->getUniverseRadius()) {
61  INCL_WARN("Tried to evaluate local energy for a particle outside the maximum radius."
62  << '\n' << p->print() << '\n'
63  << "Maximum radius = " << n->getDensity()->getMaximumRadius() << '\n'
64  << "Universe radius = " << n->getUniverseRadius() << '\n');
65  return 0.0;
66  }
67 
68  G4double pfl0 = 0.0;
69  const ParticleType t = p->getType();
70  const G4double kinE = p->getKineticEnergy();
71  if(kinE <= n->getPotential()->getFermiEnergy(t)) {
72  pfl0 = n->getPotential()->getFermiMomentum(p);
73  } else {
74  const G4double tf0 = p->getPotentialEnergy() - n->getPotential()->getSeparationEnergy(p);
75  if(tf0<0.0) return 0.0;
76  pfl0 = std::sqrt(tf0*(tf0 + 2.0*mass));
77  }
78  const G4double pReflection = p->getReflectionMomentum()/pfl0;
79  const G4double reflectionRadius = n->getDensity()->getMaxRFromP(p->getType(), pReflection);
80  const G4double pNominal = p->getMomentum().mag()/pfl0;
81  const G4double nominalReflectionRadius = n->getDensity()->getMaxRFromP(p->getType(), pNominal);
82  const G4double pl = pfl0*n->getDensity()->getMinPFromR(t, r*nominalReflectionRadius/reflectionRadius);
83  vloc = std::sqrt(pl*pl + mass*mass) - mass;
84 
85  return vloc;
86  }
87 
88  ThreeVector makeBoostVector(Particle const * const p1, Particle const * const p2){
89  const G4double totalEnergy = p1->getEnergy() + p2->getEnergy();
90  return ((p1->getMomentum() + p2->getMomentum())/totalEnergy);
91  }
92 
93  G4double totalEnergyInCM(Particle const * const p1, Particle const * const p2){
94  return std::sqrt(squareTotalEnergyInCM(p1,p2));
95  }
96 
97  G4double squareTotalEnergyInCM(Particle const * const p1, Particle const * const p2) {
98  G4double beta2 = makeBoostVector(p1, p2).mag2();
99  if(beta2 > 1.0) {
100  INCL_ERROR("squareTotalEnergyInCM: beta2 == " << beta2 << " > 1.0" << '\n');
101  beta2 = 0.0;
102  }
103  return (1.0 - beta2)*std::pow(p1->getEnergy() + p2->getEnergy(), 2);
104  }
105 
106  G4double momentumInCM(Particle const * const p1, Particle const * const p2) {
107  const G4double m1sq = std::pow(p1->getMass(),2);
108  const G4double m2sq = std::pow(p2->getMass(),2);
109  const G4double z = p1->getEnergy()*p2->getEnergy() - p1->getMomentum().dot(p2->getMomentum());
110  G4double pcm2 = (z*z-m1sq*m2sq)/(2*z+m1sq+m2sq);
111  if(pcm2 < 0.0) {
112  INCL_ERROR("momentumInCM: pcm2 == " << pcm2 << " < 0.0" << '\n');
113  pcm2 = 0.0;
114  }
115  return std::sqrt(pcm2);
116  }
117 
118  G4double momentumInCM(const G4double E, const G4double M1, const G4double M2) {
119  return 0.5*std::sqrt((E*E - std::pow(M1 + M2, 2))
120  *(E*E - std::pow(M1 - M2, 2)))/E;
121  }
122 
123  G4double momentumInLab(const G4double s, const G4double m1, const G4double m2) {
124  const G4double m1sq = m1*m1;
125  const G4double m2sq = m2*m2;
126  G4double plab2 = (s*s-2*s*(m1sq+m2sq)+(m1sq-m2sq)*(m1sq-m2sq))/(4*m2sq);
127  if(plab2 < 0.0) {
128  INCL_ERROR("momentumInLab: plab2 == " << plab2 << " < 0.0; m1sq == " << m1sq << "; m2sq == " << m2sq << "; s == " << s << '\n');
129  plab2 = 0.0;
130  }
131  return std::sqrt(plab2);
132  }
133 
134  G4double momentumInLab(Particle const * const p1, Particle const * const p2) {
135  const G4double m1 = p1->getMass();
136  const G4double m2 = p2->getMass();
137  const G4double s = squareTotalEnergyInCM(p1, p2);
138  return momentumInLab(s, m1, m2);
139  }
140 
142  G4double E = 0.0;
143  for(ParticleIter i=pl.begin(), e=pl.end(); i!=e; ++i) {
144  E += (*i)->getEnergy();
145  }
146  return E;
147  }
148 
150  ThreeVector p(0.0, 0.0, 0.0);
151  for(ParticleIter i=pl.begin(), e=pl.end(); i!=e; ++i) {
152  p += (*i)->getMomentum();
153  }
154  return p;
155  }
156 
158  return std::sqrt(p.mag2() + m*m);
159  }
160 
162  return std::sqrt(squareInvariantMass(E, p));
163  }
164 
166  return E*E - p.mag2();
167  }
168 
170  G4double mass;
171  if(p.theType==Composite)
173  else
175  return (1.+EKin/mass);
176  }
177 
178  }
179 
180 }
G4double getReflectionMomentum() const
Return the reflection momentum.
G4double getMass() const
Get the cached particle mass.
G4double getMinPFromR(const ParticleType t, const G4double r) const
G4ThreadLocal ParticleMassFn getTableParticleMass
Static pointer to the mass function for particles.
G4double dot(const ThreeVector &v) const
void transformToLocalEnergyFrame(Nucleus const *const n, Particle *const p)
G4double getFermiMomentum(const Particle *const p) const
Return the Fermi momentum for a particle.
G4double squareTotalEnergyInCM(Particle const *const p1, Particle const *const p2)
ThreeVector sumMomenta(const ParticleList &)
const char * p
Definition: xmltok.h:285
#define INCL_ERROR(x)
const G4INCL::ThreeVector & getMomentum() const
std::string print() const
G4double getEnergy() const
#define INCL_WARN(x)
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
ThreeVector makeBoostVector(Particle const *const p1, Particle const *const p2)
G4double mag2() const
G4double getPotentialEnergy() const
Get the particle potential energy.
void setEnergy(G4double energy)
G4double gammaFromKineticEnergy(const ParticleSpecies &p, const G4double EKin)
const XML_Char * s
Definition: expat.h:262
tuple pl
Definition: readPY.py:5
G4double getSeparationEnergy(const Particle *const p) const
Return the separation energy for a particle.
static constexpr double m
Definition: G4SIunits.hh:129
const G4int n
G4double invariantMass(const G4double E, const ThreeVector &p)
NuclearPotential::INuclearPotential const * getPotential() const
Getter for thePotential.
G4INCL::ParticleType getType() const
NuclearDensity const * getDensity() const
Getter for theDensity.
G4double energy(const ThreeVector &p, const G4double m)
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
const G4INCL::ThreeVector & getPosition() const
G4double squareInvariantMass(const G4double E, const ThreeVector &p)
G4double getMaximumRadius() const
tuple z
Definition: test.py:28
G4double getKineticEnergy() const
Get the particle kinetic energy.
G4double getUniverseRadius() const
Getter for theUniverseRadius.
G4double momentumInLab(Particle const *const p1, Particle const *const p2)
gives the momentum in the lab frame of two particles.
G4double mag() const
double G4double
Definition: G4Types.hh:76
G4double totalEnergyInCM(Particle const *const p1, Particle const *const p2)
const ThreeVector & adjustMomentumFromEnergy()
Rescale the momentum to match the total energy.
static constexpr double m2
Definition: G4SIunits.hh:130
G4double getMaxRFromP(const ParticleType t, const G4double p) const
Get the maximum allowed radius for a given momentum.
ParticleList::const_iterator ParticleIter
G4double sumTotalEnergies(const ParticleList &)
G4double getLocalEnergy(Nucleus const *const n, Particle *const p)