| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 n
39 //
40 // Modifications:
41 //
42 // 20120522 M. Kelsey -- Set enableAtRestDoI
43 // 20120914 M. Kelsey -- Pass subType in bas
44 // 20121004 K. Genser -- use G4HadronicProce
45 // 20121016 K. Genser -- Reverting to use on
46 // 20140818 K. Genser -- Labeled tracks usin
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........oo
62
63 G4HadronStoppingProcess::G4HadronStoppingProce
64 : G4HadronicProcess(name, fHadronAtRest),
65 fElementSelector(new G4ElementSelector()
66 fEmCascade(new G4EmCaptureCascade()), /
67 fBoundDecay(0),
68 emcID(-1),
69 ncID(-1),
70 dioID(-1)
71 {
72 // Modify G4VProcess flags to emulate G4VRes
73 enableAtRestDoIt = true;
74 enablePostStepDoIt = false;
75
76 G4HadronicProcessStore::Instance()->Register
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oo
80
81 G4HadronStoppingProcess::~G4HadronStoppingProc
82 {
83 //G4HadronicProcessStore::Instance()->DeRegi
84 delete fElementSelector;
85 // NOTE: fEmCascade and fEmBoundDecay owned
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oo
89
90 G4bool G4HadronStoppingProcess::IsApplicable(c
91 {
92 return (p.GetPDGCharge() < 0.0);
93 }
94
95 //....oooOO0OOooo........oooOO0OOooo........oo
96
97 void
98 G4HadronStoppingProcess::PreparePhysicsTable(c
99 {
100 G4HadronicProcessStore::Instance()->Register
101 emcID = G4PhysicsModelCatalog::GetModelID(G4
102 ncID = G4PhysicsModelCatalog::GetModelID(G4
103 dioID = G4PhysicsModelCatalog::GetModelID(G4
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 void G4HadronStoppingProcess::BuildPhysicsTabl
109 {
110 G4HadronicProcessStore::Instance()->PrintInf
111 }
112
113 //....oooOO0OOooo........oooOO0OOooo........oo
114
115 G4double G4HadronStoppingProcess::AtRestGetPhy
116 const G4Track&, G4ForceCondition* conditio
117 {
118 *condition = NotForced;
119 return 0.0;
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oo
123
124 G4double G4HadronStoppingProcess::PostStepGetP
125 const G4Track&, G4double, G4ForceCondition
126 {
127 *condition = NotForced;
128 return DBL_MAX;
129 }
130
131 //....oooOO0OOooo........oooOO0OOooo........oo
132
133 G4VParticleChange* G4HadronStoppingProcess::At
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->Sel
141
142 G4HadFinalState* result = 0;
143 thePro.Initialise(track);
144
145 // save track time an dstart capture from ze
146 thePro.SetGlobalTime(0.0);
147 G4double time0 = track.GetGlobalTime();
148
149 G4bool nuclearCapture = true;
150
151 // Do the electromagnetic cascade in the nuc
152 // EM cascade should keep G4HadFinalState ob
153 // because it will not be deleted at the end
154 //
155 result = fEmCascade->ApplyYourself(thePro, *
156 G4double ebound = result->GetLocalEnergyDepo
157 G4double edep = 0.0;
158 G4int nSecondaries = (G4int)result->GetNumbe
159 G4int nEmCascadeSec = nSecondaries;
160
161 // Try decay from bound level
162 // For mu- the time of projectile should be
163 // Decay should keep G4HadFinalState object,
164 // because it will not be deleted at the end
165 //
166 thePro.SetBoundEnergy(ebound);
167 if(fBoundDecay) {
168 G4HadFinalState* resultDecay =
169 fBoundDecay->ApplyYourself(thePro, *nucl
170 G4int n = (G4int)resultDecay->GetNumberOfS
171 if(0 < n) {
172 nSecondaries += n;
173 result->AddSecondaries(resultDecay);
174 }
175 if(resultDecay->GetStatusChange() == stopA
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
191 track.GetMaterial(), elm);
192 }
193 catch(G4HadronicException & aE) {
194 G4ExceptionDescription ed;
195 ed << "Target element "<<elm->GetName()<
196 << nucleus->GetZ_asInt() << " A= "
197 << nucleus->GetA_asInt() << G4endl;
198 DumpState(track,"ChooseHadronicInteracti
199 ed << " No HadronicInteraction found out
200 G4Exception("G4HadronStoppingProcess::At
201 FatalException, ed);
202 }
203
204 G4HadFinalState* resultNuc = 0;
205 G4int reentryCount = 0;
206 do {
207 // sample final state
208 // nuclear interaction should keep G4Had
209 // model should define time of each seco
210 try {
211 resultNuc = model->ApplyYourself(thePro, *nu
212 ++reentryCount;
213 }
214 catch(G4HadronicException & aR) {
215 G4ExceptionDescription ed;
216 ed << "Call for " << model->GetModelName() <
217 ed << "Target element "<<elm->GetName()<<"
218 << nucleus->GetZ_asInt()
219 << " A= " << nucleus->GetA_asInt() << G4
220 DumpState(track,"ApplyYourself",ed);
221 ed << " ApplyYourself failed" << G4endl;
222 G4Exception("G4HadronStoppingProcess::AtRest
223 FatalException, ed);
224 }
225
226 // Check the result for catastrophic ene
227 resultNuc = CheckResult(thePro, *nucleus
228
229 if(reentryCount>100) {
230 G4ExceptionDescription ed;
231 ed << "Call for " << model->GetModelName() <
232 ed << "Target element "<<elm->GetName()<<"
233 << nucleus->GetZ_asInt()
234 << " A= " << nucleus->GetA_asInt() << G4
235 DumpState(track,"ApplyYourself",ed);
236 ed << " ApplyYourself does not completed aft
237 G4Exception("G4HadronStoppingProcess::AtRest
238 FatalException, ed);
239 }
240 // Loop checking, 06-Aug-2015, Vladimir
241 } while(!resultNuc);
242
243 edep = resultNuc->GetLocalEnergyDeposit();
244 std::size_t nnuc = resultNuc->GetNumberOfS
245
246 // add delay time of capture
247 for(std::size_t i=0; i<nnuc; ++i) {
248 G4HadSecondary* sec = resultNuc->GetSeco
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(fStopAndK
260 theTotalResult->ProposeLocalEnergyDeposit(ed
261 theTotalResult->SetNumberOfSecondaries(nSeco
262 G4double w = track.GetWeight();
263 theTotalResult->ProposeWeight(w);
264 for(G4int i=0; i<nSecondaries; ++i) {
265 G4HadSecondary* sec = result->GetSecondary
266
267 // add track global time to the reaction t
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 tr
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.GetTouchableHa
288 theTotalResult->AddSecondary(t);
289 }
290 result->Clear();
291
292 if (epReportLevel != 0) {
293 CheckEnergyMomentumConservation(track, *nu
294 }
295 return theTotalResult;
296 }
297
298 //....oooOO0OOooo........oooOO0OOooo........oo
299
300 void G4HadronStoppingProcess::ProcessDescripti
301 {
302 outFile << "Base process for negatively char
303 }
304
305 //....oooOO0OOooo........oooOO0OOooo........oo
306