Geant4  9.6.p02
 All Classes 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$
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  DumpState(track,"SampleZandA",ed);
117  ed << " PostStepDoIt failed on element selection" << G4endl;
118  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had003",
119  FatalException, ed);
120  }
121  G4HadronicInteraction* hadi = 0;
122  try
123  {
124  hadi = ChooseHadronicInteraction( kineticEnergy, material, elm );
125  }
126  catch(G4HadronicException & aE)
127  {
128  G4ExceptionDescription ed;
129  ed << "Target element "<< elm->GetName()<<" Z= "
130  << targNucleus->GetZ_asInt() << " A= "
131  << targNucleus->GetA_asInt() << G4endl;
132  DumpState(track,"ChooseHadronicInteraction",ed);
133  ed << " No HadronicInteraction found out" << G4endl;
134  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had005",
135  FatalException, ed);
136  }
137 
138  size_t idx = track.GetMaterialCutsCouple()->GetIndex();
139  G4double tcut = (*(G4ProductionCutsTable::GetProductionCutsTable()
140  ->GetEnergyCutsVector(3)))[idx];
141  hadi->SetRecoilEnergyThreshold(tcut);
142 
143  // Initialize the hadronic projectile from the track
144  // G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
145  G4HadProjectile theProj(track);
146  if(verboseLevel>1) {
147  G4cout << "G4HadronElasticProcess::PostStepDoIt for "
148  << part->GetParticleName()
149  << " in " << material->GetName()
150  << " Target Z= " << targNucleus->GetZ_asInt()
151  << " A= " << targNucleus->GetA_asInt() << G4endl;
152  }
153 
154  G4HadFinalState* result = 0;
155  try
156  {
157  result = hadi->ApplyYourself( theProj, *targNucleus);
158  }
159  catch(G4HadronicException aR)
160  {
161  G4ExceptionDescription ed;
162  ed << "Call for " << hadi->GetModelName() << G4endl;
163  ed << "Target element "<< elm->GetName()<<" Z= "
164  << targNucleus->GetZ_asInt()
165  << " A= " << targNucleus->GetA_asInt() << G4endl;
166  DumpState(track,"ApplyYourself",ed);
167  ed << " ApplyYourself failed" << G4endl;
168  G4Exception("G4HadronElasticProcess::PostStepDoIt", "had006",
169  FatalException, ed);
170  }
171 
172  // Check the result for catastrophic energy non-conservation
173  // cannot be applied because is not guranteed that recoil
174  // nucleus is created
175  // result = CheckResult(theProj, targNucleus, result);
176 
177  // directions
178  G4ThreeVector indir = track.GetMomentumDirection();
179  G4double phi = CLHEP::twopi*G4UniformRand();
180  G4ThreeVector it(0., 0., 1.);
181  G4ThreeVector outdir = result->GetMomentumChange();
182 
183  if(verboseLevel>1) {
184  G4cout << "Efin= " << result->GetEnergyChange()
185  << " de= " << result->GetLocalEnergyDeposit()
186  << " nsec= " << result->GetNumberOfSecondaries()
187  << " dir= " << outdir
188  << G4endl;
189  }
190 
191  // energies
192  G4double edep = result->GetLocalEnergyDeposit();
193  G4double efinal = result->GetEnergyChange();
194  if(efinal < 0.0) { efinal = 0.0; }
195  if(edep < 0.0) { edep = 0.0; }
196 
197  // NOTE: Very low energy scatters were causing numerical (FPE) errors
198  // in earlier releases; these limits have not been changed since.
199  if(efinal <= lowestEnergy) {
200  edep += efinal;
201  efinal = 0.0;
202  }
203 
204  // primary change
205  theTotalResult->ProposeEnergy(efinal);
206 
207  G4TrackStatus status = track.GetTrackStatus();
208  if(efinal > 0.0) {
209  outdir.rotate(phi, it);
210  outdir.rotateUz(indir);
211  theTotalResult->ProposeMomentumDirection(outdir);
212  } else {
213  if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
214  { status = fStopButAlive; }
215  else { status = fStopAndKill; }
216  theTotalResult->ProposeTrackStatus(status);
217  }
218 
219  //G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
220 
221  theTotalResult->SetNumberOfSecondaries(0);
222 
223  // recoil
224  if(result->GetNumberOfSecondaries() > 0) {
225  G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
226 
227  if(p->GetKineticEnergy() > tcut) {
228  theTotalResult->SetNumberOfSecondaries(1);
229  G4ThreeVector pdir = p->GetMomentumDirection();
230  // G4cout << "recoil " << pdir << G4endl;
232  pdir.rotate(phi, it);
233  pdir.rotateUz(indir);
234  // G4cout << "recoil rotated " << pdir << G4endl;
235  p->SetMomentumDirection(pdir);
236 
237  // in elastic scattering time and weight are not changed
238  G4Track* t = new G4Track(p, track.GetGlobalTime(),
239  track.GetPosition());
240  t->SetWeight(weight);
241  t->SetTouchableHandle(track.GetTouchableHandle());
242  theTotalResult->AddSecondary(t);
243 
244  } else {
245  edep += p->GetKineticEnergy();
246  delete p;
247  }
248  }
249  theTotalResult->ProposeLocalEnergyDeposit(edep);
250  theTotalResult->ProposeNonIonizingEnergyDeposit(edep);
251  result->Clear();
252 
253  return theTotalResult;
254 }
255 
256 void
257 G4HadronElasticProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
258 {
259  if(!isInitialised) {
260  isInitialised = true;
261  if(G4Neutron::Neutron() == &part) { lowestEnergy = 1.e-6*eV; }
262  }
263  G4HadronicProcess::PreparePhysicsTable(part);
264 }
265 
266 void G4HadronElasticProcess::SetLowestEnergy(G4double val)
267 {
268  lowestEnergy = val;
269 }
270 
271 void
272 G4HadronElasticProcess::SetLowestEnergyNeutron(G4double val)
273 {
274  lowestEnergy = val;
275  G4HadronicDeprecate("G4HadronElasticProcess::SetLowestEnergyNeutron()");
276 }
277 
Generated on Sat May 25 2013 14:34:09 for Geant4 by   doxygen 1.8.4