Geant4  10.00.p03
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: G4HadronElastic.cc 66892 2013-01-17 10:57:59Z gunter $
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 
45 G4HadronElastic::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 
60 G4HadronElastic::~G4HadronElastic()
61 {}
62 
63 
64 void 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 
93 G4HadFinalState* 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
224 G4double
225 G4HadronElastic::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 
G4HadronElastic::ApplyYourself
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4HadronElastic.cc:93
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
MeV
static const double MeV
Definition: G4SIunits.hh:193
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:53
G4HadronElastic::SampleInvariantT
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4HadronElastic.cc:225
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Pow
Definition: G4Pow.hh:56
outFile
std::ofstream outFile
Definition: GammaRayTel.cc:68
G4HadronElastic::G4HadronElastic
G4HadronElastic(const G4String &name="hElasticLHEP")
Definition: G4HadronElastic.cc:45
G4HadronElastic.hh
name
G4String name
Definition: TRTMaterials.hh:40
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:462
G4HadronElastic::theNeutron
G4ParticleDefinition * theNeutron
Definition: G4HadronElastic.hh:79
G4HadronicInteraction::GetModelName
const G4String & GetModelName() const
Definition: G4HadronicInteraction.hh:119
G4HadronElastic::theProton
G4ParticleDefinition * theProton
Definition: G4HadronElastic.hh:78
G4int
int G4int
Definition: G4Types.hh:78
G4HadronicInteraction::GetRecoilEnergyThreshold
G4double GetRecoilEnergyThreshold() const
Definition: G4HadronicInteraction.hh:145
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4HadronElastic::Description
virtual void Description() const
Definition: G4HadronElastic.cc:64
G4IonTable.hh
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:90
G4HadronElastic::theAlpha
G4ParticleDefinition * theAlpha
Definition: G4HadronElastic.hh:81
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable() const
Definition: G4ParticleTable.cc:712
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4Pow::Z13
G4double Z13(G4int Z) const
Definition: G4Pow.hh:129
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
G4ParticleTable.hh
G4Neutron.hh
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
G4Proton.hh
G4Alpha.hh
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
GeV
static const double GeV
Definition: G4SIunits.hh:196
G4ParticleDefinition.hh
Randomize.hh
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
A
static const G4double A[nN]
Definition: G4ElectroNuclearCrossSection.cc:115
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:86
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4Deuteron.hh
eV
static const double eV
Definition: G4SIunits.hh:194
G4HadFinalState::SetEnergyChange
void SetEnergyChange(G4double anEnergy)
Definition: G4HadFinalState.cc:42
G4HadronElastic::ComputeMomentumCMS
G4double ComputeMomentumCMS(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4HadronElastic.hh:98
G4HadronElastic::theDeuteron
G4ParticleDefinition * theDeuteron
Definition: G4HadronElastic.hh:80
G4HadronElastic::lowestEnergyLimit
G4double lowestEnergyLimit
Definition: G4HadronElastic.hh:83
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:95
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4Pow::Z23
G4double Z23(G4int Z) const
Definition: G4Pow.hh:153
G4HadronElastic::~G4HadronElastic
virtual ~G4HadronElastic()
Definition: G4HadronElastic.cc:60
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:103
G4HadFinalState::SetLocalEnergyDeposit
void SetLocalEnergyDeposit(G4double aE)
Definition: G4HadFinalState.cc:129
G4endl
#define G4endl
Definition: G4ios.hh:61
TeV
static const double TeV
Definition: G4SIunits.hh:197
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4Pow::powZ
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:258
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4Pow.hh
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4HadFinalState::SetMomentumChange
void SetMomentumChange(const G4ThreeVector &aV)
Definition: G4HadFinalState.cc:54
G4HadFinalState::AddSecondary
void AddSecondary(G4DynamicParticle *aP)
Definition: G4HadFinalState.cc:70
G4HadProjectile::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4HadProjectile.hh:101
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35