| Geant4 Cross Reference |
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
50 #include "G4HadronStoppingProcess.hh"
51 #include "G4HadronicProcessStore.hh"
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
63 G4HadronStoppingProcess::G4HadronStoppingProcess(const G4String& name)
64 : G4HadronicProcess(name, fHadronAtRest),
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
76 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80
81 G4HadronStoppingProcess::~G4HadronStoppingProcess()
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
90 G4bool G4HadronStoppingProcess::IsApplicable(const G4ParticleDefinition& p)
91 {
92 return (p.GetPDGCharge() < 0.0);
93 }
94
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96
97 void
98 G4HadronStoppingProcess::PreparePhysicsTable(const G4ParticleDefinition& p)
99 {
100 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this,&p);
101 emcID = G4PhysicsModelCatalog::GetModelID(G4String("model_" + (GetProcessName() + "_EMCascade")));
102 ncID = G4PhysicsModelCatalog::GetModelID(G4String("model_" + (GetProcessName() + "_NuclearCapture")));
103 dioID = G4PhysicsModelCatalog::GetModelID(G4String("model_" + (GetProcessName() + "_DIO")));
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107
108 void G4HadronStoppingProcess::BuildPhysicsTable(const G4ParticleDefinition& p)
109 {
110 G4HadronicProcessStore::Instance()->PrintInfo(&p);
111 }
112
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114
115 G4double G4HadronStoppingProcess::AtRestGetPhysicalInteractionLength(
116 const G4Track&, G4ForceCondition* condition)
117 {
118 *condition = NotForced;
119 return 0.0;
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123
124 G4double G4HadronStoppingProcess::PostStepGetPhysicalInteractionLength(
125 const G4Track&, G4double, G4ForceCondition* condition)
126 {
127 *condition = NotForced;
128 return DBL_MAX;
129 }
130
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132
133 G4VParticleChange* G4HadronStoppingProcess::AtRestDoIt(const G4Track& track,
134 const G4Step&)
135 {
136 // if primary is not Alive then do nothing
137 theTotalResult->Initialize(track);
138
139 G4Nucleus* nucleus = GetTargetNucleusPointer();
140 const 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 = (G4int)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 = (G4int)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
188 G4HadronicInteraction* model = 0;
189 try {
190 model = ChooseHadronicInteraction(thePro, *nucleus,
191 track.GetMaterial(), elm);
192 }
193 catch(G4HadronicException & aE) {
194 G4ExceptionDescription ed;
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) {
215 G4ExceptionDescription ed;
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) {
230 G4ExceptionDescription ed;
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 std::size_t nnuc = resultNuc->GetNumberOfSecondaries();
245
246 // add delay time of capture
247 for(std::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 //
259 theTotalResult->ProposeTrackStatus(fStopAndKill);
260 theTotalResult->ProposeLocalEnergyDeposit(edep);
261 theTotalResult->SetNumberOfSecondaries(nSecondaries);
262 G4double w = track.GetWeight();
263 theTotalResult->ProposeWeight(w);
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 SetCreatorModelID to "label" the track
279 if (i<nEmCascadeSec) {
280 t->SetCreatorModelID(emcID);
281 } else if (nuclearCapture) {
282 t->SetCreatorModelID(ncID);
283 } else {
284 t->SetCreatorModelID(dioID);
285 }
286
287 t->SetTouchableHandle(track.GetTouchableHandle());
288 theTotalResult->AddSecondary(t);
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....
306