| 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 <<
27 // 26 //
28 // Author: Hisaya Kurashige, 23 March 1998 << 27 //
29 // ------------------------------------------- << 28 //
>> 29 // --------------------------------------------------------------
>> 30 // GEANT 4 class implementation file
>> 31 //
>> 32 //
>> 33 // ------------------------------------------------------------
>> 34 // Implemented for the new scheme 23 Mar. 1998 H.Kurahige
>> 35 // --------------------------------------------------------------
30 36
31 #include "G4VParticleChange.hh" 37 #include "G4VParticleChange.hh"
32 #include "G4SystemOfUnits.hh" 38 #include "G4SystemOfUnits.hh"
>> 39 #include "G4Track.hh"
>> 40 #include "G4Step.hh"
>> 41 #include "G4TrackFastVector.hh"
33 #include "G4ExceptionSeverity.hh" 42 #include "G4ExceptionSeverity.hh"
34 43
35 const G4double G4VParticleChange::accuracyForW << 44 const G4double G4VParticleChange::accuracyForWarning = 1.0e-9;
36 const G4double G4VParticleChange::accuracyForE 45 const G4double G4VParticleChange::accuracyForException = 0.001;
37 const G4int G4VParticleChange::maxError = 10; <<
38 46
39 // ------------------------------------------- <<
40 G4VParticleChange::G4VParticleChange() 47 G4VParticleChange::G4VParticleChange()
>> 48 :theListOfSecondaries(0),
>> 49 theNumberOfSecondaries(0),
>> 50 theSizeOftheListOfSecondaries(G4TrackFastVectorSize),
>> 51 theStatusChange(fAlive),
>> 52 theSteppingControlFlag(NormalCondition),
>> 53 theLocalEnergyDeposit(0.0),
>> 54 theNonIonizingEnergyDeposit(0.0),
>> 55 theTrueStepLength(0.0),
>> 56 theFirstStepInVolume(false),
>> 57 theLastStepInVolume(false),
>> 58 theParentWeight(1.0),
>> 59 isParentWeightProposed(false),
>> 60 fSetSecondaryWeightByProcess(false),
>> 61 theParentGlobalTime(0.0),
>> 62 verboseLevel(1),
>> 63 debugFlag(false)
41 { 64 {
42 #ifdef G4VERBOSE 65 #ifdef G4VERBOSE
43 // activate CheckIt if in VERBOSE mode << 66 // activate CHeckIt if in VERBOSE mode
44 debugFlag = true; 67 debugFlag = true;
45 #endif 68 #endif
>> 69 theListOfSecondaries = new G4TrackFastVector();
>> 70 }
>> 71
>> 72 G4VParticleChange::~G4VParticleChange() {
>> 73 // check if tracks still exist in theListOfSecondaries
>> 74 if (theNumberOfSecondaries>0) {
>> 75 #ifdef G4VERBOSE
>> 76 if (verboseLevel>0) {
>> 77 G4cout << "G4VParticleChange::~G4VParticleChange() Warning ";
>> 78 G4cout << "theListOfSecondaries is not empty ";
>> 79 }
>> 80 #endif
>> 81 for (G4int index= 0; index<theNumberOfSecondaries; index++){
>> 82 delete (*theListOfSecondaries)[index] ;
>> 83 }
>> 84 }
>> 85 delete theListOfSecondaries;
>> 86 }
>> 87
>> 88 G4VParticleChange::G4VParticleChange(const G4VParticleChange &right)
>> 89 :theListOfSecondaries(0),
>> 90 theNumberOfSecondaries(0),
>> 91 theSizeOftheListOfSecondaries(G4TrackFastVectorSize),
>> 92 theStatusChange( right.theStatusChange),
>> 93 theSteppingControlFlag(right.theSteppingControlFlag),
>> 94 theLocalEnergyDeposit(right.theLocalEnergyDeposit),
>> 95 theNonIonizingEnergyDeposit(right.theNonIonizingEnergyDeposit),
>> 96 theTrueStepLength(right.theTrueStepLength),
>> 97 theFirstStepInVolume( right.theFirstStepInVolume),
>> 98 theLastStepInVolume(right.theLastStepInVolume),
>> 99 theParentWeight(right.theParentWeight),
>> 100 isParentWeightProposed(false),
>> 101 fSetSecondaryWeightByProcess(right.fSetSecondaryWeightByProcess),
>> 102 theParentGlobalTime(0.0),
>> 103 verboseLevel(right.verboseLevel),
>> 104 debugFlag(right.debugFlag)
>> 105 {
>> 106 theListOfSecondaries = new G4TrackFastVector();
>> 107 theNumberOfSecondaries = right.theNumberOfSecondaries;
>> 108 for (G4int index = 0; index<theNumberOfSecondaries; index++){
>> 109 G4Track* newTrack = new G4Track(*((*right.theListOfSecondaries)[index] ));
>> 110 theListOfSecondaries->SetElement(index, newTrack);
>> 111 }
46 } 112 }
47 113
48 // ------------------------------------------- << 114
49 void G4VParticleChange::AddSecondary(G4Track* << 115 G4VParticleChange & G4VParticleChange::operator=(const G4VParticleChange &right)
50 { 116 {
51 if(debugFlag) << 117 if (this != &right){
52 CheckSecondary(*aTrack); << 118 if (theNumberOfSecondaries>0) {
>> 119 #ifdef G4VERBOSE
>> 120 if (verboseLevel>0) {
>> 121 G4cout << "G4VParticleChange: assignment operator Warning ";
>> 122 G4cout << "theListOfSecondaries is not empty ";
>> 123 }
>> 124 #endif
>> 125 for (G4int index = 0; index<theNumberOfSecondaries; index++){
>> 126 if ( (*theListOfSecondaries)[index] ) delete ((*theListOfSecondaries)[index]) ;
>> 127 }
>> 128 }
>> 129 delete theListOfSecondaries;
>> 130
>> 131 theListOfSecondaries = new G4TrackFastVector();
>> 132 theNumberOfSecondaries = right.theNumberOfSecondaries;
>> 133 for (G4int index = 0; index<theNumberOfSecondaries; index++){
>> 134 G4Track* newTrack = new G4Track(*((*right.theListOfSecondaries)[index] ));
>> 135 theListOfSecondaries->SetElement(index, newTrack);
>> 136 }
>> 137 theStatusChange = right.theStatusChange;
>> 138 theSteppingControlFlag = right.theSteppingControlFlag;
>> 139 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
>> 140 theNonIonizingEnergyDeposit = right.theNonIonizingEnergyDeposit;
>> 141 theTrueStepLength = right.theTrueStepLength;
>> 142
>> 143 theFirstStepInVolume = right.theFirstStepInVolume;
>> 144 theLastStepInVolume = right.theLastStepInVolume;
53 145
54 if(!fSetSecondaryWeightByProcess) << 146 theParentWeight = right.theParentWeight;
55 aTrack->SetWeight(theParentWeight); << 147 isParentWeightProposed = right.isParentWeightProposed;
>> 148 fSetSecondaryWeightByProcess = right.fSetSecondaryWeightByProcess;
56 149
57 // add a secondary after size check << 150 theParentGlobalTime = right.theParentGlobalTime;
58 if(theSizeOftheListOfSecondaries > theNumber << 151
59 { << 152 verboseLevel = right.verboseLevel;
60 theListOfSecondaries[theNumberOfSecondarie << 153 debugFlag = right.debugFlag;
>> 154
61 } 155 }
62 else << 156 return *this;
63 { << 157 }
64 theListOfSecondaries.push_back(aTrack); << 158
65 ++theSizeOftheListOfSecondaries; << 159 G4bool G4VParticleChange::operator==(const G4VParticleChange &right) const
>> 160 {
>> 161 return (this == (G4VParticleChange *) &right);
>> 162 }
>> 163
>> 164
>> 165 G4bool G4VParticleChange::operator!=(const G4VParticleChange &right) const
>> 166 {
>> 167 return (this != (G4VParticleChange *) &right);
>> 168 }
>> 169
>> 170 void G4VParticleChange::AddSecondary(G4Track *aTrack)
>> 171 {
>> 172 if (debugFlag) CheckSecondary(*aTrack);
>> 173
>> 174 // add a secondary after size check
>> 175 if (theSizeOftheListOfSecondaries > theNumberOfSecondaries) {
>> 176 // Set weight of secondary tracks
>> 177 if (!fSetSecondaryWeightByProcess) aTrack->SetWeight(theParentWeight);
>> 178 theListOfSecondaries->SetElement(theNumberOfSecondaries, aTrack);
>> 179 theNumberOfSecondaries++;
>> 180 } else {
>> 181 delete aTrack;
>> 182
>> 183 #ifdef G4VERBOSE
>> 184 if (verboseLevel>0) {
>> 185 G4cout << "G4VParticleChange::AddSecondary() Warning ";
>> 186 G4cout << "theListOfSecondaries is full !! " << G4endl;
>> 187 G4cout << " The track is deleted " << G4endl;
>> 188 }
>> 189 #endif
>> 190 G4Exception("G4VParticleChange::AddSecondary",
>> 191 "TRACK101", JustWarning,
>> 192 "Secondary Bug is full. The track is deleted");
66 } 193 }
67 ++theNumberOfSecondaries; <<
68 } 194 }
69 195
70 // ------------------------------------------- << 196
>> 197
>> 198 // Virtual methods for updating G4Step
>> 199 //
>> 200
71 G4Step* G4VParticleChange::UpdateStepInfo(G4St 201 G4Step* G4VParticleChange::UpdateStepInfo(G4Step* pStep)
72 { 202 {
73 // Update the G4Step specific attributes 203 // Update the G4Step specific attributes
74 pStep->SetStepLength(theTrueStepLength); << 204 pStep->SetStepLength( theTrueStepLength );
75 pStep->AddTotalEnergyDeposit(theLocalEnergyD << 205 pStep->AddTotalEnergyDeposit( theLocalEnergyDeposit );
76 pStep->AddNonIonizingEnergyDeposit(theNonIon << 206 pStep->AddNonIonizingEnergyDeposit( theNonIonizingEnergyDeposit );
77 pStep->SetControlFlag(theSteppingControlFlag << 207 pStep->SetControlFlag( theSteppingControlFlag );
78 << 208
79 if(theFirstStepInVolume) << 209 if (theFirstStepInVolume) {pStep->SetFirstStepFlag();}
80 { << 210 else {pStep->ClearFirstStepFlag();}
81 pStep->SetFirstStepFlag(); << 211 if (theLastStepInVolume) {pStep->SetLastStepFlag();}
82 } << 212 else {pStep->ClearLastStepFlag();}
83 else <<
84 { <<
85 pStep->ClearFirstStepFlag(); <<
86 } <<
87 if(theLastStepInVolume) <<
88 { <<
89 pStep->SetLastStepFlag(); <<
90 } <<
91 else <<
92 { <<
93 pStep->ClearLastStepFlag(); <<
94 } <<
95 213
96 return pStep; 214 return pStep;
97 } 215 }
98 216
99 // ------------------------------------------- << 217
100 G4Step* G4VParticleChange::UpdateStepForAtRest 218 G4Step* G4VParticleChange::UpdateStepForAtRest(G4Step* Step)
101 { << 219 {
102 if(isParentWeightProposed) << 220 if (isParentWeightProposed ){
103 { << 221 Step->GetPostStepPoint()->SetWeight( theParentWeight );
104 Step->GetPostStepPoint()->SetWeight(thePar <<
105 } 222 }
106 return UpdateStepInfo(Step); 223 return UpdateStepInfo(Step);
107 } 224 }
108 225
109 // ------------------------------------------- << 226
110 G4Step* G4VParticleChange::UpdateStepForAlongS 227 G4Step* G4VParticleChange::UpdateStepForAlongStep(G4Step* Step)
111 { 228 {
112 if(isParentWeightProposed) << 229 if (isParentWeightProposed ){
113 { <<
114 G4double initialWeight = Step->GetPreStepP 230 G4double initialWeight = Step->GetPreStepPoint()->GetWeight();
115 G4double currentWeight = Step->GetPostStep 231 G4double currentWeight = Step->GetPostStepPoint()->GetWeight();
116 G4double finalWeight = (theParentWeight << 232 G4double finalWeight = (theParentWeight/initialWeight)*currentWeight;
117 Step->GetPostStepPoint()->SetWeight(finalW << 233 Step->GetPostStepPoint()->SetWeight( finalWeight );
118 } << 234 }
119 return UpdateStepInfo(Step); 235 return UpdateStepInfo(Step);
120 } 236 }
121 237
122 // ------------------------------------------- <<
123 G4Step* G4VParticleChange::UpdateStepForPostSt 238 G4Step* G4VParticleChange::UpdateStepForPostStep(G4Step* Step)
124 { 239 {
125 if(isParentWeightProposed) << 240 if (isParentWeightProposed ){
126 { << 241 Step->GetPostStepPoint()->SetWeight( theParentWeight );
127 Step->GetPostStepPoint()->SetWeight(thePar <<
128 } 242 }
129 return UpdateStepInfo(Step); 243 return UpdateStepInfo(Step);
130 } 244 }
131 245
132 // ------------------------------------------- << 246
>> 247 //----------------------------------------------------------------
>> 248 // methods for printing messages
>> 249 //
>> 250
133 void G4VParticleChange::DumpInfo() const 251 void G4VParticleChange::DumpInfo() const
134 { 252 {
135 auto vol = theCurrentTrack->GetVolume(); << 253
136 G4String vname = (nullptr == vol) ? G4String << 254 // Show header
137 G4long olprc = G4cout.precision(8); << 255 G4int olprc = G4cout.precision(3);
138 G4cout << " --------------------------- << 256 G4cout << " -----------------------------------------------"
139 G4cout << " G4VParticleChange Informa << 257 << G4endl;
140 G4cout << " TrackID : " < << 258 G4cout << " G4ParticleChange Information " << std::setw(20) << G4endl;
141 << G4endl; << 259 G4cout << " -----------------------------------------------"
142 G4cout << " ParentID : " < << 260 << G4endl;
143 << G4endl; << 261
144 G4cout << " Particle : " << 262 G4cout << " # of 2ndaries : "
145 << theCurrentTrack->GetParticleDefinition() << 263 << std::setw(20) << theNumberOfSecondaries
>> 264 << G4endl;
>> 265
>> 266 if (theNumberOfSecondaries >0) {
>> 267 G4cout << " Pointer to 2ndaries : "
>> 268 << std::setw(20) << GetSecondary(0)
>> 269 << G4endl;
>> 270 G4cout << " (Showed only 1st one)"
146 << G4endl; 271 << 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 } <<
194 G4cout << G4endl; <<
195 G4cout << " TruePathLength (mm) : " < <<
196 << theTrueStepLength / mm << G4endl; <<
197 G4cout << " Stepping Control : " < <<
198 << theSteppingControlFlag << G4endl; <<
199 if(theFirstStepInVolume) <<
200 { <<
201 G4cout << " First step in volume" << <<
202 } <<
203 if(theLastStepInVolume) <<
204 { <<
205 G4cout << " Last step in volume" << <<
206 } <<
207 <<
208 #ifdef G4VERBOSE <<
209 if(nError == maxError) <<
210 { <<
211 G4cout << " ------------------------- <<
212 G4cout << " G4VParticleChange warni <<
213 G4cout << " ------------------------- <<
214 } 272 }
215 #endif << 273 G4cout << " -----------------------------------------------"
216 << 274 << G4endl;
>> 275
>> 276 G4cout << " Energy Deposit (MeV): "
>> 277 << std::setw(20) << theLocalEnergyDeposit/MeV
>> 278 << G4endl;
>> 279
>> 280 G4cout << " Non-ionizing Energy Deposit (MeV): "
>> 281 << std::setw(20) << theNonIonizingEnergyDeposit/MeV
>> 282 << G4endl;
>> 283
>> 284 G4cout << " Track Status : "
>> 285 << std::setw(20);
>> 286 if( theStatusChange == fAlive ){
>> 287 G4cout << " Alive";
>> 288 } else if( theStatusChange == fStopButAlive ){
>> 289 G4cout << " StopButAlive";
>> 290 } else if( theStatusChange == fStopAndKill ){
>> 291 G4cout << " StopAndKill";
>> 292 } else if( theStatusChange == fKillTrackAndSecondaries ){
>> 293 G4cout << " KillTrackAndSecondaries";
>> 294 } else if( theStatusChange == fSuspend ){
>> 295 G4cout << " Suspend";
>> 296 } else if( theStatusChange == fPostponeToNextEvent ){
>> 297 G4cout << " PostponeToNextEvent";
>> 298 }
>> 299 G4cout << G4endl;
>> 300 G4cout << " True Path Length (mm) : "
>> 301 << std::setw(20) << theTrueStepLength/mm
>> 302 << G4endl;
>> 303 G4cout << " Stepping Control : "
>> 304 << std::setw(20) << theSteppingControlFlag
>> 305 << G4endl;
>> 306 if (theFirstStepInVolume) {
>> 307 G4cout << " First Step In the voulme : " << G4endl;
>> 308 }
>> 309 if (theLastStepInVolume) {
>> 310 G4cout << " Last Step In the voulme : " << G4endl;
>> 311 }
217 G4cout.precision(olprc); 312 G4cout.precision(olprc);
218 } 313 }
219 314
220 // ------------------------------------------- << 315 G4bool G4VParticleChange::CheckIt(const G4Track&
221 G4bool G4VParticleChange::CheckIt([[maybe_unus << 316 #ifdef G4VERBOSE
>> 317 aTrack
>> 318 #endif
>> 319 )
222 { 320 {
223 G4bool isOK = true; << 321
>> 322 G4bool exitWithError = false;
>> 323 G4double accuracy;
>> 324 static G4ThreadLocal G4int nError = 0;
>> 325 #ifdef G4VERBOSE
>> 326 const G4int maxError = 30;
>> 327 #endif
224 328
225 // Energy deposit should not be negative 329 // Energy deposit should not be negative
226 if(theLocalEnergyDeposit < 0.0) << 330 G4bool itsOKforEnergy = true;
227 { << 331 accuracy = -1.0*theLocalEnergyDeposit/MeV;
228 isOK = false; << 332 if (accuracy > accuracyForWarning) {
229 ++nError; << 333 itsOKforEnergy = false;
230 #ifdef G4VERBOSE << 334 nError += 1;
231 if(nError < maxError) << 335 exitWithError = (accuracy > accuracyForException);
232 { << 336 #ifdef G4VERBOSE
233 G4cout << " G4VParticleChange::CheckIt << 337 if (nError < maxError) {
234 G4cout << "the energy deposit " << theLo << 338 G4cout << " G4VParticleChange::CheckIt : ";
235 << " MeV is negative !!" << G4endl; << 339 G4cout << "the energy deposit is negative !!"
>> 340 << " Difference: " << accuracy << "[MeV] " <<G4endl;
>> 341 G4cout << aTrack.GetDefinition()->GetParticleName()
>> 342 << " E=" << aTrack.GetKineticEnergy()/MeV
>> 343 << " pos=" << aTrack.GetPosition().x()/m
>> 344 << ", " << aTrack.GetPosition().y()/m
>> 345 << ", " << aTrack.GetPosition().z()/m
>> 346 <<G4endl;
236 } 347 }
237 #endif 348 #endif
238 theLocalEnergyDeposit = 0.0; <<
239 } 349 }
240 << 350
241 // true path length should not be negative 351 // true path length should not be negative
242 if(theTrueStepLength < 0.0) << 352 G4bool itsOKforStepLength = true;
243 { << 353 accuracy = -1.0*theTrueStepLength/mm;
244 isOK = false; << 354 if (accuracy > accuracyForWarning) {
245 ++nError; << 355 itsOKforStepLength = false;
246 #ifdef G4VERBOSE << 356 nError += 1;
247 if(nError < maxError) << 357 exitWithError = (accuracy > accuracyForException);
248 { << 358 #ifdef G4VERBOSE
249 G4cout << " G4VParticleChange::CheckIt << 359 if (nError < maxError) {
250 G4cout << "true path length " << theTrue << 360 G4cout << " G4VParticleChange::CheckIt : ";
251 << " mm is negative !!" << G4endl; << 361 G4cout << "the true step length is negative !!"
>> 362 << " Difference: " << accuracy << "[MeV] " <<G4endl;
>> 363 G4cout << aTrack.GetDefinition()->GetParticleName()
>> 364 << " E=" << aTrack.GetKineticEnergy()/MeV
>> 365 << " pos=" << aTrack.GetPosition().x()/m
>> 366 << ", " << aTrack.GetPosition().y()/m
>> 367 << ", " << aTrack.GetPosition().z()/m
>> 368 <<G4endl;
252 } 369 }
253 #endif 370 #endif
254 theTrueStepLength = (1.e-12) * mm; <<
255 } <<
256 371
257 if(!isOK) << 372 }
258 { <<
259 if(nError < maxError) <<
260 { <<
261 #ifdef G4VERBOSE 373 #ifdef G4VERBOSE
262 // dump out information of this particle << 374 if (!itsOKforStepLength || !itsOKforEnergy ){
263 DumpInfo(); << 375 // dump out information of this particle change
>> 376 DumpInfo();
>> 377 }
264 #endif 378 #endif
265 G4Exception("G4VParticleChange::CheckIt( << 379
266 "Step length and/or energy deposit are i << 380 // Exit with error
267 } << 381 if (exitWithError) {
>> 382 G4Exception("G4VParticleChange::CheckIt",
>> 383 "TRACK001", EventMustBeAborted,
>> 384 "Step length and/or energy deposit was illegal");
268 } 385 }
269 return isOK; << 386
>> 387 // correction
>> 388 if ( !itsOKforStepLength ) {
>> 389 theTrueStepLength = (1.e-12)*mm;
>> 390 }
>> 391 if ( !itsOKforEnergy ) {
>> 392 theLocalEnergyDeposit = 0.0;
>> 393 }
>> 394 return (itsOKforStepLength && itsOKforEnergy );
270 } 395 }
271 396
272 // ------------------------------------------- <<
273 G4bool G4VParticleChange::CheckSecondary(G4Tra 397 G4bool G4VParticleChange::CheckSecondary(G4Track& aTrack)
274 { 398 {
275 G4bool isOK = true; << 399 G4bool exitWithError = false;
>> 400 G4double accuracy;
>> 401 static G4ThreadLocal G4int nError = 0;
>> 402 #ifdef G4VERBOSE
>> 403 const G4int maxError = 30;
>> 404 #endif
276 405
277 // MomentumDirection should be unit vector 406 // MomentumDirection should be unit vector
278 G4double ekin = aTrack.GetKineticEnergy(); << 407 G4bool itsOKforMomentum = true;
279 auto dir = aTrack.GetMomentumDirection(); << 408 if (aTrack.GetKineticEnergy()>0.){
280 G4double accuracy = std::abs(dir.mag2() - 1. << 409 accuracy = std::fabs((aTrack.GetMomentumDirection()).mag2()-1.0);
281 if(accuracy > accuracyForWarning) << 410 if (accuracy > accuracyForWarning) {
282 { << 411 itsOKforMomentum = false;
283 isOK = false; << 412 nError += 1;
284 ++nError; << 413 exitWithError = exitWithError || (accuracy > accuracyForException);
285 #ifdef G4VERBOSE << 414 #ifdef G4VERBOSE
286 if(nError < maxError) << 415 if (nError < maxError) {
287 { << 416 G4cout << " G4VParticleChange::CheckSecondary : ";
288 G4String mname = aTrack.GetCreatorModelN << 417 G4cout << "the Momentum direction is not unit vector !! "
289 G4cout << " G4VParticleChange::CheckSeco << 418 << " Difference: " << accuracy << G4endl;
290 G4cout << " the momentum direction " << << 419 G4cout << aTrack.GetDefinition()->GetParticleName()
291 << " is not unit vector !!" << G4endl; << 420 << " E=" << aTrack.GetKineticEnergy()/MeV
292 G4cout << " Difference=" << accuracy << 421 << " pos=" << aTrack.GetPosition().x()/m
293 << " Ekin(MeV)=" << ekin/MeV << 422 << ", " << aTrack.GetPosition().y()/m
294 << " " << aTrack.GetParticleDefinition << 423 << ", " << aTrack.GetPosition().z()/m
295 << " created by " << mname << 424 <<G4endl;
296 << G4endl; << 425 }
297 } <<
298 #endif 426 #endif
299 aTrack.SetMomentumDirection(dir.unit()); << 427 }
300 } 428 }
301 << 429
302 // Kinetic Energy should not be negative 430 // Kinetic Energy should not be negative
303 if(ekin < 0.0) << 431 G4bool itsOKforEnergy = true;
304 { << 432 accuracy = -1.0*(aTrack.GetKineticEnergy())/MeV;
305 isOK = false; << 433 if (accuracy > accuracyForWarning) {
306 ++nError; << 434 itsOKforEnergy = false;
307 #ifdef G4VERBOSE << 435 nError += 1;
308 if(nError < maxError) << 436 exitWithError = exitWithError || (accuracy > accuracyForException);
309 { << 437 #ifdef G4VERBOSE
310 G4String mname = aTrack.GetCreatorModelN << 438 if (nError < maxError) {
311 G4cout << " G4VParticleChange::CheckSeco << 439 G4cout << " G4VParticleChange::CheckSecondary : ";
312 G4cout << " Ekin(MeV)=" << ekin << " is << 440 G4cout << "the kinetic energy is negative !!"
313 << aTrack.GetParticleDefinition()->GetP << 441 << " Difference: " << accuracy << "[MeV] " <<G4endl;
314 << " created by " << mname << 442 G4cout << " G4VParticleChange::CheckSecondary : ";
315 << G4endl; << 443 G4cout << "the global time of secondary is earlier than the parent !!"
>> 444 << " Difference: " << accuracy << "[ns] " <<G4endl;
>> 445 G4cout << aTrack.GetDefinition()->GetParticleName()
>> 446 << " E=" << aTrack.GetKineticEnergy()/MeV
>> 447 << " pos=" << aTrack.GetPosition().x()/m
>> 448 << ", " << aTrack.GetPosition().y()/m
>> 449 << ", " << aTrack.GetPosition().z()/m
>> 450 <<G4endl;
316 } 451 }
317 #endif 452 #endif
318 aTrack.SetKineticEnergy(0.0); <<
319 } 453 }
320 <<
321 // Check timing of secondaries 454 // Check timing of secondaries
322 G4double time = aTrack.GetGlobalTime(); << 455 G4bool itsOKforTiming = true;
323 if(time < theParentGlobalTime) << 456
324 { << 457 accuracy = (theParentGlobalTime - aTrack.GetGlobalTime())/ns;
325 isOK = false; << 458 if (accuracy > accuracyForWarning){
326 ++nError; << 459 itsOKforTiming = false;
327 #ifdef G4VERBOSE << 460 nError += 1;
328 if(nError < maxError) << 461 exitWithError = (accuracy > accuracyForException);
329 { << 462 #ifdef G4VERBOSE
330 G4String mname = aTrack.GetCreatorModelN << 463 if (nError < maxError) {
331 G4cout << " G4VParticleChange::CheckSeco << 464 G4cout << " G4VParticleChange::CheckSecondary : ";
332 G4cout << " The global time of secondary << 465 G4cout << "the global time of secondary goes back comapared to the parent !!"
333 G4cout << " ParentTime(ns)=" << theParen << 466 << " Difference: " << accuracy << "[ns] " <<G4endl;
334 << " SecondaryTime(ns)= " << time/ns << 467 G4cout << aTrack.GetDefinition()->GetParticleName()
335 << " Difference(ns)=" << (theParentGlob << 468 << " E=" << aTrack.GetKineticEnergy()/MeV
336 << G4endl; << 469 << " pos=" << aTrack.GetPosition().x()/m
337 G4cout << " Ekin(MeV)=" << ekin << 470 << ", " << aTrack.GetPosition().y()/m
338 << aTrack.GetParticleDefinition()->GetP << 471 << ", " << aTrack.GetPosition().z()/m
339 << " created by " << mname << G4endl; << 472 << " time=" << aTrack.GetGlobalTime()/ns
>> 473 << " parent time=" << theParentGlobalTime/ns
>> 474 <<G4endl;
340 } 475 }
341 #endif 476 #endif
342 aTrack.SetGlobalTime(theParentGlobalTime); <<
343 } 477 }
344 478
345 // Exit with error 479 // Exit with error
346 if(!isOK) << 480 if (exitWithError) {
347 { << 481 G4Exception("G4VParticleChange::CheckSecondary",
348 if(nError < maxError) << 482 "TRACK001", EventMustBeAborted,
349 { << 483 "Secondary with illegal energy/momentum ");
350 #ifdef G4VERBOSE << 484 }
351 DumpInfo(); << 485
352 #endif << 486 G4bool itsOK = itsOKforMomentum && itsOKforEnergy && itsOKforTiming;
353 G4Exception("G4VParticleChange::CheckSec << 487
354 JustWarning, "Secondary with illegal tim << 488 //correction
355 } << 489 if (!itsOKforMomentum) {
>> 490 G4double vmag = (aTrack.GetMomentumDirection()).mag();
>> 491 aTrack.SetMomentumDirection((1./vmag)*aTrack.GetMomentumDirection());
>> 492 }
>> 493 if (!itsOKforEnergy) {
>> 494 aTrack.SetKineticEnergy(0.0);
356 } 495 }
357 return isOK; << 496
>> 497 if (!itsOK) {
>> 498 this->DumpInfo();
>> 499
>> 500 }
>> 501
>> 502
>> 503 return itsOK;
358 } 504 }
359 505
360 // ------------------------------------------- << 506
361 G4double G4VParticleChange::GetAccuracyForWarn 507 G4double G4VParticleChange::GetAccuracyForWarning() const
362 { 508 {
363 return accuracyForWarning; 509 return accuracyForWarning;
364 } 510 }
365 511
366 // ------------------------------------------- <<
367 G4double G4VParticleChange::GetAccuracyForExce 512 G4double G4VParticleChange::GetAccuracyForException() const
368 { 513 {
369 return accuracyForException; 514 return accuracyForException;
370 } 515 }
371 516
372 // ------------------------------------------- << 517
373 // Obsolete methods for parent weight << 518 ///////////////////////////////////////////////////////////
374 // << 519 //Obsolete methods for parent weight
375 void G4VParticleChange::SetParentWeightByProce << 520 /////////////////////
376 G4bool G4VParticleChange::IsParentWeightSetByP << 521 void G4VParticleChange::SetParentWeightByProcess(G4bool )
>> 522 {
>> 523 }
>> 524
>> 525
>> 526 G4bool G4VParticleChange::IsParentWeightSetByProcess() const
>> 527 {
>> 528 return true;
>> 529 }
377 530