Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
39 #include "G4SystemOfUnits.hh"
40 #include "G4Nucleus.hh"
41 #include "G4ProcessManager.hh"
44 #include "G4ProductionCutsTable.hh"
45 #include "G4HadronicException.hh"
46 #include "G4HadronicDeprecate.hh"
47 
49  : G4HadronicProcess(pName, fHadronElastic), isInitialised(false)
50 {
52  lowestEnergy = 1.*keV;
53 }
54 
56 {}
57 
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 
85  const G4Step& /*step*/)
86 {
88  theTotalResult->Initialize(track);
89  G4double weight = track.GetWeight();
91 
92  // For elastic scattering, _any_ result is considered an interaction
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  {
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  {
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();
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  {
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 
211  if(efinal > 0.0) {
212  outdir.rotate(phi, it);
213  outdir.rotateUz(indir);
215  } else {
216  if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
217  { status = fStopButAlive; }
218  else { status = fStopAndKill; }
220  }
221 
222  //G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
223 
225 
226  // recoil
227  if(result->GetNumberOfSecondaries() > 0) {
228  G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
229 
230  if(p->GetKineticEnergy() > tcut) {
233  // G4cout << "recoil " << pdir << G4endl;
234  //!! is not needed for models inheriting G4HadronElastic
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);
246 
247  } else {
248  edep += p->GetKineticEnergy();
249  delete p;
250  }
251  }
254  result->Clear();
255 
256  return theTotalResult;
257 }
258 
259 void
261 {
262  if(!isInitialised) {
263  isInitialised = true;
264  if(G4Neutron::Neutron() == &part) { lowestEnergy = 1.e-6*eV; }
265  }
267 }
268 
270 {
271  lowestEnergy = val;
272 }
273 
274 void
276 {
277  lowestEnergy = val;
278  G4HadronicDeprecate("G4HadronElasticProcess::SetLowestEnergyNeutron()");
279 }
280 
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4HadronElasticProcess(const G4String &procName="hadElastic")
G4int verboseLevel
Definition: G4VProcess.hh:368
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4HadSecondary * GetSecondary(size_t i)
G4double GetKineticEnergy() const
std::ofstream outFile
Definition: GammaRayTel.cc:68
const G4DynamicParticle * GetDynamicParticle() const
const char * p
Definition: xmltok.h:285
const G4String & GetName() const
Definition: G4Material.hh:176
const G4ThreeVector & GetPosition() const
G4TrackStatus GetTrackStatus() const
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
void SetTouchableHandle(const G4TouchableHandle &apValue)
virtual void SetLowestEnergy(G4double)
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4double GetEnergyChange() const
G4ParticleDefinition * GetDefinition() const
const G4String & GetModelName() const
#define G4HadronicDeprecate(name)
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:447
G4ProcessManager * GetProcessManager() const
const G4String & GetParticleName() const
void SetWeight(G4double aValue)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
string material
Definition: eplot.py:19
void ProposeWeight(G4double finalWeight)
void DumpState(const G4Track &, const G4String &, G4ExceptionDescription &)
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4double GetKineticEnergy() const
void AddDataSet(G4VCrossSectionDataSet *aDataSet)
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4ParticleChange * theTotalResult
void ProposeNonIonizingEnergyDeposit(G4double anEnergyPart)
const G4ThreeVector & GetMomentumDirection() const
G4HadronicInteraction * ChooseHadronicInteraction(G4double kineticEnergy, G4Material *aMaterial, G4Element *anElement)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
G4CrossSectionDataStore * GetCrossSectionDataStore()
Definition: G4Step.hh:76
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetGlobalTime() const
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
const G4TouchableHandle & GetTouchableHandle() const
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
int status
Definition: tracer.cxx:24
G4int size() const
virtual void Initialize(const G4Track &)
static G4ProductionCutsTable * GetProductionCutsTable()
const G4ThreeVector & GetMomentumDirection() const
G4Nucleus * GetTargetNucleusPointer()
const G4ThreeVector & GetMomentumChange() const
void SetNumberOfSecondaries(G4int totSecondaries)
G4DynamicParticle * GetParticle()
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void ProposeEnergy(G4double finalEnergy)
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
virtual void Description() const
#define G4endl
Definition: G4ios.hh:61
void SetRecoilEnergyThreshold(G4double val)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4Element * SampleZandA(const G4DynamicParticle *, const G4Material *, G4Nucleus &target)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
const G4String & GetName() const
Definition: G4Element.hh:127
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
G4TrackStatus
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
void Report(std::ostream &aS)
G4int GetNumberOfSecondaries() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
G4double GetLocalEnergyDeposit() const
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
virtual void SetLowestEnergyNeutron(G4double)