Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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