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