| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
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 // 7 //
26 // G4VParticleChange class implementation << 8 // $Id: G4VParticleChange.cc,v 1.5 2000/02/13 15:08:28 kurasige Exp $
>> 9 // GEANT4 tag $Name: geant4-03-00 $
27 // 10 //
28 // Author: Hisaya Kurashige, 23 March 1998 << 11 //
29 // ------------------------------------------- << 12 // --------------------------------------------------------------
>> 13 // GEANT 4 class implementation file
>> 14 //
>> 15 // For information related to this code contact:
>> 16 // CERN, CN Division, ASD Group
>> 17 //
>> 18 // ------------------------------------------------------------
>> 19 // Implemented for the new scheme 23 Mar. 1998 H.Kurahige
>> 20 // --------------------------------------------------------------
30 21
31 #include "G4VParticleChange.hh" 22 #include "G4VParticleChange.hh"
32 #include "G4SystemOfUnits.hh" << 23 #include "G4Track.hh"
33 #include "G4ExceptionSeverity.hh" << 24 #include "G4Step.hh"
>> 25 #include "G4TrackFastVector.hh"
>> 26 #include "G4Mars5GeVMechanism.hh"
34 27
35 const G4double G4VParticleChange::accuracyForW << 28 const G4double G4VParticleChange::accuracyForWarning = 1.0e-9;
36 const G4double G4VParticleChange::accuracyForE 29 const G4double G4VParticleChange::accuracyForException = 0.001;
37 const G4int G4VParticleChange::maxError = 10; <<
38 30
39 // ------------------------------------------- << 31 G4VParticleChange::G4VParticleChange():
40 G4VParticleChange::G4VParticleChange() << 32 theNumberOfSecondaries(0),
>> 33 theSizeOftheListOfSecondaries(G4TrackFastVectorSize),
>> 34 theStatusChange(fAlive),
>> 35 theSteppingControlFlag(NormalCondition),
>> 36 theLocalEnergyDeposit(0.0),
>> 37 theParentWeight(1.0),
>> 38 theEBMechanism(0),
>> 39 fUseEB(false),
>> 40 verboseLevel(1)
41 { 41 {
>> 42 debugFlag = false;
42 #ifdef G4VERBOSE 43 #ifdef G4VERBOSE
43 // activate CheckIt if in VERBOSE mode << 44 // activate CHeckIt if in VERBOSE mode
44 debugFlag = true; 45 debugFlag = true;
45 #endif 46 #endif
>> 47 theListOfSecondaries = new G4TrackFastVector();
46 } 48 }
47 49
48 // ------------------------------------------- << 50 G4VParticleChange::G4VParticleChange(G4bool useEB):
49 void G4VParticleChange::AddSecondary(G4Track* << 51 theNumberOfSecondaries(0),
50 { << 52 theSizeOftheListOfSecondaries(G4TrackFastVectorSize),
51 if(debugFlag) << 53 theStatusChange(fAlive),
52 CheckSecondary(*aTrack); << 54 theSteppingControlFlag(NormalCondition),
53 << 55 theLocalEnergyDeposit(0.0),
54 if(!fSetSecondaryWeightByProcess) << 56 theParentWeight(1.0),
55 aTrack->SetWeight(theParentWeight); << 57 verboseLevel(1)
>> 58 {
>> 59 fUseEB = useEB;
>> 60 // debug flag (activate CheckIt() )
>> 61 debugFlag = false;
>> 62 #ifdef G4VERBOSE
>> 63 // activate CHeckIt if in VERBOSE mode
>> 64 debugFlag = true;
>> 65 #endif
>> 66 theListOfSecondaries = new G4TrackFastVector();
>> 67 // register G4EvtBiasMechanism as a default
>> 68 theEBMechanism = new G4Mars5GeVMechanism();
>> 69 }
56 70
57 // add a secondary after size check << 71 G4VParticleChange::~G4VParticleChange() {
58 if(theSizeOftheListOfSecondaries > theNumber << 72 // check if tracks still exist in theListOfSecondaries
59 { << 73 if (theNumberOfSecondaries>0) {
60 theListOfSecondaries[theNumberOfSecondarie << 74 #ifdef G4VERBOSE
61 } << 75 if (verboseLevel>0) {
62 else << 76 G4cerr << "G4VParticleChange::~G4VParticleChange() Warning ";
63 { << 77 G4cerr << "theListOfSecondaries is not empty ";
64 theListOfSecondaries.push_back(aTrack); << 78 }
65 ++theSizeOftheListOfSecondaries; << 79 #endif
>> 80 for (G4int index= 0; index<theNumberOfSecondaries; index++){
>> 81 if ( (*theListOfSecondaries)[index] ) delete (*theListOfSecondaries)[index] ;
>> 82 }
66 } 83 }
67 ++theNumberOfSecondaries; << 84 if (theEBMechanism !=0) delete theEBMechanism;
>> 85 delete theListOfSecondaries;
68 } 86 }
69 87
70 // ------------------------------------------- << 88 // copy and assignment operators are implemented as "shallow copy"
71 G4Step* G4VParticleChange::UpdateStepInfo(G4St << 89 G4VParticleChange::G4VParticleChange(const G4VParticleChange &right):
>> 90 theNumberOfSecondaries(0),
>> 91 theSizeOftheListOfSecondaries(G4TrackFastVectorSize),
>> 92 theStatusChange(fAlive),
>> 93 theSteppingControlFlag(NormalCondition),
>> 94 theLocalEnergyDeposit(0.0),
>> 95 theParentWeight(1.0),
>> 96 fUseEB(false),
>> 97 verboseLevel(1)
72 { 98 {
73 // Update the G4Step specific attributes << 99 debugFlag = false;
74 pStep->SetStepLength(theTrueStepLength); << 100 #ifdef G4VERBOSE
75 pStep->AddTotalEnergyDeposit(theLocalEnergyD << 101 // activate CHeckIt if in VERBOSE mode
76 pStep->AddNonIonizingEnergyDeposit(theNonIon << 102 debugFlag = true;
77 pStep->SetControlFlag(theSteppingControlFlag << 103 #endif
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 104
96 return pStep; << 105 theListOfSecondaries = right.theListOfSecondaries;
>> 106 theSizeOftheListOfSecondaries = right.theSizeOftheListOfSecondaries;
>> 107 theNumberOfSecondaries = right.theNumberOfSecondaries;
>> 108 theStatusChange = right.theStatusChange;
>> 109 theTrueStepLength = right.theTrueStepLength;
>> 110 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
>> 111 theSteppingControlFlag = right.theSteppingControlFlag;
97 } 112 }
98 113
99 // ------------------------------------------- << 114
100 G4Step* G4VParticleChange::UpdateStepForAtRest << 115 G4VParticleChange & G4VParticleChange::operator=(const G4VParticleChange &right)
101 { 116 {
102 if(isParentWeightProposed) << 117 debugFlag = false;
103 { << 118 #ifdef G4VERBOSE
104 Step->GetPostStepPoint()->SetWeight(thePar << 119 // activate CHeckIt if in VERBOSE mode
105 } << 120 debugFlag = true;
106 return UpdateStepInfo(Step); << 121 #endif
>> 122 if (this != &right)
>> 123 {
>> 124 theListOfSecondaries = right.theListOfSecondaries;
>> 125 theSizeOftheListOfSecondaries = right.theSizeOftheListOfSecondaries;
>> 126 theNumberOfSecondaries = right.theNumberOfSecondaries;
>> 127 theStatusChange = right.theStatusChange;
>> 128 theTrueStepLength = right.theTrueStepLength;
>> 129 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
>> 130 theSteppingControlFlag = right.theSteppingControlFlag;
>> 131 }
>> 132 return *this;
107 } 133 }
108 134
109 // ------------------------------------------- << 135 G4bool G4VParticleChange::operator==(const G4VParticleChange &right) const
110 G4Step* G4VParticleChange::UpdateStepForAlongS <<
111 { 136 {
112 if(isParentWeightProposed) << 137 return (this == (G4VParticleChange *) &right);
113 { <<
114 G4double initialWeight = Step->GetPreStepP <<
115 G4double currentWeight = Step->GetPostStep <<
116 G4double finalWeight = (theParentWeight <<
117 Step->GetPostStepPoint()->SetWeight(finalW <<
118 } <<
119 return UpdateStepInfo(Step); <<
120 } 138 }
121 139
122 // ------------------------------------------- << 140
123 G4Step* G4VParticleChange::UpdateStepForPostSt << 141 G4bool G4VParticleChange::operator!=(const G4VParticleChange &right) const
124 { 142 {
125 if(isParentWeightProposed) << 143 return (this != (G4VParticleChange *) &right);
126 { <<
127 Step->GetPostStepPoint()->SetWeight(thePar <<
128 } <<
129 return UpdateStepInfo(Step); <<
130 } 144 }
131 145
132 // ------------------------------------------- << 146 //----------------------------------------------------------------
>> 147 // methods for printing messages
>> 148 //
>> 149
133 void G4VParticleChange::DumpInfo() const 150 void G4VParticleChange::DumpInfo() const
134 { 151 {
135 auto vol = theCurrentTrack->GetVolume(); << 152
136 G4String vname = (nullptr == vol) ? G4String << 153 // Show header
137 G4long olprc = G4cout.precision(8); << 154 G4cout.precision(3);
138 G4cout << " --------------------------- << 155 G4cout << " -----------------------------------------------"
139 G4cout << " G4VParticleChange Informa << 156 << G4endl;
140 G4cout << " TrackID : " < << 157 G4cout << " G4ParticleChange Information " << G4std::setw(20) << G4endl;
141 << G4endl; << 158 G4cout << " -----------------------------------------------"
142 G4cout << " ParentID : " < << 159 << G4endl;
143 << G4endl; << 160
144 G4cout << " Particle : " << 161 G4cout << " # of 2ndaries : "
145 << theCurrentTrack->GetParticleDefinition() << 162 << G4std::setw(20) << theNumberOfSecondaries
>> 163 << G4endl;
>> 164
>> 165 if (theNumberOfSecondaries >0) {
>> 166 G4cout << " Pointer to 2ndaries : "
>> 167 << G4std::setw(20) << GetSecondary(0)
>> 168 << G4endl;
>> 169 G4cout << " (Showed only 1st one)"
146 << G4endl; 170 << G4endl;
147 G4cout << " Kinetic energy (MeV): " <<
148 << theCurrentTrack->GetKineticEnergy( <<
149 G4cout << " Position (mm) : " <<
150 << theCurrentTrack->GetPosition() << <<
151 G4cout << " Direction : " <<
152 << theCurrentTrack->GetMomentumDirect <<
153 G4cout << " PhysicsVolume : " < <<
154 G4cout << " Material : " <<
155 << theCurrentTrack->GetMaterial()->GetName( <<
156 G4cout << " --------------------------- <<
157 <<
158 G4cout << " # of secondaries : " < <<
159 << theNumberOfSecondaries << G4endl; <<
160 <<
161 G4cout << " --------------------------- <<
162 <<
163 G4cout << " Energy Deposit (MeV): " < <<
164 << theLocalEnergyDeposit / MeV << G4e <<
165 <<
166 G4cout << " NIEL Energy Deposit (MeV): " < <<
167 << theNonIonizingEnergyDeposit / MeV <<
168 <<
169 G4cout << " Track Status : " < <<
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 == fKillTrackAndSeco <<
183 { <<
184 G4cout << " KillTrackAndSecondaries"; <<
185 } <<
186 else if(theStatusChange == fSuspend) <<
187 { <<
188 G4cout << " Suspend"; <<
189 } <<
190 else if(theStatusChange == fPostponeToNextEv <<
191 { <<
192 G4cout << " PostponeToNextEvent"; <<
193 } 171 }
194 G4cout << G4endl; << 172 G4cout << " -----------------------------------------------"
195 G4cout << " TruePathLength (mm) : " < << 173 << G4endl;
196 << theTrueStepLength / mm << G4endl; << 174
197 G4cout << " Stepping Control : " < << 175 G4cout << " Energy Deposit (MeV): "
198 << theSteppingControlFlag << G4endl; << 176 << G4std::setw(20) << theLocalEnergyDeposit/MeV
199 if(theFirstStepInVolume) << 177 << G4endl;
200 { << 178
201 G4cout << " First step in volume" << << 179 G4cout << " Track Status : "
>> 180 << G4std::setw(20);
>> 181 if( theStatusChange == fAlive ){
>> 182 G4cout << " Alive";
>> 183 } else if( theStatusChange == fStopButAlive ){
>> 184 G4cout << " StopButAlive";
>> 185 } else if( theStatusChange == fStopAndKill ){
>> 186 G4cout << " StopAndKill";
>> 187 } else if( theStatusChange == fKillTrackAndSecondaries ){
>> 188 G4cout << " KillTrackAndSecondaries";
>> 189 } else if( theStatusChange == fSuspend ){
>> 190 G4cout << " Suspend";
>> 191 } else if( theStatusChange == fPostponeToNextEvent ){
>> 192 G4cout << " PostponeToNextEvent";
>> 193 }
>> 194 G4cout << G4endl;
>> 195 G4cout << " True Path Length (mm) : "
>> 196 << G4std::setw(20) << theTrueStepLength/mm
>> 197 << G4endl;
>> 198 G4cout << " Stepping Control : "
>> 199 << G4std::setw(20) << theSteppingControlFlag
>> 200 << G4endl;
>> 201 G4cout << " Event Biasing : ";
>> 202 if (fUseEB) {
>> 203 G4cout << G4std::setw(20) << theEBMechanism->GetName();
>> 204 } else {
>> 205 G4cout << " not used ";
>> 206 }
>> 207 G4cout << G4endl;
>> 208 }
>> 209
>> 210 G4bool G4VParticleChange::CheckIt(const G4Track& aTrack)
>> 211 {
>> 212
>> 213 G4bool exitWithError = false;
>> 214 G4double accuracy;
>> 215 G4double newEnergyDeposit;
>> 216
>> 217 // Energy deposit should not be negative
>> 218 G4bool itsOKforEnergy = true;
>> 219 accuracy = -1.0*theLocalEnergyDeposit/MeV;
>> 220 if (accuracy > accuracyForWarning) {
>> 221 G4cout << " G4VParticleChange::CheckIt : ";
>> 222 G4cout << "the energy deposit is negative !!" << G4endl;
>> 223 G4cout << " Difference: " << accuracy << "[MeV] " <<G4endl;
>> 224 itsOKforEnergy = false;
>> 225 if (accuracy > accuracyForException) exitWithError = true;
202 } 226 }
203 if(theLastStepInVolume) << 227
204 { << 228 // true path length should not be negative
205 G4cout << " Last step in volume" << << 229 G4bool itsOKforStepLength = true;
>> 230 accuracy = -1.0*theTrueStepLength/mm;
>> 231 if (accuracy > accuracyForWarning) {
>> 232 G4cout << " G4VParticleChange::CheckIt : ";
>> 233 G4cout << "the true step length is negative !!" << G4endl;
>> 234 G4cout << " Difference: " << accuracy << "[MeV] " <<G4endl;
>> 235 itsOKforStepLength = false;
>> 236 if (accuracy > accuracyForException) exitWithError = true;
206 } 237 }
207 238
208 #ifdef G4VERBOSE << 239 G4bool itsOK = itsOKforStepLength && itsOKforEnergy ;
209 if(nError == maxError) << 240 // dump out information of this particle change
210 { << 241 if (! itsOK ){
211 G4cout << " ------------------------- << 242 G4cout << " G4VParticleChange::CheckIt " <<G4endl;
212 G4cout << " G4VParticleChange warni << 243 DumpInfo();
213 G4cout << " ------------------------- <<
214 } 244 }
215 #endif <<
216 <<
217 G4cout.precision(olprc); <<
218 } <<
219 245
220 // ------------------------------------------- << 246 // Exit with error
221 G4bool G4VParticleChange::CheckIt([[maybe_unus << 247 if (exitWithError) G4Exception("G4VParticleChange::CheckIt");
222 { <<
223 G4bool isOK = true; <<
224 248
225 // Energy deposit should not be negative << 249 // correction
226 if(theLocalEnergyDeposit < 0.0) << 250 if ( !itsOKforStepLength ) {
227 { << 251 theTrueStepLength = (1.e-12)*mm;
228 isOK = false; << 252 }
229 ++nError; << 253 if ( !itsOKforEnergy ) {
230 #ifdef G4VERBOSE <<
231 if(nError < maxError) <<
232 { <<
233 G4cout << " G4VParticleChange::CheckIt <<
234 G4cout << "the energy deposit " << theLo <<
235 << " MeV is negative !!" << G4endl; <<
236 } <<
237 #endif <<
238 theLocalEnergyDeposit = 0.0; 254 theLocalEnergyDeposit = 0.0;
239 } 255 }
>> 256 return itsOK;
>> 257 }
240 258
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 " << theTrue <<
251 << " mm is negative !!" << G4endl; <<
252 } <<
253 #endif <<
254 theTrueStepLength = (1.e-12) * mm; <<
255 } <<
256 259
257 if(!isOK) <<
258 { <<
259 if(nError < maxError) <<
260 { <<
261 #ifdef G4VERBOSE <<
262 // dump out information of this particle <<
263 DumpInfo(); <<
264 #endif <<
265 G4Exception("G4VParticleChange::CheckIt( <<
266 "Step length and/or energy deposit are i <<
267 } <<
268 } <<
269 return isOK; <<
270 } <<
271 260
272 // ------------------------------------------- <<
273 G4bool G4VParticleChange::CheckSecondary(G4Tra <<
274 { <<
275 G4bool isOK = true; <<
276 261
277 // MomentumDirection should be unit vector <<
278 G4double ekin = aTrack.GetKineticEnergy(); <<
279 auto dir = aTrack.GetMomentumDirection(); <<
280 G4double accuracy = std::abs(dir.mag2() - 1. <<
281 if(accuracy > accuracyForWarning) <<
282 { <<
283 isOK = false; <<
284 ++nError; <<
285 #ifdef G4VERBOSE <<
286 if(nError < maxError) <<
287 { <<
288 G4String mname = aTrack.GetCreatorModelN <<
289 G4cout << " G4VParticleChange::CheckSeco <<
290 G4cout << " the momentum direction " << <<
291 << " is not unit vector !!" << G4endl; <<
292 G4cout << " Difference=" << accuracy <<
293 << " Ekin(MeV)=" << ekin/MeV <<
294 << " " << aTrack.GetParticleDefinition <<
295 << " created by " << mname <<
296 << G4endl; <<
297 } <<
298 #endif <<
299 aTrack.SetMomentumDirection(dir.unit()); <<
300 } <<
301 262
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.GetCreatorModelN <<
311 G4cout << " G4VParticleChange::CheckSeco <<
312 G4cout << " Ekin(MeV)=" << ekin << " is <<
313 << aTrack.GetParticleDefinition()->GetP <<
314 << " created by " << mname <<
315 << G4endl; <<
316 } <<
317 #endif <<
318 aTrack.SetKineticEnergy(0.0); <<
319 } <<
320 263
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.GetCreatorModelN <<
331 G4cout << " G4VParticleChange::CheckSeco <<
332 G4cout << " The global time of secondary <<
333 G4cout << " ParentTime(ns)=" << theParen <<
334 << " SecondaryTime(ns)= " << time/ns <<
335 << " Difference(ns)=" << (theParentGlob <<
336 << G4endl; <<
337 G4cout << " Ekin(MeV)=" << ekin <<
338 << aTrack.GetParticleDefinition()->GetP <<
339 << " created by " << mname << G4endl; <<
340 } <<
341 #endif <<
342 aTrack.SetGlobalTime(theParentGlobalTime); <<
343 } <<
344 264
345 // Exit with error <<
346 if(!isOK) <<
347 { <<
348 if(nError < maxError) <<
349 { <<
350 #ifdef G4VERBOSE <<
351 DumpInfo(); <<
352 #endif <<
353 G4Exception("G4VParticleChange::CheckSec <<
354 JustWarning, "Secondary with illegal tim <<
355 } <<
356 } <<
357 return isOK; <<
358 } <<
359 265
360 // ------------------------------------------- <<
361 G4double G4VParticleChange::GetAccuracyForWarn <<
362 { <<
363 return accuracyForWarning; <<
364 } <<
365 266
366 // ------------------------------------------- <<
367 G4double G4VParticleChange::GetAccuracyForExce <<
368 { <<
369 return accuracyForException; <<
370 } <<
371 267
372 // ------------------------------------------- <<
373 // Obsolete methods for parent weight <<
374 // <<
375 void G4VParticleChange::SetParentWeightByProce <<
376 G4bool G4VParticleChange::IsParentWeightSetByP <<
377 268