Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4WHadronElasticProcess.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 21 April 2006 V.Ivanchenko
31 //
32 // Modified:
33 // 24.04.06 V.Ivanchenko add neutron scattering on hydrogen from CHIPS
34 // 07.06.06 V.Ivanchenko fix problem of rotation of final state
35 // 25.07.06 V.Ivanchenko add 19 MeV low energy for CHIPS
36 // 26.09.06 V.Ivanchenko add lowestEnergy
37 // 20.10.06 V.Ivanchenko initialise lowestEnergy=0 for neitrals, eV for charged
38 // 23.01.07 V.Ivanchenko add cross section interfaces with Z and A
39 // 02.05.07 V.Ivanchenko add He3
40 // 13.01.10 M.Kosov: Commented not used G4QElasticCrossSection & G4QCHIPSWorld
41 // 24.02.11 V.Ivanchenko use particle name in IfApplicable,
42 // added anti particles for light ions
43 // 07.09.11 M.Kelsey: Follow chanhe to G4HadFinalState interface
44 // 14.09.12 M.Kelsey: Pass subType code to base ctor
45 //
46 
47 #include <iostream>
48 #include <typeinfo>
49 
50 #include "G4WHadronElasticProcess.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4Nucleus.hh"
53 #include "G4ProcessManager.hh"
54 #include "G4CrossSectionDataStore.hh"
55 #include "G4HadronElasticDataSet.hh"
56 #include "G4ProductionCutsTable.hh"
57 #include "G4HadronicException.hh"
58 #include "G4HadronicDeprecate.hh"
59 
60 G4WHadronElasticProcess::G4WHadronElasticProcess(const G4String& pName)
61  : G4HadronicProcess(pName, fHadronElastic) {
62  AddDataSet(new G4HadronElasticDataSet);
63  theNeutron = G4Neutron::Neutron();
64  lowestEnergy = 1.*keV;
65  lowestEnergyNeutron = 1.e-6*eV;
66  G4HadronicDeprecate("G4WHadronElasticProcess");
67 }
68 
69 G4WHadronElasticProcess::~G4WHadronElasticProcess()
70 {}
71 
72 void G4WHadronElasticProcess::Description() const
73 {
74  char* dirName = getenv("G4PhysListDocDir");
75  if (dirName) {
76  std::ofstream outFile;
77  G4String outFileName = GetProcessName() + ".html";
78  G4String pathName = G4String(dirName) + "/" + outFileName;
79  outFile.open(pathName);
80  outFile << "<html>\n";
81  outFile << "<head>\n";
82 
83  outFile << "<title>Description of G4WHadronElasticProcess</title>\n";
84  outFile << "</head>\n";
85  outFile << "<body>\n";
86 
87  outFile << "G4WHadronElasticProcess handles the elastic scattering of\n"
88  << "hadrons by invoking one or more hadronic models and one or\n"
89  << "more hadronic cross sections.\n";
90 
91  outFile << "</body>\n";
92  outFile << "</html>\n";
93  outFile.close();
94  }
95 }
96 
97 
98 G4VParticleChange*
99 G4WHadronElasticProcess::PostStepDoIt(const G4Track& track,
100  const G4Step& /*step*/)
101 {
102  theTotalResult->Clear();
103  theTotalResult->Initialize(track);
104  G4double weight = track.GetWeight();
105  theTotalResult->ProposeWeight(weight);
106 
107  // For elastic scattering, _any_ result is considered an interaction
108  ClearNumberOfInteractionLengthLeft();
109 
110  G4double kineticEnergy = track.GetKineticEnergy();
111  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
112  const G4ParticleDefinition* part = dynParticle->GetDefinition();
113 
114  // NOTE: Very low energy scatters were causing numerical (FPE) errors
115  // in earlier releases; these limits have not been changed since.
116  if (part == theNeutron) {
117  if(kineticEnergy <= lowestEnergyNeutron) return theTotalResult;
118  } else if(kineticEnergy <= lowestEnergy) return theTotalResult;
119 
120  G4Material* material = track.GetMaterial();
121  G4Nucleus* targNucleus = GetTargetNucleusPointer();
122 
123  // Select element
124  G4Element* elm = 0;
125  try
126  {
127  elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,
128  *targNucleus);
129  }
130  catch(G4HadronicException & aR)
131  {
132  G4ExceptionDescription ed;
133  DumpState(track,"SampleZandA",ed);
134  ed << " PostStepDoIt failed on element selection" << G4endl;
135  G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had003",
136  FatalException, ed);
137  }
138  G4HadronicInteraction* hadi = 0;
139  try
140  {
141  hadi = ChooseHadronicInteraction( kineticEnergy, material, elm );
142  }
143  catch(G4HadronicException & aE)
144  {
145  G4ExceptionDescription ed;
146  ed << "Target element "<< elm->GetName()<<" Z= "
147  << targNucleus->GetZ_asInt() << " A= "
148  << targNucleus->GetA_asInt() << G4endl;
149  DumpState(track,"ChooseHadronicInteraction",ed);
150  ed << " No HadronicInteraction found out" << G4endl;
151  G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had005",
152  FatalException, ed);
153  }
154 
155  size_t idx = track.GetMaterialCutsCouple()->GetIndex();
156  G4double tcut = (*(G4ProductionCutsTable::GetProductionCutsTable()
157  ->GetEnergyCutsVector(3)))[idx];
158  hadi->SetRecoilEnergyThreshold(tcut);
159 
160  // Initialize the hadronic projectile from the track
161  // G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
162  G4HadProjectile theProj(track);
163  if(verboseLevel>1) {
164  G4cout << "G4WHadronElasticProcess::PostStepDoIt for "
165  << part->GetParticleName()
166  << " in " << material->GetName()
167  << " Target Z= " << targNucleus->GetZ_asInt()
168  << " A= " << targNucleus->GetA_asInt() << G4endl;
169  }
170 
171  G4HadFinalState* result = 0;
172  try
173  {
174  result = hadi->ApplyYourself( theProj, *targNucleus);
175  }
176  catch(G4HadronicException aR)
177  {
178  G4ExceptionDescription ed;
179  ed << "Call for " << hadi->GetModelName() << G4endl;
180  ed << "Target element "<< elm->GetName()<<" Z= "
181  << targNucleus->GetZ_asInt()
182  << " A= " << targNucleus->GetA_asInt() << G4endl;
183  DumpState(track,"ApplyYourself",ed);
184  ed << " ApplyYourself failed" << G4endl;
185  G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had006",
186  FatalException, ed);
187  }
188 
189  // Check the result for catastrophic energy non-conservation
190  // cannot be applied because is not guranteed that recoil
191  // nucleus is created
192  // result = CheckResult(theProj, targNucleus, result);
193 
194  // directions
195  G4ThreeVector indir = track.GetMomentumDirection();
196  G4double phi = CLHEP::twopi*G4UniformRand();
197  G4ThreeVector it(0., 0., 1.);
198  G4ThreeVector outdir = result->GetMomentumChange();
199 
200  if(verboseLevel>1) {
201  G4cout << "Efin= " << result->GetEnergyChange()
202  << " de= " << result->GetLocalEnergyDeposit()
203  << " nsec= " << result->GetNumberOfSecondaries()
204  << " dir= " << outdir
205  << G4endl;
206  }
207 
208  // energies
209  G4double edep = result->GetLocalEnergyDeposit();
210  G4double efinal = result->GetEnergyChange();
211  if(efinal < 0.0) { efinal = 0.0; }
212  if(edep < 0.0) { edep = 0.0; }
213 
214  // NOTE: Very low energy scatters were causing numerical (FPE) errors
215  // in earlier releases; these limits have not been changed since.
216  if((part == theNeutron && efinal <= lowestEnergyNeutron) ||
217  (part != theNeutron && efinal <= lowestEnergy)) {
218  edep += efinal;
219  efinal = 0.0;
220  }
221 
222  // primary change
223  theTotalResult->ProposeEnergy(efinal);
224 
225  G4TrackStatus status = track.GetTrackStatus();
226  if(efinal > 0.0) {
227  outdir.rotate(phi, it);
228  outdir.rotateUz(indir);
229  theTotalResult->ProposeMomentumDirection(outdir);
230  } else {
231  if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
232  { status = fStopButAlive; }
233  else { status = fStopAndKill; }
234  theTotalResult->ProposeTrackStatus(status);
235  }
236 
237  //G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
238 
239  theTotalResult->SetNumberOfSecondaries(0);
240 
241  // recoil
242  if(result->GetNumberOfSecondaries() > 0) {
243  G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
244 
245  if(p->GetKineticEnergy() > lowestEnergy) {
246  theTotalResult->SetNumberOfSecondaries(1);
247  G4ThreeVector pdir = p->GetMomentumDirection();
248  // G4cout << "recoil " << pdir << G4endl;
250  pdir.rotate(phi, it);
251  pdir.rotateUz(indir);
252  // G4cout << "recoil rotated " << pdir << G4endl;
253  p->SetMomentumDirection(pdir);
254 
255  // in elastic scattering time and weight are not changed
256  G4Track* t = new G4Track(p, track.GetGlobalTime(),
257  track.GetPosition());
258  t->SetWeight(weight);
259  t->SetTouchableHandle(track.GetTouchableHandle());
260  theTotalResult->AddSecondary(t);
261 
262  } else {
263  edep += p->GetKineticEnergy();
264  delete p;
265  }
266  }
267  theTotalResult->ProposeLocalEnergyDeposit(edep);
268  theTotalResult->ProposeNonIonizingEnergyDeposit(edep);
269  result->Clear();
270 
271  return theTotalResult;
272 }
Generated on Sat May 25 2013 14:34:09 for Geant4 by   doxygen 1.8.4