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