| 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 // G4VParticleChange class implementation 26 // G4VParticleChange class implementation
27 // 27 //
28 // Author: Hisaya Kurashige, 23 March 1998 28 // Author: Hisaya Kurashige, 23 March 1998
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4VParticleChange.hh" 31 #include "G4VParticleChange.hh"
32 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
33 #include "G4ExceptionSeverity.hh" 33 #include "G4ExceptionSeverity.hh"
34 34
35 const G4double G4VParticleChange::accuracyForW 35 const G4double G4VParticleChange::accuracyForWarning = 1.0e-9;
36 const G4double G4VParticleChange::accuracyForE 36 const G4double G4VParticleChange::accuracyForException = 0.001;
37 const G4int G4VParticleChange::maxError = 10; 37 const G4int G4VParticleChange::maxError = 10;
38 38
39 // ------------------------------------------- 39 // --------------------------------------------------------------------
40 G4VParticleChange::G4VParticleChange() 40 G4VParticleChange::G4VParticleChange()
41 { 41 {
42 #ifdef G4VERBOSE 42 #ifdef G4VERBOSE
43 // activate CheckIt if in VERBOSE mode 43 // activate CheckIt if in VERBOSE mode
44 debugFlag = true; 44 debugFlag = true;
45 #endif 45 #endif
46 } 46 }
47 47
48 // ------------------------------------------- 48 // --------------------------------------------------------------------
49 void G4VParticleChange::AddSecondary(G4Track* 49 void G4VParticleChange::AddSecondary(G4Track* aTrack)
50 { 50 {
51 if(debugFlag) 51 if(debugFlag)
52 CheckSecondary(*aTrack); 52 CheckSecondary(*aTrack);
53 53
54 if(!fSetSecondaryWeightByProcess) 54 if(!fSetSecondaryWeightByProcess)
55 aTrack->SetWeight(theParentWeight); 55 aTrack->SetWeight(theParentWeight);
56 56
57 // add a secondary after size check 57 // add a secondary after size check
58 if(theSizeOftheListOfSecondaries > theNumber 58 if(theSizeOftheListOfSecondaries > theNumberOfSecondaries)
59 { 59 {
60 theListOfSecondaries[theNumberOfSecondarie 60 theListOfSecondaries[theNumberOfSecondaries] = aTrack;
61 } 61 }
62 else 62 else
63 { 63 {
64 theListOfSecondaries.push_back(aTrack); 64 theListOfSecondaries.push_back(aTrack);
65 ++theSizeOftheListOfSecondaries; 65 ++theSizeOftheListOfSecondaries;
66 } 66 }
67 ++theNumberOfSecondaries; 67 ++theNumberOfSecondaries;
68 } 68 }
69 69
70 // ------------------------------------------- 70 // --------------------------------------------------------------------
71 G4Step* G4VParticleChange::UpdateStepInfo(G4St 71 G4Step* G4VParticleChange::UpdateStepInfo(G4Step* pStep)
72 { 72 {
73 // Update the G4Step specific attributes 73 // Update the G4Step specific attributes
74 pStep->SetStepLength(theTrueStepLength); 74 pStep->SetStepLength(theTrueStepLength);
75 pStep->AddTotalEnergyDeposit(theLocalEnergyD 75 pStep->AddTotalEnergyDeposit(theLocalEnergyDeposit);
76 pStep->AddNonIonizingEnergyDeposit(theNonIon 76 pStep->AddNonIonizingEnergyDeposit(theNonIonizingEnergyDeposit);
77 pStep->SetControlFlag(theSteppingControlFlag 77 pStep->SetControlFlag(theSteppingControlFlag);
78 78
79 if(theFirstStepInVolume) 79 if(theFirstStepInVolume)
80 { 80 {
81 pStep->SetFirstStepFlag(); 81 pStep->SetFirstStepFlag();
82 } 82 }
83 else 83 else
84 { 84 {
85 pStep->ClearFirstStepFlag(); 85 pStep->ClearFirstStepFlag();
86 } 86 }
87 if(theLastStepInVolume) 87 if(theLastStepInVolume)
88 { 88 {
89 pStep->SetLastStepFlag(); 89 pStep->SetLastStepFlag();
90 } 90 }
91 else 91 else
92 { 92 {
93 pStep->ClearLastStepFlag(); 93 pStep->ClearLastStepFlag();
94 } 94 }
95 95
96 return pStep; 96 return pStep;
97 } 97 }
98 98
99 // ------------------------------------------- 99 // --------------------------------------------------------------------
100 G4Step* G4VParticleChange::UpdateStepForAtRest 100 G4Step* G4VParticleChange::UpdateStepForAtRest(G4Step* Step)
101 { 101 {
102 if(isParentWeightProposed) 102 if(isParentWeightProposed)
103 { 103 {
104 Step->GetPostStepPoint()->SetWeight(thePar 104 Step->GetPostStepPoint()->SetWeight(theParentWeight);
105 } 105 }
106 return UpdateStepInfo(Step); 106 return UpdateStepInfo(Step);
107 } 107 }
108 108
109 // ------------------------------------------- 109 // --------------------------------------------------------------------
110 G4Step* G4VParticleChange::UpdateStepForAlongS 110 G4Step* G4VParticleChange::UpdateStepForAlongStep(G4Step* Step)
111 { 111 {
112 if(isParentWeightProposed) 112 if(isParentWeightProposed)
113 { 113 {
114 G4double initialWeight = Step->GetPreStepP 114 G4double initialWeight = Step->GetPreStepPoint()->GetWeight();
115 G4double currentWeight = Step->GetPostStep 115 G4double currentWeight = Step->GetPostStepPoint()->GetWeight();
116 G4double finalWeight = (theParentWeight 116 G4double finalWeight = (theParentWeight / initialWeight) * currentWeight;
117 Step->GetPostStepPoint()->SetWeight(finalW 117 Step->GetPostStepPoint()->SetWeight(finalWeight);
118 } 118 }
119 return UpdateStepInfo(Step); 119 return UpdateStepInfo(Step);
120 } 120 }
121 121
122 // ------------------------------------------- 122 // --------------------------------------------------------------------
123 G4Step* G4VParticleChange::UpdateStepForPostSt 123 G4Step* G4VParticleChange::UpdateStepForPostStep(G4Step* Step)
124 { 124 {
125 if(isParentWeightProposed) 125 if(isParentWeightProposed)
126 { 126 {
127 Step->GetPostStepPoint()->SetWeight(thePar 127 Step->GetPostStepPoint()->SetWeight(theParentWeight);
128 } 128 }
129 return UpdateStepInfo(Step); 129 return UpdateStepInfo(Step);
130 } 130 }
131 131
132 // ------------------------------------------- 132 // --------------------------------------------------------------------
133 void G4VParticleChange::DumpInfo() const 133 void G4VParticleChange::DumpInfo() const
134 { 134 {
135 auto vol = theCurrentTrack->GetVolume(); 135 auto vol = theCurrentTrack->GetVolume();
136 G4String vname = (nullptr == vol) ? G4String << 136 G4String vname = (nullptr == vol) ? "" : vol->GetName();
137 G4long olprc = G4cout.precision(8); 137 G4long olprc = G4cout.precision(8);
138 G4cout << " --------------------------- 138 G4cout << " -----------------------------------------------" << G4endl;
139 G4cout << " G4VParticleChange Informa 139 G4cout << " G4VParticleChange Information " << G4endl;
140 G4cout << " TrackID : " < 140 G4cout << " TrackID : " << theCurrentTrack->GetTrackID()
141 << G4endl; 141 << G4endl;
142 G4cout << " ParentID : " < 142 G4cout << " ParentID : " << theCurrentTrack->GetParentID()
143 << G4endl; 143 << G4endl;
144 G4cout << " Particle : " 144 G4cout << " Particle : "
145 << theCurrentTrack->GetParticleDefinition() 145 << theCurrentTrack->GetParticleDefinition()->GetParticleName()
146 << G4endl; 146 << G4endl;
147 G4cout << " Kinetic energy (MeV): " 147 G4cout << " Kinetic energy (MeV): "
148 << theCurrentTrack->GetKineticEnergy( 148 << theCurrentTrack->GetKineticEnergy() << G4endl;
149 G4cout << " Position (mm) : " 149 G4cout << " Position (mm) : "
150 << theCurrentTrack->GetPosition() << 150 << theCurrentTrack->GetPosition() << G4endl;
151 G4cout << " Direction : " 151 G4cout << " Direction : "
152 << theCurrentTrack->GetMomentumDirect 152 << theCurrentTrack->GetMomentumDirection() << G4endl;
153 G4cout << " PhysicsVolume : " < 153 G4cout << " PhysicsVolume : " << vname << G4endl;
154 G4cout << " Material : " 154 G4cout << " Material : "
155 << theCurrentTrack->GetMaterial()->GetName( 155 << theCurrentTrack->GetMaterial()->GetName() << G4endl;
156 G4cout << " --------------------------- 156 G4cout << " -----------------------------------------------" << G4endl;
157 157
158 G4cout << " # of secondaries : " < 158 G4cout << " # of secondaries : " << std::setw(20)
159 << theNumberOfSecondaries << G4endl; 159 << theNumberOfSecondaries << G4endl;
160 160
161 G4cout << " --------------------------- 161 G4cout << " -----------------------------------------------" << G4endl;
162 162
163 G4cout << " Energy Deposit (MeV): " < 163 G4cout << " Energy Deposit (MeV): " << std::setw(20)
164 << theLocalEnergyDeposit / MeV << G4e 164 << theLocalEnergyDeposit / MeV << G4endl;
165 165
166 G4cout << " NIEL Energy Deposit (MeV): " < 166 G4cout << " NIEL Energy Deposit (MeV): " << std::setw(20)
167 << theNonIonizingEnergyDeposit / MeV 167 << theNonIonizingEnergyDeposit / MeV << G4endl;
168 168
169 G4cout << " Track Status : " < 169 G4cout << " Track Status : " << std::setw(20);
170 if(theStatusChange == fAlive) 170 if(theStatusChange == fAlive)
171 { 171 {
172 G4cout << " Alive"; 172 G4cout << " Alive";
173 } 173 }
174 else if(theStatusChange == fStopButAlive) 174 else if(theStatusChange == fStopButAlive)
175 { 175 {
176 G4cout << " StopButAlive"; 176 G4cout << " StopButAlive";
177 } 177 }
178 else if(theStatusChange == fStopAndKill) 178 else if(theStatusChange == fStopAndKill)
179 { 179 {
180 G4cout << " StopAndKill"; 180 G4cout << " StopAndKill";
181 } 181 }
182 else if(theStatusChange == fKillTrackAndSeco 182 else if(theStatusChange == fKillTrackAndSecondaries)
183 { 183 {
184 G4cout << " KillTrackAndSecondaries"; 184 G4cout << " KillTrackAndSecondaries";
185 } 185 }
186 else if(theStatusChange == fSuspend) 186 else if(theStatusChange == fSuspend)
187 { 187 {
188 G4cout << " Suspend"; 188 G4cout << " Suspend";
189 } 189 }
190 else if(theStatusChange == fPostponeToNextEv 190 else if(theStatusChange == fPostponeToNextEvent)
191 { 191 {
192 G4cout << " PostponeToNextEvent"; 192 G4cout << " PostponeToNextEvent";
193 } 193 }
194 G4cout << G4endl; 194 G4cout << G4endl;
195 G4cout << " TruePathLength (mm) : " < 195 G4cout << " TruePathLength (mm) : " << std::setw(20)
196 << theTrueStepLength / mm << G4endl; 196 << theTrueStepLength / mm << G4endl;
197 G4cout << " Stepping Control : " < 197 G4cout << " Stepping Control : " << std::setw(20)
198 << theSteppingControlFlag << G4endl; 198 << theSteppingControlFlag << G4endl;
199 if(theFirstStepInVolume) 199 if(theFirstStepInVolume)
200 { 200 {
201 G4cout << " First step in volume" << 201 G4cout << " First step in volume" << G4endl;
202 } 202 }
203 if(theLastStepInVolume) 203 if(theLastStepInVolume)
204 { 204 {
205 G4cout << " Last step in volume" << 205 G4cout << " Last step in volume" << G4endl;
206 } 206 }
207 207
208 #ifdef G4VERBOSE 208 #ifdef G4VERBOSE
209 if(nError == maxError) 209 if(nError == maxError)
210 { 210 {
211 G4cout << " ------------------------- 211 G4cout << " -----------------------------------------------" << G4endl;
212 G4cout << " G4VParticleChange warni 212 G4cout << " G4VParticleChange warnings closed " << G4endl;
213 G4cout << " ------------------------- 213 G4cout << " -----------------------------------------------" << G4endl;
214 } 214 }
215 #endif 215 #endif
216 216
217 G4cout.precision(olprc); 217 G4cout.precision(olprc);
218 } 218 }
219 219
220 // ------------------------------------------- 220 // --------------------------------------------------------------------
221 G4bool G4VParticleChange::CheckIt([[maybe_unus 221 G4bool G4VParticleChange::CheckIt([[maybe_unused]] const G4Track& aTrack)
222 { 222 {
223 G4bool isOK = true; 223 G4bool isOK = true;
224 224
225 // Energy deposit should not be negative 225 // Energy deposit should not be negative
226 if(theLocalEnergyDeposit < 0.0) 226 if(theLocalEnergyDeposit < 0.0)
227 { 227 {
228 isOK = false; 228 isOK = false;
229 ++nError; 229 ++nError;
230 #ifdef G4VERBOSE 230 #ifdef G4VERBOSE
231 if(nError < maxError) 231 if(nError < maxError)
232 { 232 {
233 G4cout << " G4VParticleChange::CheckIt 233 G4cout << " G4VParticleChange::CheckIt : ";
234 G4cout << "the energy deposit " << theLo 234 G4cout << "the energy deposit " << theLocalEnergyDeposit/MeV
235 << " MeV is negative !!" << G4endl; 235 << " MeV is negative !!" << G4endl;
236 } 236 }
237 #endif 237 #endif
238 theLocalEnergyDeposit = 0.0; 238 theLocalEnergyDeposit = 0.0;
239 } 239 }
240 240
241 // true path length should not be negative 241 // true path length should not be negative
242 if(theTrueStepLength < 0.0) 242 if(theTrueStepLength < 0.0)
243 { 243 {
244 isOK = false; 244 isOK = false;
245 ++nError; 245 ++nError;
246 #ifdef G4VERBOSE 246 #ifdef G4VERBOSE
247 if(nError < maxError) 247 if(nError < maxError)
248 { 248 {
249 G4cout << " G4VParticleChange::CheckIt 249 G4cout << " G4VParticleChange::CheckIt : ";
250 G4cout << "true path length " << theTrue 250 G4cout << "true path length " << theTrueStepLength/mm
251 << " mm is negative !!" << G4endl; 251 << " mm is negative !!" << G4endl;
252 } 252 }
253 #endif 253 #endif
254 theTrueStepLength = (1.e-12) * mm; 254 theTrueStepLength = (1.e-12) * mm;
255 } 255 }
256 256
257 if(!isOK) 257 if(!isOK)
258 { 258 {
259 if(nError < maxError) 259 if(nError < maxError)
260 { 260 {
261 #ifdef G4VERBOSE 261 #ifdef G4VERBOSE
262 // dump out information of this particle 262 // dump out information of this particle change
263 DumpInfo(); 263 DumpInfo();
264 #endif 264 #endif
265 G4Exception("G4VParticleChange::CheckIt( 265 G4Exception("G4VParticleChange::CheckIt()", "TRACK001", JustWarning,
266 "Step length and/or energy deposit are i 266 "Step length and/or energy deposit are illegal");
267 } 267 }
268 } 268 }
269 return isOK; 269 return isOK;
270 } 270 }
271 271
272 // ------------------------------------------- 272 // --------------------------------------------------------------------
273 G4bool G4VParticleChange::CheckSecondary(G4Tra 273 G4bool G4VParticleChange::CheckSecondary(G4Track& aTrack)
274 { 274 {
275 G4bool isOK = true; 275 G4bool isOK = true;
276 276
277 // MomentumDirection should be unit vector 277 // MomentumDirection should be unit vector
278 G4double ekin = aTrack.GetKineticEnergy(); 278 G4double ekin = aTrack.GetKineticEnergy();
279 auto dir = aTrack.GetMomentumDirection(); 279 auto dir = aTrack.GetMomentumDirection();
280 G4double accuracy = std::abs(dir.mag2() - 1. 280 G4double accuracy = std::abs(dir.mag2() - 1.0);
281 if(accuracy > accuracyForWarning) 281 if(accuracy > accuracyForWarning)
282 { 282 {
283 isOK = false; 283 isOK = false;
284 ++nError; 284 ++nError;
285 #ifdef G4VERBOSE 285 #ifdef G4VERBOSE
286 if(nError < maxError) 286 if(nError < maxError)
287 { 287 {
288 G4String mname = aTrack.GetCreatorModelN 288 G4String mname = aTrack.GetCreatorModelName();
289 G4cout << " G4VParticleChange::CheckSeco 289 G4cout << " G4VParticleChange::CheckSecondary : " << G4endl;
290 G4cout << " the momentum direction " << 290 G4cout << " the momentum direction " << dir
291 << " is not unit vector !!" << G4endl; 291 << " is not unit vector !!" << G4endl;
292 G4cout << " Difference=" << accuracy 292 G4cout << " Difference=" << accuracy
293 << " Ekin(MeV)=" << ekin/MeV 293 << " Ekin(MeV)=" << ekin/MeV
294 << " " << aTrack.GetParticleDefinition 294 << " " << aTrack.GetParticleDefinition()->GetParticleName()
295 << " created by " << mname 295 << " created by " << mname
296 << G4endl; 296 << G4endl;
297 } 297 }
298 #endif 298 #endif
299 aTrack.SetMomentumDirection(dir.unit()); 299 aTrack.SetMomentumDirection(dir.unit());
300 } 300 }
301 301
302 // Kinetic Energy should not be negative 302 // Kinetic Energy should not be negative
303 if(ekin < 0.0) 303 if(ekin < 0.0)
304 { 304 {
305 isOK = false; 305 isOK = false;
306 ++nError; 306 ++nError;
307 #ifdef G4VERBOSE 307 #ifdef G4VERBOSE
308 if(nError < maxError) 308 if(nError < maxError)
309 { 309 {
310 G4String mname = aTrack.GetCreatorModelN 310 G4String mname = aTrack.GetCreatorModelName();
311 G4cout << " G4VParticleChange::CheckSeco 311 G4cout << " G4VParticleChange::CheckSecondary : " << G4endl;
312 G4cout << " Ekin(MeV)=" << ekin << " is 312 G4cout << " Ekin(MeV)=" << ekin << " is negative !! "
313 << aTrack.GetParticleDefinition()->GetP 313 << aTrack.GetParticleDefinition()->GetParticleName()
314 << " created by " << mname 314 << " created by " << mname
315 << G4endl; 315 << G4endl;
316 } 316 }
317 #endif 317 #endif
318 aTrack.SetKineticEnergy(0.0); 318 aTrack.SetKineticEnergy(0.0);
319 } 319 }
320 320
321 // Check timing of secondaries 321 // Check timing of secondaries
322 G4double time = aTrack.GetGlobalTime(); 322 G4double time = aTrack.GetGlobalTime();
323 if(time < theParentGlobalTime) 323 if(time < theParentGlobalTime)
324 { 324 {
325 isOK = false; 325 isOK = false;
326 ++nError; 326 ++nError;
327 #ifdef G4VERBOSE 327 #ifdef G4VERBOSE
328 if(nError < maxError) 328 if(nError < maxError)
329 { 329 {
330 G4String mname = aTrack.GetCreatorModelN 330 G4String mname = aTrack.GetCreatorModelName();
331 G4cout << " G4VParticleChange::CheckSeco 331 G4cout << " G4VParticleChange::CheckSecondary : " << G4endl;
332 G4cout << " The global time of secondary 332 G4cout << " The global time of secondary goes back compared to the parent !!" << G4endl;
333 G4cout << " ParentTime(ns)=" << theParen 333 G4cout << " ParentTime(ns)=" << theParentGlobalTime/ns
334 << " SecondaryTime(ns)= " << time/ns 334 << " SecondaryTime(ns)= " << time/ns
335 << " Difference(ns)=" << (theParentGlob 335 << " Difference(ns)=" << (theParentGlobalTime - time)/ns
336 << G4endl; 336 << G4endl;
337 G4cout << " Ekin(MeV)=" << ekin 337 G4cout << " Ekin(MeV)=" << ekin
338 << aTrack.GetParticleDefinition()->GetP 338 << aTrack.GetParticleDefinition()->GetParticleName()
339 << " created by " << mname << G4endl; 339 << " created by " << mname << G4endl;
340 } 340 }
341 #endif 341 #endif
342 aTrack.SetGlobalTime(theParentGlobalTime); 342 aTrack.SetGlobalTime(theParentGlobalTime);
343 } 343 }
344 344
345 // Exit with error 345 // Exit with error
346 if(!isOK) 346 if(!isOK)
347 { 347 {
348 if(nError < maxError) 348 if(nError < maxError)
349 { 349 {
350 #ifdef G4VERBOSE 350 #ifdef G4VERBOSE
351 DumpInfo(); 351 DumpInfo();
352 #endif 352 #endif
353 G4Exception("G4VParticleChange::CheckSec 353 G4Exception("G4VParticleChange::CheckSecondary()", "TRACK001",
354 JustWarning, "Secondary with illegal tim 354 JustWarning, "Secondary with illegal time and/or energy and/or momentum");
355 } 355 }
356 } 356 }
357 return isOK; 357 return isOK;
358 } 358 }
359 359
360 // ------------------------------------------- 360 // --------------------------------------------------------------------
361 G4double G4VParticleChange::GetAccuracyForWarn 361 G4double G4VParticleChange::GetAccuracyForWarning() const
362 { 362 {
363 return accuracyForWarning; 363 return accuracyForWarning;
364 } 364 }
365 365
366 // ------------------------------------------- 366 // --------------------------------------------------------------------
367 G4double G4VParticleChange::GetAccuracyForExce 367 G4double G4VParticleChange::GetAccuracyForException() const
368 { 368 {
369 return accuracyForException; 369 return accuracyForException;
370 } 370 }
371 371
372 // ------------------------------------------- 372 // --------------------------------------------------------------------
373 // Obsolete methods for parent weight 373 // Obsolete methods for parent weight
374 // 374 //
375 void G4VParticleChange::SetParentWeightByProce 375 void G4VParticleChange::SetParentWeightByProcess(G4bool) {}
376 G4bool G4VParticleChange::IsParentWeightSetByP 376 G4bool G4VParticleChange::IsParentWeightSetByProcess() const { return true; }
377 377