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
G4WHadronElasticProcess.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 Hadron Elastic Scattering Process
29//
30// Created 21 April 2006 V.Ivanchenko
31//
32// Modified:
33// 24.04.06 V.Ivanchenko add neutron scattering on hydrogen from CHIPS
34// 07.06.06 V.Ivanchenko fix problem of rotation of final state
35// 25.07.06 V.Ivanchenko add 19 MeV low energy for CHIPS
36// 26.09.06 V.Ivanchenko add lowestEnergy
37// 20.10.06 V.Ivanchenko initialise lowestEnergy=0 for neitrals, eV for charged
38// 23.01.07 V.Ivanchenko add cross section interfaces with Z and A
39// 02.05.07 V.Ivanchenko add He3
40// 13.01.10 M.Kosov: Commented not used G4QElasticCrossSection & G4QCHIPSWorld
41// 24.02.11 V.Ivanchenko use particle name in IfApplicable,
42// added anti particles for light ions
43// 07.09.11 M.Kelsey: Follow chanhe to G4HadFinalState interface
44// 14.09.12 M.Kelsey: Pass subType code to base ctor
45//
46
47#include <iostream>
48#include <typeinfo>
49
50#include "G4WHadronElasticProcess.hh"
51#include "G4SystemOfUnits.hh"
52#include "G4Nucleus.hh"
53#include "G4ProcessManager.hh"
54#include "G4CrossSectionDataStore.hh"
55#include "G4HadronElasticDataSet.hh"
56#include "G4ProductionCutsTable.hh"
57#include "G4HadronicException.hh"
58#include "G4HadronicDeprecate.hh"
59
60G4WHadronElasticProcess::G4WHadronElasticProcess(const G4String& pName)
61 : G4HadronicProcess(pName, fHadronElastic) {
62 AddDataSet(new G4HadronElasticDataSet);
63 theNeutron = G4Neutron::Neutron();
64 lowestEnergy = 1.*keV;
65 lowestEnergyNeutron = 1.e-6*eV;
66 G4HadronicDeprecate("G4WHadronElasticProcess");
67}
68
69G4WHadronElasticProcess::~G4WHadronElasticProcess()
70{}
71
72void G4WHadronElasticProcess::Description() const
73{
74 char* dirName = getenv("G4PhysListDocDir");
75 if (dirName) {
76 std::ofstream outFile;
77 G4String outFileName = GetProcessName() + ".html";
78 G4String pathName = G4String(dirName) + "/" + outFileName;
79 outFile.open(pathName);
80 outFile << "<html>\n";
81 outFile << "<head>\n";
82
83 outFile << "<title>Description of G4WHadronElasticProcess</title>\n";
84 outFile << "</head>\n";
85 outFile << "<body>\n";
86
87 outFile << "G4WHadronElasticProcess handles the elastic scattering of\n"
88 << "hadrons by invoking one or more hadronic models and one or\n"
89 << "more hadronic cross sections.\n";
90
91 outFile << "</body>\n";
92 outFile << "</html>\n";
93 outFile.close();
94 }
95}
96
97
98G4VParticleChange*
99G4WHadronElasticProcess::PostStepDoIt(const G4Track& track,
100 const G4Step& /*step*/)
101{
102 theTotalResult->Clear();
103 theTotalResult->Initialize(track);
104 G4double weight = track.GetWeight();
105 theTotalResult->ProposeWeight(weight);
106
107 // For elastic scattering, _any_ result is considered an interaction
108 ClearNumberOfInteractionLengthLeft();
109
110 G4double kineticEnergy = track.GetKineticEnergy();
111 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
112 const G4ParticleDefinition* part = dynParticle->GetDefinition();
113
114 // NOTE: Very low energy scatters were causing numerical (FPE) errors
115 // in earlier releases; these limits have not been changed since.
116 if (part == theNeutron) {
117 if(kineticEnergy <= lowestEnergyNeutron) return theTotalResult;
118 } else if(kineticEnergy <= lowestEnergy) return theTotalResult;
119
120 G4Material* material = track.GetMaterial();
121 G4Nucleus* targNucleus = GetTargetNucleusPointer();
122
123 // Select element
124 G4Element* elm = 0;
125 try
126 {
127 elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,
128 *targNucleus);
129 }
130 catch(G4HadronicException & aR)
131 {
132 G4ExceptionDescription ed;
133 DumpState(track,"SampleZandA",ed);
134 ed << " PostStepDoIt failed on element selection" << G4endl;
135 G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had003",
136 FatalException, ed);
137 }
138 G4HadronicInteraction* hadi = 0;
139 try
140 {
141 hadi = ChooseHadronicInteraction( kineticEnergy, material, elm );
142 }
143 catch(G4HadronicException & aE)
144 {
145 G4ExceptionDescription ed;
146 ed << "Target element "<< elm->GetName()<<" Z= "
147 << targNucleus->GetZ_asInt() << " A= "
148 << targNucleus->GetA_asInt() << G4endl;
149 DumpState(track,"ChooseHadronicInteraction",ed);
150 ed << " No HadronicInteraction found out" << G4endl;
151 G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had005",
152 FatalException, ed);
153 }
154
155 size_t idx = track.GetMaterialCutsCouple()->GetIndex();
156 G4double tcut = (*(G4ProductionCutsTable::GetProductionCutsTable()
157 ->GetEnergyCutsVector(3)))[idx];
158 hadi->SetRecoilEnergyThreshold(tcut);
159
160 // Initialize the hadronic projectile from the track
161 // G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
162 G4HadProjectile theProj(track);
163 if(verboseLevel>1) {
164 G4cout << "G4WHadronElasticProcess::PostStepDoIt for "
165 << part->GetParticleName()
166 << " in " << material->GetName()
167 << " Target Z= " << targNucleus->GetZ_asInt()
168 << " A= " << targNucleus->GetA_asInt() << G4endl;
169 }
170
171 G4HadFinalState* result = 0;
172 try
173 {
174 result = hadi->ApplyYourself( theProj, *targNucleus);
175 }
176 catch(G4HadronicException aR)
177 {
178 G4ExceptionDescription ed;
179 ed << "Call for " << hadi->GetModelName() << G4endl;
180 ed << "Target element "<< elm->GetName()<<" Z= "
181 << targNucleus->GetZ_asInt()
182 << " A= " << targNucleus->GetA_asInt() << G4endl;
183 DumpState(track,"ApplyYourself",ed);
184 ed << " ApplyYourself failed" << G4endl;
185 G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had006",
186 FatalException, ed);
187 }
188
189 // Check the result for catastrophic energy non-conservation
190 // cannot be applied because is not guranteed that recoil
191 // nucleus is created
192 // result = CheckResult(theProj, targNucleus, result);
193
194 // directions
195 G4ThreeVector indir = track.GetMomentumDirection();
196 G4double phi = CLHEP::twopi*G4UniformRand();
197 G4ThreeVector it(0., 0., 1.);
198 G4ThreeVector outdir = result->GetMomentumChange();
199
200 if(verboseLevel>1) {
201 G4cout << "Efin= " << result->GetEnergyChange()
202 << " de= " << result->GetLocalEnergyDeposit()
203 << " nsec= " << result->GetNumberOfSecondaries()
204 << " dir= " << outdir
205 << G4endl;
206 }
207
208 // energies
209 G4double edep = result->GetLocalEnergyDeposit();
210 G4double efinal = result->GetEnergyChange();
211 if(efinal < 0.0) { efinal = 0.0; }
212 if(edep < 0.0) { edep = 0.0; }
213
214 // NOTE: Very low energy scatters were causing numerical (FPE) errors
215 // in earlier releases; these limits have not been changed since.
216 if((part == theNeutron && efinal <= lowestEnergyNeutron) ||
217 (part != theNeutron && efinal <= lowestEnergy)) {
218 edep += efinal;
219 efinal = 0.0;
220 }
221
222 // primary change
223 theTotalResult->ProposeEnergy(efinal);
224
225 G4TrackStatus status = track.GetTrackStatus();
226 if(efinal > 0.0) {
227 outdir.rotate(phi, it);
228 outdir.rotateUz(indir);
229 theTotalResult->ProposeMomentumDirection(outdir);
230 } else {
231 if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
232 { status = fStopButAlive; }
233 else { status = fStopAndKill; }
234 theTotalResult->ProposeTrackStatus(status);
235 }
236
237 //G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
238
239 theTotalResult->SetNumberOfSecondaries(0);
240
241 // recoil
242 if(result->GetNumberOfSecondaries() > 0) {
243 G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
244
245 if(p->GetKineticEnergy() > lowestEnergy) {
246 theTotalResult->SetNumberOfSecondaries(1);
247 G4ThreeVector pdir = p->GetMomentumDirection();
248 // G4cout << "recoil " << pdir << G4endl;
249 //!! is not needed for models inheriting G4HadronElastic
250 pdir.rotate(phi, it);
251 pdir.rotateUz(indir);
252 // G4cout << "recoil rotated " << pdir << G4endl;
253 p->SetMomentumDirection(pdir);
254
255 // in elastic scattering time and weight are not changed
256 G4Track* t = new G4Track(p, track.GetGlobalTime(),
257 track.GetPosition());
258 t->SetWeight(weight);
259 t->SetTouchableHandle(track.GetTouchableHandle());
260 theTotalResult->AddSecondary(t);
261
262 } else {
263 edep += p->GetKineticEnergy();
264 delete p;
265 }
266 }
267 theTotalResult->ProposeLocalEnergyDeposit(edep);
268 theTotalResult->ProposeNonIonizingEnergyDeposit(edep);
269 result->Clear();
270
271 return theTotalResult;
272}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4HadronicDeprecate
#define G4HadronicDeprecate(name)
Definition: G4HadronicDeprecate.hh:33
fHadronElastic
@ fHadronElastic
Definition: G4HadronicProcessType.hh:44
G4TrackStatus
G4TrackStatus
Definition: G4TrackStatus.hh:51
fStopAndKill
@ fStopAndKill
Definition: G4TrackStatus.hh:56
fStopButAlive
@ fStopButAlive
Definition: G4TrackStatus.hh:54
G4double
double G4double
Definition: G4Types.hh:64
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
CLHEP::Hep3Vector::rotate
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
G4CrossSectionDataStore::SampleZandA
G4Element * SampleZandA(const G4DynamicParticle *, const G4Material *, G4Nucleus &target)
Definition: G4CrossSectionDataStore.cc:212
G4DynamicParticle::SetMomentumDirection
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Element::GetName
const G4String & GetName() const
Definition: G4Element.hh:127
G4HadFinalState::GetEnergyChange
G4double GetEnergyChange() const
Definition: G4HadFinalState.cc:52
G4HadFinalState::GetLocalEnergyDeposit
G4double GetLocalEnergyDeposit() const
Definition: G4HadFinalState.cc:131
G4HadFinalState::GetNumberOfSecondaries
G4int GetNumberOfSecondaries() const
Definition: G4HadFinalState.cc:40
G4HadFinalState::GetMomentumChange
const G4ThreeVector & GetMomentumChange() const
Definition: G4HadFinalState.cc:68
G4HadFinalState::GetSecondary
G4HadSecondary * GetSecondary(size_t i)
Definition: G4HadFinalState.cc:109
G4HadSecondary::GetParticle
G4DynamicParticle * GetParticle()
Definition: G4HadSecondary.hh:39
G4HadronicInteraction::ApplyYourself
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
G4HadronicInteraction::GetModelName
const G4String & GetModelName() const
Definition: G4HadronicInteraction.hh:119
G4HadronicInteraction::SetRecoilEnergyThreshold
void SetRecoilEnergyThreshold(G4double val)
Definition: G4HadronicInteraction.hh:142
G4HadronicProcess::GetTargetNucleusPointer
G4Nucleus * GetTargetNucleusPointer()
Definition: G4HadronicProcess.hh:149
G4HadronicProcess::theTotalResult
G4ParticleChange * theTotalResult
Definition: G4HadronicProcess.hh:222
G4HadronicProcess::AddDataSet
void AddDataSet(G4VCrossSectionDataSet *aDataSet)
Definition: G4HadronicProcess.hh:117
G4HadronicProcess::GetCrossSectionDataStore
G4CrossSectionDataStore * GetCrossSectionDataStore()
Definition: G4HadronicProcess.hh:170
G4HadronicProcess::ChooseHadronicInteraction
G4HadronicInteraction * ChooseHadronicInteraction(G4double kineticEnergy, G4Material *aMaterial, G4Element *anElement)
Definition: G4HadronicProcess.hh:142
G4HadronicProcess::DumpState
void DumpState(const G4Track &, const G4String &, G4ExceptionDescription &)
Definition: G4HadronicProcess.cc:823
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:177
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4ParticleChange::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChange.cc:218
G4ParticleChange::ProposeEnergy
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4ParticleChange::Initialize
virtual void Initialize(const G4Track &)
Definition: G4ParticleChange.cc:228
G4ParticleDefinition::GetProcessManager
G4ProcessManager * GetProcessManager() const
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:115
G4ProcessManager::GetAtRestProcessVector
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4ProcessVector::size
G4int size() const
G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
Definition: G4ProductionCutsTable.hh:250
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4Step
Definition: G4Step.hh:78
G4Track::GetTrackStatus
G4TrackStatus GetTrackStatus() const
G4Track::GetWeight
G4double GetWeight() const
G4Track::SetWeight
void SetWeight(G4double aValue)
G4Track::GetPosition
const G4ThreeVector & GetPosition() const
G4Track::SetTouchableHandle
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4Track::GetGlobalTime
G4double GetGlobalTime() const
G4Track::GetMaterial
G4Material * GetMaterial() const
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4Track::GetTouchableHandle
const G4TouchableHandle & GetTouchableHandle() const
G4Track::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4Track::GetKineticEnergy
G4double GetKineticEnergy() const
G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VParticleChange::ProposeNonIonizingEnergyDeposit
void ProposeNonIonizingEnergyDeposit(G4double anEnergyPart)
G4VParticleChange::ProposeWeight
void ProposeWeight(G4double finalWeight)
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4VParticleChange::SetNumberOfSecondaries
void SetNumberOfSecondaries(G4int totSecondaries)
G4VProcess::ClearNumberOfInteractionLengthLeft
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:418
G4VProcess::verboseLevel
G4int verboseLevel
Definition: G4VProcess.hh:368
G4VProcess::GetProcessName
const G4String & GetProcessName() const
Definition: G4VProcess.hh:379
G4WHadronElasticProcess::PostStepDoIt
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Definition: G4WHadronElasticProcess.cc:99
G4WHadronElasticProcess::G4WHadronElasticProcess
G4WHadronElasticProcess(const G4String &procName="hadElastic")
Definition: G4WHadronElasticProcess.cc:60
G4WHadronElasticProcess::Description
virtual void Description() const
Definition: G4WHadronElasticProcess.cc:72
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76