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