Geant4  10.02.p01
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 94351 2015-11-12 15:35:32Z 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 // 20140818 K. Genser -- Labeled tracks using G4PhysicsModelCatalog
48 //
49 //------------------------------------------------------------------------
50 
53 #include "G4HadronicProcessType.hh"
54 #include "G4EmCaptureCascade.hh"
55 #include "G4Nucleus.hh"
56 #include "G4HadFinalState.hh"
57 #include "G4HadProjectile.hh"
58 #include "G4HadSecondary.hh"
59 #include "G4Material.hh"
60 #include "G4Element.hh"
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 
66  fElementSelector(new G4ElementSelector()),
67  fEmCascade(new G4EmCaptureCascade()), // Owned by InteractionRegistry
68  fBoundDecay(0),
69  emcID(-1),
70  ncID(-1),
71  dioID(-1)
72 {
73  // Modify G4VProcess flags to emulate G4VRest instead of G4VDiscrete
74  enableAtRestDoIt = true;
75  enablePostStepDoIt = false;
76 
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
83 {
84  //G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
85  delete fElementSelector;
86  // NOTE: fEmCascade and fEmBoundDecay owned by registry, not locally
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
92 {
93  return (p.GetPDGCharge() < 0.0);
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
98 void
100 {
103  ncID = G4PhysicsModelCatalog::Register(G4String((GetProcessName() + "_NuclearCapture")));
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108 
110 {
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115 
118 {
119  *condition = NotForced;
120  return 0.0;
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 
127 {
128  *condition = NotForced;
129  return DBL_MAX;
130 }
131 
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133 
135  const G4Step&)
136 {
137  // if primary is not Alive then do nothing
138  theTotalResult->Initialize(track);
139 
140  G4Nucleus* nucleus = GetTargetNucleusPointer();
141  G4Element* elm = fElementSelector->SelectZandA(track, nucleus);
142 
143  G4HadFinalState* result = 0;
144  thePro.Initialise(track);
145 
146  // save track time an dstart capture from zero time
147  thePro.SetGlobalTime(0.0);
148  G4double time0 = track.GetGlobalTime();
149 
150  G4bool nuclearCapture = true;
151 
152  // Do the electromagnetic cascade in the nuclear field.
153  // EM cascade should keep G4HadFinalState object,
154  // because it will not be deleted at the end of this method
155  //
156  result = fEmCascade->ApplyYourself(thePro, *nucleus);
157  G4double ebound = result->GetLocalEnergyDeposit();
158  G4double edep = 0.0;
159  G4int nSecondaries = result->GetNumberOfSecondaries();
160  G4int nEmCascadeSec = nSecondaries;
161 
162  // Try decay from bound level
163  // For mu- the time of projectile should be changed.
164  // Decay should keep G4HadFinalState object,
165  // because it will not be deleted at the end of this method.
166  //
167  thePro.SetBoundEnergy(ebound);
168  if(fBoundDecay) {
169  G4HadFinalState* resultDecay =
170  fBoundDecay->ApplyYourself(thePro, *nucleus);
171  G4int n = resultDecay->GetNumberOfSecondaries();
172  if(0 < n) {
173  nSecondaries += n;
174  result->AddSecondaries(resultDecay);
175  }
176  if(resultDecay->GetStatusChange() == stopAndKill) {
177  nuclearCapture = false;
178  }
179  resultDecay->Clear();
180  }
181 
182  if(nuclearCapture) {
183 
184  // delay of capture
185  G4double capTime = thePro.GetGlobalTime();
186  thePro.SetGlobalTime(0.0);
187 
188  // select model
189  G4HadronicInteraction* model = 0;
190  try {
191  model = ChooseHadronicInteraction(thePro, *nucleus,
192  track.GetMaterial(), elm);
193  }
194  catch(G4HadronicException & aE) {
196  ed << "Target element "<<elm->GetName()<<" Z= "
197  << nucleus->GetZ_asInt() << " A= "
198  << nucleus->GetA_asInt() << G4endl;
199  DumpState(track,"ChooseHadronicInteraction",ed);
200  ed << " No HadronicInteraction found out" << G4endl;
201  G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had005",
202  FatalException, ed);
203  }
204 
205  G4HadFinalState* resultNuc = 0;
206  G4int reentryCount = 0;
207  do {
208  // sample final state
209  // nuclear interaction should keep G4HadFinalState object
210  // model should define time of each secondary particle
211  try {
212  resultNuc = model->ApplyYourself(thePro, *nucleus);
213  ++reentryCount;
214  }
215  catch(G4HadronicException aR) {
217  ed << "Call for " << model->GetModelName() << G4endl;
218  ed << "Target element "<<elm->GetName()<<" Z= "
219  << nucleus->GetZ_asInt()
220  << " A= " << nucleus->GetA_asInt() << G4endl;
221  DumpState(track,"ApplyYourself",ed);
222  ed << " ApplyYourself failed" << G4endl;
223  G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
224  FatalException, ed);
225  }
226 
227  // Check the result for catastrophic energy non-conservation
228  resultNuc = CheckResult(thePro, *nucleus, resultNuc);
229 
230  if(reentryCount>100) {
232  ed << "Call for " << model->GetModelName() << G4endl;
233  ed << "Target element "<<elm->GetName()<<" Z= "
234  << nucleus->GetZ_asInt()
235  << " A= " << nucleus->GetA_asInt() << G4endl;
236  DumpState(track,"ApplyYourself",ed);
237  ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
238  G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
239  FatalException, ed);
240  }
241  // Loop checking, 06-Aug-2015, Vladimir Ivanchenko
242  } while(!resultNuc);
243 
244  edep = resultNuc->GetLocalEnergyDeposit();
245  size_t nnuc = resultNuc->GetNumberOfSecondaries();
246 
247  // add delay time of capture
248  for(size_t i=0; i<nnuc; ++i) {
249  G4HadSecondary* sec = resultNuc->GetSecondary(i);
250  sec->SetTime(capTime + sec->GetTime());
251  }
252 
253  nSecondaries += nnuc;
254  result->AddSecondaries(resultNuc);
255  resultNuc->Clear();
256  }
257 
258  // Fill results
259  //
262  theTotalResult->SetNumberOfSecondaries(nSecondaries);
263  G4double w = track.GetWeight();
265  for(G4int i=0; i<nSecondaries; ++i) {
266  G4HadSecondary* sec = result->GetSecondary(i);
267 
268  // add track global time to the reaction time
269  G4double time = sec->GetTime();
270  if(time < 0.0) { time = 0.0; }
271  time += time0;
272 
273  // create secondary track
274  G4Track* t = new G4Track(sec->GetParticle(),
275  time,
276  track.GetPosition());
277  t->SetWeight(w*sec->GetWeight());
278 
279  // use SetCreatorModelIndex to "label" the track
280  if (i<nEmCascadeSec) {
282  } else if (nuclearCapture) {
284  } else {
286  }
287 
290  }
291  result->Clear();
292 
293  if (epReportLevel != 0) {
294  CheckEnergyMomentumConservation(track, *nucleus);
295  }
296  return theTotalResult;
297 }
298 
299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
300 
301 void G4HadronStoppingProcess::ProcessDescription(std::ostream& outFile) const
302 {
303  outFile << "Base process for negatively charged particle capture at rest.\n";
304 }
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double condition(const G4ErrorSymMatrix &m)
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()
G4String name
Definition: TRTMaterials.hh:40
const G4ThreeVector & GetPosition() const
const G4double w[NPOINTSGL]
virtual G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
void SetTouchableHandle(const G4TouchableHandle &apValue)
void SetGlobalTime(G4double t)
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)
void SetCreatorModelIndex(G4int idx)
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
G4double GetGlobalTime() const
Definition: G4Step.hh:76
G4bool enablePostStepDoIt
Definition: G4VProcess.hh:352
const G4int n
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
G4double GetGlobalTime() const
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
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 &)
void SetTime(G4double aT)
G4Nucleus * GetTargetNucleusPointer()
void SetNumberOfSecondaries(G4int totSecondaries)
static G4int Register(const G4String &)
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
G4HadronicInteraction * ChooseHadronicInteraction(const G4HadProjectile &aHadProjectile, G4Nucleus &aTargetNucleus, G4Material *aMaterial, G4Element *anElement)
#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
G4HadFinalState * CheckResult(const G4HadProjectile &thePro, const G4Nucleus &targetNucleus, G4HadFinalState *result)
G4HadFinalStateStatus GetStatusChange() const
void PrintInfo(const G4ParticleDefinition *)
G4double GetWeight() const