Geant4  10.01.p03
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 86393 2014-11-10 16:34:27Z 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 
62 {
63  char* dirName = getenv("G4PhysListDocDir");
64  if (dirName) {
65  std::ofstream outFile;
66  G4String outFileName = GetProcessName() + ".html";
67  G4String pathName = G4String(dirName) + "/" + outFileName;
68  outFile.open(pathName);
69  outFile << "<html>\n";
70  outFile << "<head>\n";
71 
72  outFile << "<title>Description of G4HadronElasticProcess</title>\n";
73  outFile << "</head>\n";
74  outFile << "<body>\n";
75 
76  outFile << "G4HadronElasticProcess handles the elastic scattering of\n"
77  << "hadrons by invoking one or more hadronic models and one or\n"
78  << "more hadronic cross sections.\n";
79 
80  outFile << "</body>\n";
81  outFile << "</html>\n";
82  outFile.close();
83  }
84 }
85 
88  const G4Step& /*step*/)
89 {
91  theTotalResult->Initialize(track);
92  G4double weight = track.GetWeight();
94 
95  // For elastic scattering, _any_ result is considered an interaction
97 
98  G4double kineticEnergy = track.GetKineticEnergy();
99  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
100  const G4ParticleDefinition* part = dynParticle->GetDefinition();
101 
102  // NOTE: Very low energy scatters were causing numerical (FPE) errors
103  // in earlier releases; these limits have not been changed since.
104  if (kineticEnergy <= lowestEnergy) return theTotalResult;
105 
106  G4Material* material = track.GetMaterial();
107  G4Nucleus* targNucleus = GetTargetNucleusPointer();
108 
109  // Select element
110  G4Element* elm = 0;
111  try
112  {
113  elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,
114  *targNucleus);
115  }
116  catch(G4HadronicException & aR)
117  {
119  aR.Report(ed);
120  DumpState(track,"SampleZandA",ed);
121  ed << " PostStepDoIt failed on element selection" << G4endl;
122  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had003",
123  FatalException, ed);
124  }
125 
126  // Initialize the hadronic projectile from the track
127  // G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
128  G4HadProjectile theProj(track);
129  G4HadronicInteraction* hadi = 0;
130  G4HadFinalState* result = 0;
131 
132  if(fDiffraction) {
133  G4double ratio =
134  fDiffractionRatio->ComputeRatio(part, kineticEnergy,
135  targNucleus->GetZ_asInt(),
136  targNucleus->GetA_asInt());
137  // diffraction is chosen
138  if(ratio > 0.0 && G4UniformRand() < ratio) {
139  try
140  {
141  result = fDiffraction->ApplyYourself(theProj, *targNucleus);
142  }
143  catch(G4HadronicException aR)
144  {
146  aR.Report(ed);
147  ed << "Call for " << fDiffraction->GetModelName() << G4endl;
148  ed << "Target element "<< elm->GetName()<<" Z= "
149  << targNucleus->GetZ_asInt()
150  << " A= " << targNucleus->GetA_asInt() << G4endl;
151  DumpState(track,"ApplyYourself",ed);
152  ed << " ApplyYourself failed" << G4endl;
153  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had006",
154  FatalException, ed);
155  }
156  // Check the result for catastrophic energy non-conservation
157  result = CheckResult(theProj, *targNucleus, result);
158 
159  result->SetTrafoToLab(theProj.GetTrafoToLab());
161 
162  FillResult(result, track);
163 
164  if (epReportLevel != 0) {
165  CheckEnergyMomentumConservation(track, *targNucleus);
166  }
167  return theTotalResult;
168  }
169  }
170 
171  // ordinary elastic scattering
172  try
173  {
174  hadi = ChooseHadronicInteraction( theProj, *targNucleus, material, elm );
175  }
176  catch(G4HadronicException & aE)
177  {
179  aE.Report(ed);
180  ed << "Target element "<< elm->GetName()<<" Z= "
181  << targNucleus->GetZ_asInt() << " A= "
182  << targNucleus->GetA_asInt() << G4endl;
183  DumpState(track,"ChooseHadronicInteraction",ed);
184  ed << " No HadronicInteraction found out" << G4endl;
185  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had005",
186  FatalException, ed);
187  }
188 
189  size_t idx = track.GetMaterialCutsCouple()->GetIndex();
191  ->GetEnergyCutsVector(3)))[idx];
192  hadi->SetRecoilEnergyThreshold(tcut);
193 
194  if(verboseLevel>1) {
195  G4cout << "G4HadronElasticProcess::PostStepDoIt for "
196  << part->GetParticleName()
197  << " in " << material->GetName()
198  << " Target Z= " << targNucleus->GetZ_asInt()
199  << " A= " << targNucleus->GetA_asInt() << G4endl;
200  }
201 
202  try
203  {
204  result = hadi->ApplyYourself( theProj, *targNucleus);
205  }
206  catch(G4HadronicException aR)
207  {
209  aR.Report(ed);
210  ed << "Call for " << hadi->GetModelName() << G4endl;
211  ed << "Target element "<< elm->GetName()<<" Z= "
212  << targNucleus->GetZ_asInt()
213  << " A= " << targNucleus->GetA_asInt() << G4endl;
214  DumpState(track,"ApplyYourself",ed);
215  ed << " ApplyYourself failed" << G4endl;
216  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had006",
217  FatalException, ed);
218  }
219 
220  // Check the result for catastrophic energy non-conservation
221  // cannot be applied because is not guranteed that recoil
222  // nucleus is created
223  // result = CheckResult(theProj, targNucleus, result);
224 
225  // directions
226  G4ThreeVector indir = track.GetMomentumDirection();
227  G4double phi = CLHEP::twopi*G4UniformRand();
228  G4ThreeVector it(0., 0., 1.);
229  G4ThreeVector outdir = result->GetMomentumChange();
230 
231  if(verboseLevel>1) {
232  G4cout << "Efin= " << result->GetEnergyChange()
233  << " de= " << result->GetLocalEnergyDeposit()
234  << " nsec= " << result->GetNumberOfSecondaries()
235  << " dir= " << outdir
236  << G4endl;
237  }
238 
239  // energies
240  G4double edep = result->GetLocalEnergyDeposit();
241  G4double efinal = result->GetEnergyChange();
242  if(efinal < 0.0) { efinal = 0.0; }
243  if(edep < 0.0) { edep = 0.0; }
244 
245  // NOTE: Very low energy scatters were causing numerical (FPE) errors
246  // in earlier releases; these limits have not been changed since.
247  if(efinal <= lowestEnergy) {
248  edep += efinal;
249  efinal = 0.0;
250  }
251 
252  // primary change
253  theTotalResult->ProposeEnergy(efinal);
254 
255  G4TrackStatus status = track.GetTrackStatus();
256  if(efinal > 0.0) {
257  outdir.rotate(phi, it);
258  outdir.rotateUz(indir);
260  } else {
261  if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
262  { status = fStopButAlive; }
263  else { status = fStopAndKill; }
265  }
266 
267  //G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
268 
270 
271  // recoil
272  if(result->GetNumberOfSecondaries() > 0) {
273  G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
274 
275  if(p->GetKineticEnergy() > tcut) {
278  // G4cout << "recoil " << pdir << G4endl;
280  pdir.rotate(phi, it);
281  pdir.rotateUz(indir);
282  // G4cout << "recoil rotated " << pdir << G4endl;
283  p->SetMomentumDirection(pdir);
284 
285  // in elastic scattering time and weight are not changed
286  G4Track* t = new G4Track(p, track.GetGlobalTime(),
287  track.GetPosition());
288  t->SetWeight(weight);
291 
292  } else {
293  edep += p->GetKineticEnergy();
294  delete p;
295  }
296  }
299  result->Clear();
300 
301  return theTotalResult;
302 }
303 
304 void
306 {
307  if(!isInitialised) {
308  isInitialised = true;
309  if(G4Neutron::Neutron() == &part) { lowestEnergy = 1.e-6*eV; }
310  }
312 }
313 
315 {
316  lowestEnergy = val;
317 }
318 
319 void
321 {
322  lowestEnergy = val;
323  G4HadronicDeprecate("G4HadronElasticProcess::SetLowestEnergyNeutron()");
324 }
325 
328 {
329  if(hi && xsr) {
330  fDiffraction = hi;
331  fDiffractionRatio = xsr;
332  }
333 }
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:93
G4GLOB_DLL std::ostream G4cout
G4ParticleChange * theTotalResult
void ProposeNonIonizingEnergyDeposit(G4double anEnergyPart)
const G4ThreeVector & GetMomentumDirection() const
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
G4int size() const
virtual void Initialize(const G4Track &)
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:194
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)
G4double GetWeight() const
virtual void Description() 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:195
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)