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