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