| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4HadronStoppingProcess.cc 88993 2015-03-17 11:17:13Z gcosmo $
26 // 27 //
27 //-------------------------------------------- 28 //---------------------------------------------------------------------
28 // 29 //
29 // GEANT4 Class 30 // GEANT4 Class
30 // 31 //
31 // File name: G4HadronStoppingProcess 32 // File name: G4HadronStoppingProcess
32 // 33 //
33 // Author V.Ivanchenko 21 April 2012 34 // Author V.Ivanchenko 21 April 2012
34 // 35 //
35 // 36 //
36 // Class Description: 37 // Class Description:
37 // 38 //
38 // Base process class for nuclear capture of n 39 // Base process class for nuclear capture of negatively charged particles
39 // 40 //
40 // Modifications: 41 // Modifications:
41 // 42 //
42 // 20120522 M. Kelsey -- Set enableAtRestDoI 43 // 20120522 M. Kelsey -- Set enableAtRestDoIt flag for G4ProcessManager
43 // 20120914 M. Kelsey -- Pass subType in bas 44 // 20120914 M. Kelsey -- Pass subType in base ctor, remove enable flags
44 // 20121004 K. Genser -- use G4HadronicProce 45 // 20121004 K. Genser -- use G4HadronicProcessType in the constructor
45 // 20121016 K. Genser -- Reverting to use on 46 // 20121016 K. Genser -- Reverting to use one argument c'tor
46 // 20140818 K. Genser -- Labeled tracks usin 47 // 20140818 K. Genser -- Labeled tracks using G4PhysicsModelCatalog
47 // 48 //
48 //-------------------------------------------- 49 //------------------------------------------------------------------------
49 50
50 #include "G4HadronStoppingProcess.hh" 51 #include "G4HadronStoppingProcess.hh"
51 #include "G4HadronicProcessStore.hh" 52 #include "G4HadronicProcessStore.hh"
52 #include "G4HadronicProcessType.hh" 53 #include "G4HadronicProcessType.hh"
53 #include "G4EmCaptureCascade.hh" 54 #include "G4EmCaptureCascade.hh"
54 #include "G4Nucleus.hh" 55 #include "G4Nucleus.hh"
55 #include "G4HadFinalState.hh" 56 #include "G4HadFinalState.hh"
56 #include "G4HadProjectile.hh" 57 #include "G4HadProjectile.hh"
57 #include "G4HadSecondary.hh" 58 #include "G4HadSecondary.hh"
58 #include "G4Material.hh" 59 #include "G4Material.hh"
59 #include "G4Element.hh" 60 #include "G4Element.hh"
60 61
61 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 63
63 G4HadronStoppingProcess::G4HadronStoppingProce 64 G4HadronStoppingProcess::G4HadronStoppingProcess(const G4String& name)
64 : G4HadronicProcess(name, fHadronAtRest), 65 : G4HadronicProcess(name, fHadronAtRest),
65 fElementSelector(new G4ElementSelector() 66 fElementSelector(new G4ElementSelector()),
66 fEmCascade(new G4EmCaptureCascade()), / 67 fEmCascade(new G4EmCaptureCascade()), // Owned by InteractionRegistry
67 fBoundDecay(0), 68 fBoundDecay(0),
68 emcID(-1), 69 emcID(-1),
69 ncID(-1), 70 ncID(-1),
70 dioID(-1) 71 dioID(-1)
71 { 72 {
72 // Modify G4VProcess flags to emulate G4VRes 73 // Modify G4VProcess flags to emulate G4VRest instead of G4VDiscrete
73 enableAtRestDoIt = true; 74 enableAtRestDoIt = true;
74 enablePostStepDoIt = false; 75 enablePostStepDoIt = false;
75 76
76 G4HadronicProcessStore::Instance()->Register 77 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
77 } 78 }
78 79
79 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 81
81 G4HadronStoppingProcess::~G4HadronStoppingProc 82 G4HadronStoppingProcess::~G4HadronStoppingProcess()
82 { 83 {
83 //G4HadronicProcessStore::Instance()->DeRegi << 84 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
84 delete fElementSelector; 85 delete fElementSelector;
85 // NOTE: fEmCascade and fEmBoundDecay owned 86 // NOTE: fEmCascade and fEmBoundDecay owned by registry, not locally
86 } 87 }
87 88
88 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 90
90 G4bool G4HadronStoppingProcess::IsApplicable(c 91 G4bool G4HadronStoppingProcess::IsApplicable(const G4ParticleDefinition& p)
91 { 92 {
92 return (p.GetPDGCharge() < 0.0); 93 return (p.GetPDGCharge() < 0.0);
93 } 94 }
94 95
95 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 97
97 void 98 void
98 G4HadronStoppingProcess::PreparePhysicsTable(c 99 G4HadronStoppingProcess::PreparePhysicsTable(const G4ParticleDefinition& p)
99 { 100 {
100 G4HadronicProcessStore::Instance()->Register 101 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this,&p);
101 emcID = G4PhysicsModelCatalog::GetModelID(G4 << 102 emcID = G4PhysicsModelCatalog::Register(G4String((GetProcessName() + "_EMCascade")));
102 ncID = G4PhysicsModelCatalog::GetModelID(G4 << 103 ncID = G4PhysicsModelCatalog::Register(G4String((GetProcessName() + "_NuclearCapture")));
103 dioID = G4PhysicsModelCatalog::GetModelID(G4 << 104 dioID = G4PhysicsModelCatalog::Register(G4String((GetProcessName() + "_DIO")));
104 } 105 }
105 106
106 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 108
108 void G4HadronStoppingProcess::BuildPhysicsTabl 109 void G4HadronStoppingProcess::BuildPhysicsTable(const G4ParticleDefinition& p)
109 { 110 {
110 G4HadronicProcessStore::Instance()->PrintInf 111 G4HadronicProcessStore::Instance()->PrintInfo(&p);
111 } 112 }
112 113
113 //....oooOO0OOooo........oooOO0OOooo........oo 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114 115
115 G4double G4HadronStoppingProcess::AtRestGetPhy 116 G4double G4HadronStoppingProcess::AtRestGetPhysicalInteractionLength(
116 const G4Track&, G4ForceCondition* conditio 117 const G4Track&, G4ForceCondition* condition)
117 { 118 {
118 *condition = NotForced; 119 *condition = NotForced;
119 return 0.0; 120 return 0.0;
120 } 121 }
121 122
122 //....oooOO0OOooo........oooOO0OOooo........oo 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 124
124 G4double G4HadronStoppingProcess::PostStepGetP 125 G4double G4HadronStoppingProcess::PostStepGetPhysicalInteractionLength(
125 const G4Track&, G4double, G4ForceCondition 126 const G4Track&, G4double, G4ForceCondition* condition)
126 { 127 {
127 *condition = NotForced; 128 *condition = NotForced;
128 return DBL_MAX; 129 return DBL_MAX;
129 } 130 }
130 131
131 //....oooOO0OOooo........oooOO0OOooo........oo 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132 133
133 G4VParticleChange* G4HadronStoppingProcess::At 134 G4VParticleChange* G4HadronStoppingProcess::AtRestDoIt(const G4Track& track,
134 const G4Step&) 135 const G4Step&)
135 { 136 {
136 // if primary is not Alive then do nothing 137 // if primary is not Alive then do nothing
137 theTotalResult->Initialize(track); 138 theTotalResult->Initialize(track);
138 139
139 G4Nucleus* nucleus = GetTargetNucleusPointer 140 G4Nucleus* nucleus = GetTargetNucleusPointer();
140 const G4Element* elm = fElementSelector->Sel << 141 G4Element* elm = fElementSelector->SelectZandA(track, nucleus);
141 142
142 G4HadFinalState* result = 0; 143 G4HadFinalState* result = 0;
143 thePro.Initialise(track); 144 thePro.Initialise(track);
144 145
145 // save track time an dstart capture from ze 146 // save track time an dstart capture from zero time
146 thePro.SetGlobalTime(0.0); 147 thePro.SetGlobalTime(0.0);
147 G4double time0 = track.GetGlobalTime(); 148 G4double time0 = track.GetGlobalTime();
148 149
149 G4bool nuclearCapture = true; 150 G4bool nuclearCapture = true;
150 151
151 // Do the electromagnetic cascade in the nuc 152 // Do the electromagnetic cascade in the nuclear field.
152 // EM cascade should keep G4HadFinalState ob 153 // EM cascade should keep G4HadFinalState object,
153 // because it will not be deleted at the end 154 // because it will not be deleted at the end of this method
154 // 155 //
155 result = fEmCascade->ApplyYourself(thePro, * 156 result = fEmCascade->ApplyYourself(thePro, *nucleus);
156 G4double ebound = result->GetLocalEnergyDepo 157 G4double ebound = result->GetLocalEnergyDeposit();
157 G4double edep = 0.0; 158 G4double edep = 0.0;
158 G4int nSecondaries = (G4int)result->GetNumbe << 159 G4int nSecondaries = result->GetNumberOfSecondaries();
159 G4int nEmCascadeSec = nSecondaries; 160 G4int nEmCascadeSec = nSecondaries;
160 161
161 // Try decay from bound level 162 // Try decay from bound level
162 // For mu- the time of projectile should be 163 // For mu- the time of projectile should be changed.
163 // Decay should keep G4HadFinalState object, 164 // Decay should keep G4HadFinalState object,
164 // because it will not be deleted at the end 165 // because it will not be deleted at the end of this method.
165 // 166 //
166 thePro.SetBoundEnergy(ebound); 167 thePro.SetBoundEnergy(ebound);
167 if(fBoundDecay) { 168 if(fBoundDecay) {
168 G4HadFinalState* resultDecay = 169 G4HadFinalState* resultDecay =
169 fBoundDecay->ApplyYourself(thePro, *nucl 170 fBoundDecay->ApplyYourself(thePro, *nucleus);
170 G4int n = (G4int)resultDecay->GetNumberOfS << 171 G4int n = resultDecay->GetNumberOfSecondaries();
171 if(0 < n) { 172 if(0 < n) {
172 nSecondaries += n; 173 nSecondaries += n;
173 result->AddSecondaries(resultDecay); 174 result->AddSecondaries(resultDecay);
174 } 175 }
175 if(resultDecay->GetStatusChange() == stopA 176 if(resultDecay->GetStatusChange() == stopAndKill) {
176 nuclearCapture = false; 177 nuclearCapture = false;
177 } 178 }
178 resultDecay->Clear(); 179 resultDecay->Clear();
179 } 180 }
180 181
181 if(nuclearCapture) { 182 if(nuclearCapture) {
182 183
183 // delay of capture 184 // delay of capture
184 G4double capTime = thePro.GetGlobalTime(); 185 G4double capTime = thePro.GetGlobalTime();
185 thePro.SetGlobalTime(0.0); 186 thePro.SetGlobalTime(0.0);
186 187
187 // select model 188 // select model
188 G4HadronicInteraction* model = 0; 189 G4HadronicInteraction* model = 0;
189 try { 190 try {
190 model = ChooseHadronicInteraction(thePro 191 model = ChooseHadronicInteraction(thePro, *nucleus,
191 track.GetMaterial(), elm); 192 track.GetMaterial(), elm);
192 } 193 }
193 catch(G4HadronicException & aE) { 194 catch(G4HadronicException & aE) {
194 G4ExceptionDescription ed; 195 G4ExceptionDescription ed;
195 ed << "Target element "<<elm->GetName()< 196 ed << "Target element "<<elm->GetName()<<" Z= "
196 << nucleus->GetZ_asInt() << " A= " 197 << nucleus->GetZ_asInt() << " A= "
197 << nucleus->GetA_asInt() << G4endl; 198 << nucleus->GetA_asInt() << G4endl;
198 DumpState(track,"ChooseHadronicInteracti 199 DumpState(track,"ChooseHadronicInteraction",ed);
199 ed << " No HadronicInteraction found out 200 ed << " No HadronicInteraction found out" << G4endl;
200 G4Exception("G4HadronStoppingProcess::At 201 G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had005",
201 FatalException, ed); 202 FatalException, ed);
202 } 203 }
203 204
204 G4HadFinalState* resultNuc = 0; 205 G4HadFinalState* resultNuc = 0;
205 G4int reentryCount = 0; 206 G4int reentryCount = 0;
206 do { 207 do {
207 // sample final state 208 // sample final state
208 // nuclear interaction should keep G4Had 209 // nuclear interaction should keep G4HadFinalState object
209 // model should define time of each seco 210 // model should define time of each secondary particle
210 try { 211 try {
211 resultNuc = model->ApplyYourself(thePro, *nu 212 resultNuc = model->ApplyYourself(thePro, *nucleus);
212 ++reentryCount; 213 ++reentryCount;
213 } 214 }
214 catch(G4HadronicException & aR) { << 215 catch(G4HadronicException aR) {
215 G4ExceptionDescription ed; 216 G4ExceptionDescription ed;
216 ed << "Call for " << model->GetModelName() < 217 ed << "Call for " << model->GetModelName() << G4endl;
217 ed << "Target element "<<elm->GetName()<<" 218 ed << "Target element "<<elm->GetName()<<" Z= "
218 << nucleus->GetZ_asInt() 219 << nucleus->GetZ_asInt()
219 << " A= " << nucleus->GetA_asInt() << G4 220 << " A= " << nucleus->GetA_asInt() << G4endl;
220 DumpState(track,"ApplyYourself",ed); 221 DumpState(track,"ApplyYourself",ed);
221 ed << " ApplyYourself failed" << G4endl; 222 ed << " ApplyYourself failed" << G4endl;
222 G4Exception("G4HadronStoppingProcess::AtRest 223 G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
223 FatalException, ed); 224 FatalException, ed);
224 } 225 }
225 226
226 // Check the result for catastrophic ene 227 // Check the result for catastrophic energy non-conservation
227 resultNuc = CheckResult(thePro, *nucleus 228 resultNuc = CheckResult(thePro, *nucleus, resultNuc);
228 229
229 if(reentryCount>100) { 230 if(reentryCount>100) {
230 G4ExceptionDescription ed; 231 G4ExceptionDescription ed;
231 ed << "Call for " << model->GetModelName() < 232 ed << "Call for " << model->GetModelName() << G4endl;
232 ed << "Target element "<<elm->GetName()<<" 233 ed << "Target element "<<elm->GetName()<<" Z= "
233 << nucleus->GetZ_asInt() 234 << nucleus->GetZ_asInt()
234 << " A= " << nucleus->GetA_asInt() << G4 235 << " A= " << nucleus->GetA_asInt() << G4endl;
235 DumpState(track,"ApplyYourself",ed); 236 DumpState(track,"ApplyYourself",ed);
236 ed << " ApplyYourself does not completed aft 237 ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
237 G4Exception("G4HadronStoppingProcess::AtRest 238 G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
238 FatalException, ed); 239 FatalException, ed);
239 } 240 }
240 // Loop checking, 06-Aug-2015, Vladimir << 241 }
241 } while(!resultNuc); << 242 while(!resultNuc);
242 243
243 edep = resultNuc->GetLocalEnergyDeposit(); 244 edep = resultNuc->GetLocalEnergyDeposit();
244 std::size_t nnuc = resultNuc->GetNumberOfS << 245 size_t nnuc = resultNuc->GetNumberOfSecondaries();
245 246
246 // add delay time of capture 247 // add delay time of capture
247 for(std::size_t i=0; i<nnuc; ++i) { << 248 for(size_t i=0; i<nnuc; ++i) {
248 G4HadSecondary* sec = resultNuc->GetSeco 249 G4HadSecondary* sec = resultNuc->GetSecondary(i);
249 sec->SetTime(capTime + sec->GetTime()); 250 sec->SetTime(capTime + sec->GetTime());
250 } 251 }
251 252
252 nSecondaries += nnuc; 253 nSecondaries += nnuc;
253 result->AddSecondaries(resultNuc); 254 result->AddSecondaries(resultNuc);
254 resultNuc->Clear(); 255 resultNuc->Clear();
255 } 256 }
256 257
257 // Fill results 258 // Fill results
258 // 259 //
259 theTotalResult->ProposeTrackStatus(fStopAndK 260 theTotalResult->ProposeTrackStatus(fStopAndKill);
260 theTotalResult->ProposeLocalEnergyDeposit(ed 261 theTotalResult->ProposeLocalEnergyDeposit(edep);
261 theTotalResult->SetNumberOfSecondaries(nSeco 262 theTotalResult->SetNumberOfSecondaries(nSecondaries);
262 G4double w = track.GetWeight(); 263 G4double w = track.GetWeight();
263 theTotalResult->ProposeWeight(w); 264 theTotalResult->ProposeWeight(w);
264 for(G4int i=0; i<nSecondaries; ++i) { 265 for(G4int i=0; i<nSecondaries; ++i) {
265 G4HadSecondary* sec = result->GetSecondary 266 G4HadSecondary* sec = result->GetSecondary(i);
266 267
267 // add track global time to the reaction t 268 // add track global time to the reaction time
268 G4double time = sec->GetTime(); 269 G4double time = sec->GetTime();
269 if(time < 0.0) { time = 0.0; } 270 if(time < 0.0) { time = 0.0; }
270 time += time0; 271 time += time0;
271 272
272 // create secondary track 273 // create secondary track
273 G4Track* t = new G4Track(sec->GetParticle( 274 G4Track* t = new G4Track(sec->GetParticle(),
274 time, 275 time,
275 track.GetPosition()); 276 track.GetPosition());
276 t->SetWeight(w*sec->GetWeight()); 277 t->SetWeight(w*sec->GetWeight());
277 278
278 // use SetCreatorModelID to "label" the tr << 279 // use SetCreatorModelIndex to "label" the track
279 if (i<nEmCascadeSec) { 280 if (i<nEmCascadeSec) {
280 t->SetCreatorModelID(emcID); << 281 t->SetCreatorModelIndex(emcID);
281 } else if (nuclearCapture) { 282 } else if (nuclearCapture) {
282 t->SetCreatorModelID(ncID); << 283 t->SetCreatorModelIndex(ncID);
283 } else { 284 } else {
284 t->SetCreatorModelID(dioID); << 285 t->SetCreatorModelIndex(dioID);
285 } 286 }
286 287
287 t->SetTouchableHandle(track.GetTouchableHa 288 t->SetTouchableHandle(track.GetTouchableHandle());
288 theTotalResult->AddSecondary(t); 289 theTotalResult->AddSecondary(t);
289 } 290 }
290 result->Clear(); 291 result->Clear();
291 292
292 if (epReportLevel != 0) { 293 if (epReportLevel != 0) {
293 CheckEnergyMomentumConservation(track, *nu 294 CheckEnergyMomentumConservation(track, *nucleus);
294 } 295 }
295 return theTotalResult; 296 return theTotalResult;
296 } 297 }
297 298
298 //....oooOO0OOooo........oooOO0OOooo........oo 299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299 300
300 void G4HadronStoppingProcess::ProcessDescripti 301 void G4HadronStoppingProcess::ProcessDescription(std::ostream& outFile) const
301 { 302 {
302 outFile << "Base process for negatively char 303 outFile << "Base process for negatively charged particle capture at rest.\n";
303 } 304 }
304 305
305 //....oooOO0OOooo........oooOO0OOooo........oo 306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
306 307