ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HadronStoppingProcess.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4HadronStoppingProcess.cc
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 //
27 //---------------------------------------------------------------------
28 //
29 // GEANT4 Class
30 //
31 // File name: G4HadronStoppingProcess
32 //
33 // Author V.Ivanchenko 21 April 2012
34 //
35 //
36 // Class Description:
37 //
38 // Base process class for nuclear capture of negatively charged particles
39 //
40 // Modifications:
41 //
42 // 20120522 M. Kelsey -- Set enableAtRestDoIt flag for G4ProcessManager
43 // 20120914 M. Kelsey -- Pass subType in base ctor, remove enable flags
44 // 20121004 K. Genser -- use G4HadronicProcessType in the constructor
45 // 20121016 K. Genser -- Reverting to use one argument c'tor
46 // 20140818 K. Genser -- Labeled tracks using G4PhysicsModelCatalog
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  fElementSelector(new G4ElementSelector()),
66  fEmCascade(new G4EmCaptureCascade()), // Owned by InteractionRegistry
67  fBoundDecay(0),
68  emcID(-1),
69  ncID(-1),
70  dioID(-1)
71 {
72  // Modify G4VProcess flags to emulate G4VRest instead of G4VDiscrete
73  enableAtRestDoIt = true;
74  enablePostStepDoIt = false;
75 
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
82 {
83  //G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
84  delete fElementSelector;
85  // NOTE: fEmCascade and fEmBoundDecay owned by registry, not locally
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
91 {
92  return (p.GetPDGCharge() < 0.0);
93 }
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 
97 void
99 {
102  ncID = G4PhysicsModelCatalog::Register(G4String((GetProcessName() + "_NuclearCapture")));
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 
109 {
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114 
117 {
118  *condition = NotForced;
119  return 0.0;
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 
126 {
127  *condition = NotForced;
128  return DBL_MAX;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132 
134  const G4Step&)
135 {
136  // if primary is not Alive then do nothing
137  theTotalResult->Initialize(track);
138 
139  G4Nucleus* nucleus = GetTargetNucleusPointer();
140  G4Element* elm = fElementSelector->SelectZandA(track, nucleus);
141 
142  G4HadFinalState* result = 0;
143  thePro.Initialise(track);
144 
145  // save track time an dstart capture from zero time
146  thePro.SetGlobalTime(0.0);
147  G4double time0 = track.GetGlobalTime();
148 
149  G4bool nuclearCapture = true;
150 
151  // Do the electromagnetic cascade in the nuclear field.
152  // EM cascade should keep G4HadFinalState object,
153  // because it will not be deleted at the end of this method
154  //
155  result = fEmCascade->ApplyYourself(thePro, *nucleus);
156  G4double ebound = result->GetLocalEnergyDeposit();
157  G4double edep = 0.0;
158  G4int nSecondaries = result->GetNumberOfSecondaries();
159  G4int nEmCascadeSec = nSecondaries;
160 
161  // Try decay from bound level
162  // For mu- the time of projectile should be changed.
163  // Decay should keep G4HadFinalState object,
164  // because it will not be deleted at the end of this method.
165  //
166  thePro.SetBoundEnergy(ebound);
167  if(fBoundDecay) {
168  G4HadFinalState* resultDecay =
169  fBoundDecay->ApplyYourself(thePro, *nucleus);
170  G4int n = resultDecay->GetNumberOfSecondaries();
171  if(0 < n) {
172  nSecondaries += n;
173  result->AddSecondaries(resultDecay);
174  }
175  if(resultDecay->GetStatusChange() == stopAndKill) {
176  nuclearCapture = false;
177  }
178  resultDecay->Clear();
179  }
180 
181  if(nuclearCapture) {
182 
183  // delay of capture
184  G4double capTime = thePro.GetGlobalTime();
185  thePro.SetGlobalTime(0.0);
186 
187  // select model
189  try {
190  model = ChooseHadronicInteraction(thePro, *nucleus,
191  track.GetMaterial(), elm);
192  }
193  catch(G4HadronicException & aE) {
195  ed << "Target element "<<elm->GetName()<<" Z= "
196  << nucleus->GetZ_asInt() << " A= "
197  << nucleus->GetA_asInt() << G4endl;
198  DumpState(track,"ChooseHadronicInteraction",ed);
199  ed << " No HadronicInteraction found out" << G4endl;
200  G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had005",
201  FatalException, ed);
202  }
203 
204  G4HadFinalState* resultNuc = 0;
205  G4int reentryCount = 0;
206  do {
207  // sample final state
208  // nuclear interaction should keep G4HadFinalState object
209  // model should define time of each secondary particle
210  try {
211  resultNuc = model->ApplyYourself(thePro, *nucleus);
212  ++reentryCount;
213  }
214  catch(G4HadronicException & aR) {
216  ed << "Call for " << model->GetModelName() << G4endl;
217  ed << "Target element "<<elm->GetName()<<" Z= "
218  << nucleus->GetZ_asInt()
219  << " A= " << nucleus->GetA_asInt() << G4endl;
220  DumpState(track,"ApplyYourself",ed);
221  ed << " ApplyYourself failed" << G4endl;
222  G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
223  FatalException, ed);
224  }
225 
226  // Check the result for catastrophic energy non-conservation
227  resultNuc = CheckResult(thePro, *nucleus, resultNuc);
228 
229  if(reentryCount>100) {
231  ed << "Call for " << model->GetModelName() << G4endl;
232  ed << "Target element "<<elm->GetName()<<" Z= "
233  << nucleus->GetZ_asInt()
234  << " A= " << nucleus->GetA_asInt() << G4endl;
235  DumpState(track,"ApplyYourself",ed);
236  ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
237  G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
238  FatalException, ed);
239  }
240  // Loop checking, 06-Aug-2015, Vladimir Ivanchenko
241  } while(!resultNuc);
242 
243  edep = resultNuc->GetLocalEnergyDeposit();
244  size_t nnuc = resultNuc->GetNumberOfSecondaries();
245 
246  // add delay time of capture
247  for(size_t i=0; i<nnuc; ++i) {
248  G4HadSecondary* sec = resultNuc->GetSecondary(i);
249  sec->SetTime(capTime + sec->GetTime());
250  }
251 
252  nSecondaries += nnuc;
253  result->AddSecondaries(resultNuc);
254  resultNuc->Clear();
255  }
256 
257  // Fill results
258  //
261  theTotalResult->SetNumberOfSecondaries(nSecondaries);
262  G4double w = track.GetWeight();
264  for(G4int i=0; i<nSecondaries; ++i) {
265  G4HadSecondary* sec = result->GetSecondary(i);
266 
267  // add track global time to the reaction time
268  G4double time = sec->GetTime();
269  if(time < 0.0) { time = 0.0; }
270  time += time0;
271 
272  // create secondary track
273  G4Track* t = new G4Track(sec->GetParticle(),
274  time,
275  track.GetPosition());
276  t->SetWeight(w*sec->GetWeight());
277 
278  // use SetCreatorModelIndex to "label" the track
279  if (i<nEmCascadeSec) {
281  } else if (nuclearCapture) {
283  } else {
285  }
286 
289  }
290  result->Clear();
291 
292  if (epReportLevel != 0) {
293  CheckEnergyMomentumConservation(track, *nucleus);
294  }
295  return theTotalResult;
296 }
297 
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299 
300 void G4HadronStoppingProcess::ProcessDescription(std::ostream& outFile) const
301 {
302  outFile << "Base process for negatively charged particle capture at rest.\n";
303 }
304 
305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....