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