Geant4  10.02.p02
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 92396 2015-08-31 14:12:40Z 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 #include "G4HadronicInteraction.hh"
48 #include "G4VCrossSectionRatio.hh"
49 
51  : G4HadronicProcess(pName, fHadronElastic), isInitialised(false),
52  fDiffraction(0), fDiffractionRatio(0)
53 {
55  lowestEnergy = 1.*keV;
56 }
57 
59 {}
60 
61 void G4HadronElasticProcess::ProcessDescription(std::ostream& outFile) const
62 {
63 
64  outFile << "G4HadronElasticProcess handles the elastic scattering of \n"
65  << "hadrons by invoking the following hadronic model(s) and \n"
66  << "hadronic cross section(s).\n";
67 
68 }
69 
72  const G4Step& /*step*/)
73 {
75  theTotalResult->Initialize(track);
76  G4double weight = track.GetWeight();
78 
79  // For elastic scattering, _any_ result is considered an interaction
81 
82  G4double kineticEnergy = track.GetKineticEnergy();
83  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
84  const G4ParticleDefinition* part = dynParticle->GetDefinition();
85 
86  // NOTE: Very low energy scatters were causing numerical (FPE) errors
87  // in earlier releases; these limits have not been changed since.
88  if (kineticEnergy <= lowestEnergy) return theTotalResult;
89 
90  G4Material* material = track.GetMaterial();
91  G4Nucleus* targNucleus = GetTargetNucleusPointer();
92 
93  // Select element
94  G4Element* elm = 0;
95  try
96  {
97  elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,
98  *targNucleus);
99  }
100  catch(G4HadronicException & aR)
101  {
103  aR.Report(ed);
104  DumpState(track,"SampleZandA",ed);
105  ed << " PostStepDoIt failed on element selection" << G4endl;
106  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had003",
107  FatalException, ed);
108  }
109 
110  // Initialize the hadronic projectile from the track
111  // G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
112  G4HadProjectile theProj(track);
113  G4HadronicInteraction* hadi = 0;
114  G4HadFinalState* result = 0;
115 
116  if(fDiffraction)
117  {
118  G4double ratio =
119  fDiffractionRatio->ComputeRatio(part, kineticEnergy,
120  targNucleus->GetZ_asInt(),
121  targNucleus->GetA_asInt());
122  // diffraction is chosen
123  if(ratio > 0.0 && G4UniformRand() < ratio)
124  {
125  try
126  {
127 // if(part->GetParticleName() == "pi-")
128 // G4cout<<part->GetParticleName()<<"; "<<kineticEnergy/CLHEP::GeV<<" GeV; r = "<<ratio<<G4endl;
129 
130  result = fDiffraction->ApplyYourself(theProj, *targNucleus);
131  }
132  catch(G4HadronicException aR)
133  {
135  aR.Report(ed);
136  ed << "Call for " << fDiffraction->GetModelName() << G4endl;
137  ed << "Target element "<< elm->GetName()<<" Z= "
138  << targNucleus->GetZ_asInt()
139  << " A= " << targNucleus->GetA_asInt() << G4endl;
140  DumpState(track,"ApplyYourself",ed);
141  ed << " ApplyYourself failed" << G4endl;
142  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had006",
143  FatalException, ed);
144  }
145  // Check the result for catastrophic energy non-conservation
146  result = CheckResult(theProj, *targNucleus, result);
147 
148  result->SetTrafoToLab(theProj.GetTrafoToLab());
150 
151  FillResult(result, track);
152 
153  if (epReportLevel != 0) {
154  CheckEnergyMomentumConservation(track, *targNucleus);
155  }
156  return theTotalResult;
157  }
158  }
159 
160  // ordinary elastic scattering
161  try
162  {
163  hadi = ChooseHadronicInteraction( theProj, *targNucleus, material, elm );
164  }
165  catch(G4HadronicException & aE)
166  {
168  aE.Report(ed);
169  ed << "Target element "<< elm->GetName()<<" Z= "
170  << targNucleus->GetZ_asInt() << " A= "
171  << targNucleus->GetA_asInt() << G4endl;
172  DumpState(track,"ChooseHadronicInteraction",ed);
173  ed << " No HadronicInteraction found out" << G4endl;
174  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had005",
175  FatalException, ed);
176  }
177 
178  size_t idx = track.GetMaterialCutsCouple()->GetIndex();
180  ->GetEnergyCutsVector(3)))[idx];
181  hadi->SetRecoilEnergyThreshold(tcut);
182 
183  if(verboseLevel>1) {
184  G4cout << "G4HadronElasticProcess::PostStepDoIt for "
185  << part->GetParticleName()
186  << " in " << material->GetName()
187  << " Target Z= " << targNucleus->GetZ_asInt()
188  << " A= " << targNucleus->GetA_asInt() << G4endl;
189  }
190 
191  try
192  {
193  result = hadi->ApplyYourself( theProj, *targNucleus);
194  }
195  catch(G4HadronicException aR)
196  {
198  aR.Report(ed);
199  ed << "Call for " << hadi->GetModelName() << G4endl;
200  ed << "Target element "<< elm->GetName()<<" Z= "
201  << targNucleus->GetZ_asInt()
202  << " A= " << targNucleus->GetA_asInt() << G4endl;
203  DumpState(track,"ApplyYourself",ed);
204  ed << " ApplyYourself failed" << G4endl;
205  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had006",
206  FatalException, ed);
207  }
208 
209  // Check the result for catastrophic energy non-conservation
210  // cannot be applied because is not guranteed that recoil
211  // nucleus is created
212  // result = CheckResult(theProj, targNucleus, result);
213 
214  // directions
215  G4ThreeVector indir = track.GetMomentumDirection();
217  G4ThreeVector it(0., 0., 1.);
218  G4ThreeVector outdir = result->GetMomentumChange();
219 
220  if(verboseLevel>1) {
221  G4cout << "Efin= " << result->GetEnergyChange()
222  << " de= " << result->GetLocalEnergyDeposit()
223  << " nsec= " << result->GetNumberOfSecondaries()
224  << " dir= " << outdir
225  << G4endl;
226  }
227 
228  // energies
229  G4double edep = result->GetLocalEnergyDeposit();
230  G4double efinal = result->GetEnergyChange();
231  if(efinal < 0.0) { efinal = 0.0; }
232  if(edep < 0.0) { edep = 0.0; }
233 
234  // NOTE: Very low energy scatters were causing numerical (FPE) errors
235  // in earlier releases; these limits have not been changed since.
236  if(efinal <= lowestEnergy) {
237  edep += efinal;
238  efinal = 0.0;
239  }
240 
241  // primary change
242  theTotalResult->ProposeEnergy(efinal);
243 
244  G4TrackStatus status = track.GetTrackStatus();
245  if(efinal > 0.0) {
246  outdir.rotate(phi, it);
247  outdir.rotateUz(indir);
249  } else {
250  if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
251  { status = fStopButAlive; }
252  else { status = fStopAndKill; }
254  }
255 
256  //G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
257 
259 
260  // recoil
261  if(result->GetNumberOfSecondaries() > 0) {
262  G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
263 
264  if(p->GetKineticEnergy() > tcut) {
267  // G4cout << "recoil " << pdir << G4endl;
269  pdir.rotate(phi, it);
270  pdir.rotateUz(indir);
271  // G4cout << "recoil rotated " << pdir << G4endl;
272  p->SetMomentumDirection(pdir);
273 
274  // in elastic scattering time and weight are not changed
275  G4Track* t = new G4Track(p, track.GetGlobalTime(),
276  track.GetPosition());
277  t->SetWeight(weight);
280 
281  } else {
282  edep += p->GetKineticEnergy();
283  delete p;
284  }
285  }
288  result->Clear();
289 
290  return theTotalResult;
291 }
292 
293 void
295 {
296  if(!isInitialised) {
297  isInitialised = true;
298  if(G4Neutron::Neutron() == &part) { lowestEnergy = 1.e-6*eV; }
299  }
301 }
302 
304 {
305  lowestEnergy = val;
306 }
307 
308 void
310 {
311  lowestEnergy = val;
312  G4HadronicDeprecate("G4HadronElasticProcess::SetLowestEnergyNeutron()");
313 }
314 
317 {
318  if(hi && xsr) {
319  fDiffraction = hi;
320  fDiffractionRatio = xsr;
321  }
322 }
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
virtual G4double ComputeRatio(const G4ParticleDefinition *, G4double kinEnergy, G4int, G4int)=0
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
CLHEP::Hep3Vector G4ThreeVector
G4LorentzRotation & GetTrafoToLab()
const G4DynamicParticle * GetDynamicParticle() const
const G4String & GetName() const
Definition: G4Material.hh:178
const G4ThreeVector & GetMomentumChange() const
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
void CheckEnergyMomentumConservation(const G4Track &, const G4Nucleus &)
const G4String & GetModelName() const
void SetDiffraction(G4HadronicInteraction *, G4VCrossSectionRatio *)
#define G4HadronicDeprecate(name)
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:447
G4ProcessManager * GetProcessManager() const
G4VCrossSectionRatio * fDiffractionRatio
const G4String & GetParticleName() const
void SetWeight(G4double aValue)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void ProposeWeight(G4double finalWeight)
void DumpState(const G4Track &, const G4String &, G4ExceptionDescription &)
void SetTrafoToLab(const G4LorentzRotation &aT)
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4double GetKineticEnergy() const
void AddDataSet(G4VCrossSectionDataSet *aDataSet)
void FillResult(G4HadFinalState *aR, const G4Track &aT)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4ParticleChange * theTotalResult
void ProposeNonIonizingEnergyDeposit(G4double anEnergyPart)
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
Definition: G4SIunits.hh:75
G4CrossSectionDataStore * GetCrossSectionDataStore()
Definition: G4Step.hh:76
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetGlobalTime() const
const G4TouchableHandle & GetTouchableHandle() const
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4int size() const
virtual void Initialize(const G4Track &)
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:212
const G4ThreeVector & GetMomentumDirection() const
G4Nucleus * GetTargetNucleusPointer()
void SetNumberOfSecondaries(G4int totSecondaries)
G4DynamicParticle * GetParticle()
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void ProposeEnergy(G4double finalEnergy)
G4HadronicInteraction * fDiffraction
void AddSecondary(G4Track *aSecondary)
virtual void ProcessDescription(std::ostream &outFile) const
G4double GetWeight() const
G4HadronicInteraction * ChooseHadronicInteraction(const G4HadProjectile &aHadProjectile, G4Nucleus &aTargetNucleus, G4Material *aMaterial, G4Element *anElement)
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:213
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
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)
G4HadFinalState * CheckResult(const G4HadProjectile &thePro, const G4Nucleus &targetNucleus, G4HadFinalState *result)