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