Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4HadronElastic.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// $Id$
27//
28// Geant4 Header : G4HadronElastic
29//
30// Author : V.Ivanchenko 29 June 2009 (redesign old elastic model)
31//
32
33#include "G4HadronElastic.hh"
34#include "G4SystemOfUnits.hh"
35#include "G4ParticleTable.hh"
36#include "G4ParticleDefinition.hh"
37#include "G4IonTable.hh"
38#include "Randomize.hh"
39#include "G4Proton.hh"
40#include "G4Neutron.hh"
41#include "G4Deuteron.hh"
42#include "G4Alpha.hh"
43#include "G4Pow.hh"
44
45G4HadronElastic::G4HadronElastic(const G4String& name)
46 : G4HadronicInteraction(name)
47{
48 SetMinEnergy( 0.0*GeV );
49 SetMaxEnergy( 100.*TeV );
50 lowestEnergyLimit= 1.e-6*eV;
51
52 theProton = G4Proton::Proton();
53 theNeutron = G4Neutron::Neutron();
54 theDeuteron = G4Deuteron::Deuteron();
55 theAlpha = G4Alpha::Alpha();
56 //Description();
57}
58
59
60G4HadronElastic::~G4HadronElastic()
61{}
62
63
64void G4HadronElastic::Description() const
65{
66 char* dirName = getenv("G4PhysListDocDir");
67 if (dirName) {
68 std::ofstream outFile;
69 G4String outFileName = GetModelName() + ".html";
70 G4String pathName = G4String(dirName) + "/" + outFileName;
71 outFile.open(pathName);
72 outFile << "<html>\n";
73 outFile << "<head>\n";
74
75 outFile << "<title>Description of G4HadronElastic Model</title>\n";
76 outFile << "</head>\n";
77 outFile << "<body>\n";
78
79 outFile << "G4HadronElastic is a hadron-nucleus elastic scattering\n"
80 << "model which uses the Gheisha two-exponential momentum\n"
81 << "transfer parameterization. The model is fully relativistic\n"
82 << "as opposed to the original Gheisha model which was not.\n"
83 << "This model may be used for all long-lived hadrons at all\n"
84 << "incident energies.\n";
85
86 outFile << "</body>\n";
87 outFile << "</html>\n";
88 outFile.close();
89 }
90}
91
92
93G4HadFinalState* G4HadronElastic::ApplyYourself(
94 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
95{
96 theParticleChange.Clear();
97
98 const G4HadProjectile* aParticle = &aTrack;
99 G4double ekin = aParticle->GetKineticEnergy();
100 if(ekin <= lowestEnergyLimit) {
101 theParticleChange.SetEnergyChange(ekin);
102 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
103 return &theParticleChange;
104 }
105
106 G4int A = targetNucleus.GetA_asInt();
107 G4int Z = targetNucleus.GetZ_asInt();
108
109 G4double plab = aParticle->GetTotalMomentum();
110
111 // Scattered particle referred to axis of incident particle
112 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
113 G4double m1 = theParticle->GetPDGMass();
114
115 if (verboseLevel>1) {
116 G4cout << "G4HadronElastic: "
117 << aParticle->GetDefinition()->GetParticleName()
118 << " Plab(GeV/c)= " << plab/GeV
119 << " Ekin(MeV) = " << ekin/MeV
120 << " scattered off Z= " << Z
121 << " A= " << A
122 << G4endl;
123 }
124
125 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
126 G4LorentzVector lv1 = aParticle->Get4Momentum();
127 G4LorentzVector lv(0.0,0.0,0.0,mass2);
128 lv += lv1;
129
130 G4ThreeVector bst = lv.boostVector();
131 lv1.boost(-bst);
132
133 G4ThreeVector p1 = lv1.vect();
134 G4double momentumCMS = p1.mag();
135 G4double tmax = 4.0*momentumCMS*momentumCMS;
136
137 // Sampling in CM system
138 G4double t = SampleInvariantT(theParticle, plab, Z, A);
139 G4double phi = G4UniformRand()*CLHEP::twopi;
140 G4double cost = 1. - 2.0*t/tmax;
141 G4double sint;
142
143 // problem in sampling
144 if(cost > 1.0 || cost < -1.0) {
145 //if(verboseLevel > 0) {
146 G4cout << "G4HadronElastic WARNING (1 - cost)= " << 1 - cost
147 << " after scattering of "
148 << aParticle->GetDefinition()->GetParticleName()
149 << " p(GeV/c)= " << plab/GeV
150 << " on an ion Z= " << Z << " A= " << A
151 << G4endl;
152 //}
153 cost = 1.0;
154 sint = 0.0;
155
156 // normal situation
157 } else {
158 sint = std::sqrt((1.0-cost)*(1.0+cost));
159 }
160 if (verboseLevel>1) {
161 G4cout << " t= " << t << " tmax(GeV^2)= " << tmax/(GeV*GeV)
162 << " Pcms(GeV)= " << momentumCMS/GeV << " cos(t)=" << cost
163 << " sin(t)=" << sint << G4endl;
164 }
165 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
166 v1 *= momentumCMS;
167 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),
168 std::sqrt(momentumCMS*momentumCMS + m1*m1));
169
170 nlv1.boost(bst);
171
172 G4double eFinal = nlv1.e() - m1;
173 if (verboseLevel > 1) {
174 G4cout <<" m= " << m1 << " Efin(MeV)= " << eFinal
175 << " Proj: 4-mom " << lv1 << " Final: " << nlv1
176 << G4endl;
177 }
178 if(eFinal <= lowestEnergyLimit) {
179 if(eFinal < 0.0 && verboseLevel > 0) {
180 G4cout << "G4HadronElastic WARNING Efinal= " << eFinal
181 << " after scattering of "
182 << aParticle->GetDefinition()->GetParticleName()
183 << " p(GeV/c)= " << plab/GeV
184 << " on an ion Z= " << Z << " A= " << A
185 << G4endl;
186 }
187 theParticleChange.SetEnergyChange(0.0);
188 nlv1 = G4LorentzVector(0.0,0.0,0.0,m1);
189
190 } else {
191 theParticleChange.SetMomentumChange(nlv1.vect().unit());
192 theParticleChange.SetEnergyChange(eFinal);
193 }
194
195 lv -= nlv1;
196 G4double erec = lv.e() - mass2;
197 if (verboseLevel > 1) {
198 G4cout << "Recoil: " <<" m= " << mass2 << " Erec(MeV)= " << erec
199 << " 4-mom: " << lv
200 << G4endl;
201 }
202
203 if(erec > GetRecoilEnergyThreshold()) {
204 G4ParticleDefinition * theDef = 0;
205 if(Z == 1 && A == 1) { theDef = theProton; }
206 else if (Z == 1 && A == 2) { theDef = theDeuteron; }
207 else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
208 else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
209 else if (Z == 2 && A == 4) { theDef = theAlpha; }
210 else {
211 theDef =
212 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
213 }
214 G4DynamicParticle * aSec = new G4DynamicParticle(theDef, lv);
215 theParticleChange.AddSecondary(aSec);
216 } else if(erec > 0.0) {
217 theParticleChange.SetLocalEnergyDeposit(erec);
218 }
219
220 return &theParticleChange;
221}
222
223// sample momentum transfer in the CMS system
224G4double
225G4HadronElastic::SampleInvariantT(const G4ParticleDefinition* p,
226 G4double plab,
227 G4int Z, G4int A)
228{
229 static const G4double GeV2 = GeV*GeV;
230 G4double momentumCMS = ComputeMomentumCMS(p,plab,Z,A);
231 G4double tmax = 4.0*momentumCMS*momentumCMS/GeV2;
232 G4double aa, bb, cc;
233 G4double dd = 10.;
234 G4Pow* g4pow = G4Pow::GetInstance();
235 if (A <= 62) {
236 bb = 14.5*g4pow->Z23(A);
237 aa = g4pow->powZ(A, 1.63)/bb;
238 cc = 1.4*g4pow->Z13(A)/dd;
239 } else {
240 bb = 60.*g4pow->Z13(A);
241 aa = g4pow->powZ(A, 1.33)/bb;
242 cc = 0.4*g4pow->powZ(A, 0.4)/dd;
243 }
244 G4double q1 = 1.0 - std::exp(-bb*tmax);
245 G4double q2 = 1.0 - std::exp(-dd*tmax);
246 G4double s1 = q1*aa;
247 G4double s2 = q2*cc;
248 if((s1 + s2)*G4UniformRand() < s2) {
249 q1 = q2;
250 bb = dd;
251 }
252 return -GeV2*std::log(1.0 - G4UniformRand()*q1)/bb;
253}
254
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::mag
double mag() const
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4HadFinalState::AddSecondary
void AddSecondary(G4DynamicParticle *aP)
Definition: G4HadFinalState.cc:70
G4HadFinalState::SetEnergyChange
void SetEnergyChange(G4double anEnergy)
Definition: G4HadFinalState.cc:42
G4HadFinalState::SetMomentumChange
void SetMomentumChange(const G4ThreeVector &aV)
Definition: G4HadFinalState.cc:54
G4HadFinalState::SetLocalEnergyDeposit
void SetLocalEnergyDeposit(G4double aE)
Definition: G4HadFinalState.cc:129
G4HadProjectile::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4HadProjectile.hh:101
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:86
G4HadronElastic::ApplyYourself
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4HadronElastic.cc:93
G4HadronElastic::Description
virtual void Description() const
Definition: G4HadronElastic.cc:64
G4HadronElastic::ComputeMomentumCMS
G4double ComputeMomentumCMS(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4HadronElastic.hh:98
G4HadronElastic::SampleInvariantT
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4HadronElastic.cc:225
G4HadronElastic::G4HadronElastic
G4HadronElastic(const G4String &name="hElasticLHEP")
Definition: G4HadronElastic.cc:45
G4HadronElastic::~G4HadronElastic
virtual ~G4HadronElastic()
Definition: G4HadronElastic.cc:60
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:90
G4HadronicInteraction::GetRecoilEnergyThreshold
G4double GetRecoilEnergyThreshold() const
Definition: G4HadronicInteraction.hh:145
G4HadronicInteraction::GetModelName
const G4String & GetModelName() const
Definition: G4HadronicInteraction.hh:119
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:103
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int J=0)
Definition: G4IonTable.cc:267
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:53
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:115
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:63
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable()
G4Pow
Definition: G4Pow.hh:54
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:50
G4Pow::Z23
G4double Z23(G4int Z)
Definition: G4Pow.hh:134
G4Pow::Z13
G4double Z13(G4int Z)
Definition: G4Pow.hh:110
G4Pow::powZ
G4double powZ(G4int Z, G4double y)
Definition: G4Pow.hh:180
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95