Geant4  10.00.p02
G4HadronStoppingProcess.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: G4HadronStoppingProcess.cc 66367 2012-12-18 09:18:08Z gcosmo $
27 //
28 //---------------------------------------------------------------------
29 //
30 // GEANT4 Class
31 //
32 // File name: G4HadronStoppingProcess
33 //
34 // Author V.Ivanchenko 21 April 2012
35 //
36 //
37 // Class Description:
38 //
39 // Base process class for nuclear capture of negatively charged particles
40 //
41 // Modifications:
42 //
43 // 20120522 M. Kelsey -- Set enableAtRestDoIt flag for G4ProcessManager
44 // 20120914 M. Kelsey -- Pass subType in base ctor, remove enable flags
45 // 20121004 K. Genser -- use G4HadronicProcessType in the constructor
46 // 20121016 K. Genser -- Reverting to use one argument c'tor
47 //
48 //------------------------------------------------------------------------
49 
52 #include "G4HadronicProcessType.hh"
53 #include "G4EmCaptureCascade.hh"
54 #include "G4Nucleus.hh"
55 #include "G4HadFinalState.hh"
56 #include "G4HadProjectile.hh"
57 #include "G4HadSecondary.hh"
58 #include "G4Material.hh"
59 #include "G4Element.hh"
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 
65 {
66  // Modify G4VProcess flags to emulate G4VRest instead of G4VDiscrete
67  enableAtRestDoIt = true;
68  enablePostStepDoIt = false;
69 
71  fEmCascade = new G4EmCaptureCascade(); // Owned by InteractionRegistry
72  fBoundDecay = 0;
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79 {
81  delete fElementSelector;
82  // NOTE: fEmCascade and fEmBoundDecay owned by registry, not locally
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86 
88 {
89  return (p.GetPDGCharge() < 0.0);
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93 
95 {
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
102 {
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 
110 {
111  *condition = NotForced;
112  return 0.0;
113 }
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 
119 {
120  *condition = NotForced;
121  return DBL_MAX;
122 }
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125 
127  const G4Step&)
128 {
129  // if primary is not Alive then do nothing
130  theTotalResult->Initialize(track);
131 
132  G4Nucleus* nucleus = GetTargetNucleusPointer();
133  G4Element* elm = fElementSelector->SelectZandA(track, nucleus);
134 
135  G4HadFinalState* result = 0;
136  thePro.Initialise(track);
137  G4double time0 = track.GetGlobalTime();
138  G4bool nuclearCapture = true;
139 
140  // Do the electromagnetic cascade in the nuclear field.
141  // EM cascade should keep G4HadFinalState object,
142  // because it will not be deleted at the end of this method
143  //
144  result = fEmCascade->ApplyYourself(thePro, *nucleus);
145  G4double ebound = result->GetLocalEnergyDeposit();
146  G4double edep = 0.0;
147  G4int nSecondaries = result->GetNumberOfSecondaries();
148 
149  // Try decay from bound level
150  // For mu- the time of projectile should be changed.
151  // Decay should keep G4HadFinalState object,
152  // because it will not be deleted at the end of this method
153  //
154  thePro.SetBoundEnergy(ebound);
155  if(fBoundDecay) {
156  G4HadFinalState* resultDecay =
157  fBoundDecay->ApplyYourself(thePro, *nucleus);
158  G4int n = resultDecay->GetNumberOfSecondaries();
159  if(0 < n) {
160  nSecondaries += n;
161  result->AddSecondaries(resultDecay);
162  }
163  if(resultDecay->GetStatusChange() == stopAndKill) {
164  nuclearCapture = false;
165  }
166  resultDecay->Clear();
167  }
168 
169  if(nuclearCapture) {
170 
171  // select model
172  G4HadronicInteraction* model = 0;
173  try {
174  model = ChooseHadronicInteraction(0.0, track.GetMaterial(), elm);
175  }
176  catch(G4HadronicException & aE) {
178  ed << "Target element "<<elm->GetName()<<" Z= "
179  << nucleus->GetZ_asInt() << " A= "
180  << nucleus->GetA_asInt() << G4endl;
181  DumpState(track,"ChooseHadronicInteraction",ed);
182  ed << " No HadronicInteraction found out" << G4endl;
183  G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had005",
184  FatalException, ed);
185  }
186 
187  G4HadFinalState* resultNuc = 0;
188  G4int reentryCount = 0;
189  do {
190  // sample final state
191  // nuclear interaction should keep G4HadFinalState object
192  // model should define time of each secondary particle
193  try {
194  resultNuc = model->ApplyYourself(thePro, *nucleus);
195  ++reentryCount;
196  }
197  catch(G4HadronicException aR) {
199  ed << "Call for " << model->GetModelName() << G4endl;
200  ed << "Target element "<<elm->GetName()<<" Z= "
201  << nucleus->GetZ_asInt()
202  << " A= " << nucleus->GetA_asInt() << G4endl;
203  DumpState(track,"ApplyYourself",ed);
204  ed << " ApplyYourself failed" << G4endl;
205  G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
206  FatalException, ed);
207  }
208 
209  // Check the result for catastrophic energy non-conservation
210  resultNuc = CheckResult(thePro, *nucleus, resultNuc);
211 
212  if(reentryCount>100) {
214  ed << "Call for " << model->GetModelName() << G4endl;
215  ed << "Target element "<<elm->GetName()<<" Z= "
216  << nucleus->GetZ_asInt()
217  << " A= " << nucleus->GetA_asInt() << G4endl;
218  DumpState(track,"ApplyYourself",ed);
219  ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
220  G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
221  FatalException, ed);
222  }
223  }
224  while(!resultNuc);
225 
226  edep = resultNuc->GetLocalEnergyDeposit();
227  nSecondaries += resultNuc->GetNumberOfSecondaries();
228  result->AddSecondaries(resultNuc);
229  resultNuc->Clear();
230  }
231 
232  // Fill results
233  //
236  theTotalResult->SetNumberOfSecondaries(nSecondaries);
237  G4double w = track.GetWeight();
239  for(G4int i=0; i<nSecondaries; ++i) {
240  G4HadSecondary* sec = result->GetSecondary(i);
241 
242  // add track global time to the reaction time
243  G4double time = sec->GetTime();
244  if(time < 0.0) { time = 0.0; }
245  time += time0;
246 
247  // create secondary track
248  G4Track* t = new G4Track(sec->GetParticle(),
249  time,
250  track.GetPosition());
251  t->SetWeight(w*sec->GetWeight());
254  }
255  result->Clear();
256 
257  if (epReportLevel != 0) {
258  CheckEnergyMomentumConservation(track, *nucleus);
259  }
260  return theTotalResult;
261 }
262 
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
264 
266 {
267  outFile << "Base process for negatively charged particle capture at rest.\n";
268 }
269 
270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double condition(const G4ErrorSymMatrix &m)
void DeRegisterExtraProcess(G4VProcess *)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
virtual G4Element * SelectZandA(const G4Track &track, G4Nucleus *)
void AddSecondaries(const std::vector< G4HadSecondary > &addSecs)
virtual G4bool IsApplicable(const G4ParticleDefinition &)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4HadSecondary * GetSecondary(size_t i)
static G4HadronicProcessStore * Instance()
std::ofstream outFile
Definition: GammaRayTel.cc:68
G4String name
Definition: TRTMaterials.hh:40
const G4ThreeVector & GetPosition() const
virtual G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4HadProjectile thePro
G4HadronicInteraction * fEmCascade
void CheckEnergyMomentumConservation(const G4Track &, const G4Nucleus &)
const G4String & GetModelName() const
G4double GetTime() const
G4bool enableAtRestDoIt
Definition: G4VProcess.hh:350
int G4int
Definition: G4Types.hh:78
void SetWeight(G4double aValue)
virtual G4double AtRestGetPhysicalInteractionLength(const G4Track &track, G4ForceCondition *condition)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4HadronStoppingProcess(const G4String &name="hadronCaptureAtRest")
void ProposeWeight(G4double finalWeight)
void DumpState(const G4Track &, const G4String &, G4ExceptionDescription &)
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
void RegisterParticleForExtraProcess(G4VProcess *, const G4ParticleDefinition *)
G4ParticleChange * theTotalResult
bool G4bool
Definition: G4Types.hh:79
G4HadronicInteraction * ChooseHadronicInteraction(G4double kineticEnergy, G4Material *aMaterial, G4Element *anElement)
G4HadFinalState * CheckResult(const G4HadProjectile &thePro, const G4Nucleus &targetNucleus, G4HadFinalState *result) const
Definition: G4Step.hh:76
G4bool enablePostStepDoIt
Definition: G4VProcess.hh:352
const G4int n
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
G4double GetGlobalTime() const
const G4TouchableHandle & GetTouchableHandle() const
G4HadronicInteraction * fBoundDecay
G4Material * GetMaterial() const
void Initialise(const G4Track &aT)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ElementSelector * fElementSelector
void RegisterExtraProcess(G4VProcess *)
virtual void Initialize(const G4Track &)
G4Nucleus * GetTargetNucleusPointer()
void SetNumberOfSecondaries(G4int totSecondaries)
G4DynamicParticle * GetParticle()
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void SetBoundEnergy(G4double e)
virtual void ProcessDescription(std::ostream &outFile) const
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
#define G4endl
Definition: G4ios.hh:61
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
const G4String & GetName() const
Definition: G4Element.hh:127
G4ForceCondition
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
G4int GetNumberOfSecondaries() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
G4double GetLocalEnergyDeposit() const
G4HadFinalStateStatus GetStatusChange() const
void PrintInfo(const G4ParticleDefinition *)
G4double GetWeight() const