00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 #include "G4VParticleChange.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4Track.hh"
00041 #include "G4Step.hh"
00042 #include "G4TrackFastVector.hh"
00043 #include "G4ExceptionSeverity.hh"
00044
00045 const G4double G4VParticleChange::accuracyForWarning = 1.0e-9;
00046 const G4double G4VParticleChange::accuracyForException = 0.001;
00047
00048 G4VParticleChange::G4VParticleChange()
00049 :theListOfSecondaries(0),
00050 theNumberOfSecondaries(0),
00051 theSizeOftheListOfSecondaries(G4TrackFastVectorSize),
00052 theStatusChange(fAlive),
00053 theSteppingControlFlag(NormalCondition),
00054 theLocalEnergyDeposit(0.0),
00055 theNonIonizingEnergyDeposit(0.0),
00056 theTrueStepLength(0.0),
00057 theFirstStepInVolume(false),
00058 theLastStepInVolume(false),
00059 theParentWeight(1.0),
00060 isParentWeightProposed(false),
00061 fSetSecondaryWeightByProcess(false),
00062 theParentGlobalTime(0.0),
00063 verboseLevel(1),
00064 debugFlag(false)
00065 {
00066 #ifdef G4VERBOSE
00067
00068 debugFlag = true;
00069 #endif
00070 theListOfSecondaries = new G4TrackFastVector();
00071 }
00072
00073 G4VParticleChange::~G4VParticleChange() {
00074
00075 if (theNumberOfSecondaries>0) {
00076 #ifdef G4VERBOSE
00077 if (verboseLevel>0) {
00078 G4cout << "G4VParticleChange::~G4VParticleChange() Warning ";
00079 G4cout << "theListOfSecondaries is not empty ";
00080 }
00081 #endif
00082 for (G4int index= 0; index<theNumberOfSecondaries; index++){
00083 delete (*theListOfSecondaries)[index] ;
00084 }
00085 }
00086 delete theListOfSecondaries;
00087 }
00088
00089 G4VParticleChange::G4VParticleChange(const G4VParticleChange &right)
00090 :theListOfSecondaries(0),
00091 theNumberOfSecondaries(0),
00092 theSizeOftheListOfSecondaries(G4TrackFastVectorSize),
00093 theStatusChange( right.theStatusChange),
00094 theSteppingControlFlag(right.theSteppingControlFlag),
00095 theLocalEnergyDeposit(right.theLocalEnergyDeposit),
00096 theNonIonizingEnergyDeposit(right.theNonIonizingEnergyDeposit),
00097 theTrueStepLength(right.theTrueStepLength),
00098 theFirstStepInVolume( right.theFirstStepInVolume),
00099 theLastStepInVolume(right.theLastStepInVolume),
00100 theParentWeight(right.theParentWeight),
00101 isParentWeightProposed(false),
00102 fSetSecondaryWeightByProcess(right.fSetSecondaryWeightByProcess),
00103 theParentGlobalTime(0.0),
00104 verboseLevel(right.verboseLevel),
00105 debugFlag(right.debugFlag)
00106 {
00107 theListOfSecondaries = new G4TrackFastVector();
00108 theNumberOfSecondaries = right.theNumberOfSecondaries;
00109 for (G4int index = 0; index<theNumberOfSecondaries; index++){
00110 G4Track* newTrack = new G4Track(*((*right.theListOfSecondaries)[index] ));
00111 theListOfSecondaries->SetElement(index, newTrack);
00112 }
00113 }
00114
00115
00116 G4VParticleChange & G4VParticleChange::operator=(const G4VParticleChange &right)
00117 {
00118 if (this != &right){
00119 if (theNumberOfSecondaries>0) {
00120 #ifdef G4VERBOSE
00121 if (verboseLevel>0) {
00122 G4cout << "G4VParticleChange: assignment operator Warning ";
00123 G4cout << "theListOfSecondaries is not empty ";
00124 }
00125 #endif
00126 for (G4int index = 0; index<theNumberOfSecondaries; index++){
00127 if ( (*theListOfSecondaries)[index] ) delete ((*theListOfSecondaries)[index]) ;
00128 }
00129 }
00130 delete theListOfSecondaries;
00131
00132 theListOfSecondaries = new G4TrackFastVector();
00133 theNumberOfSecondaries = right.theNumberOfSecondaries;
00134 for (G4int index = 0; index<theNumberOfSecondaries; index++){
00135 G4Track* newTrack = new G4Track(*((*right.theListOfSecondaries)[index] ));
00136 theListOfSecondaries->SetElement(index, newTrack);
00137 }
00138 theStatusChange = right.theStatusChange;
00139 theSteppingControlFlag = right.theSteppingControlFlag;
00140 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
00141 theNonIonizingEnergyDeposit = right.theNonIonizingEnergyDeposit;
00142 theTrueStepLength = right.theTrueStepLength;
00143
00144 theFirstStepInVolume = right.theFirstStepInVolume;
00145 theLastStepInVolume = right.theLastStepInVolume;
00146
00147 theParentWeight = right.theParentWeight;
00148 isParentWeightProposed = right.isParentWeightProposed;
00149 fSetSecondaryWeightByProcess = right.fSetSecondaryWeightByProcess;
00150
00151 theParentGlobalTime = right.theParentGlobalTime;
00152
00153 verboseLevel = right.verboseLevel;
00154 debugFlag = right.debugFlag;
00155
00156 }
00157 return *this;
00158 }
00159
00160 G4bool G4VParticleChange::operator==(const G4VParticleChange &right) const
00161 {
00162 return (this == (G4VParticleChange *) &right);
00163 }
00164
00165
00166 G4bool G4VParticleChange::operator!=(const G4VParticleChange &right) const
00167 {
00168 return (this != (G4VParticleChange *) &right);
00169 }
00170
00171 void G4VParticleChange::AddSecondary(G4Track *aTrack)
00172 {
00173 if (debugFlag) CheckSecondary(*aTrack);
00174
00175
00176 if (theSizeOftheListOfSecondaries > theNumberOfSecondaries) {
00177
00178 if (!fSetSecondaryWeightByProcess) aTrack->SetWeight(theParentWeight);
00179 theListOfSecondaries->SetElement(theNumberOfSecondaries, aTrack);
00180 theNumberOfSecondaries++;
00181 } else {
00182 delete aTrack;
00183
00184 #ifdef G4VERBOSE
00185 if (verboseLevel>0) {
00186 G4cout << "G4VParticleChange::AddSecondary() Warning ";
00187 G4cout << "theListOfSecondaries is full !! " << G4endl;
00188 G4cout << " The track is deleted " << G4endl;
00189 }
00190 #endif
00191 G4Exception("G4VParticleChange::AddSecondary",
00192 "TRACK101", JustWarning,
00193 "Secondary Bug is full. The track is deleted");
00194 }
00195 }
00196
00197
00198
00199
00200
00201
00202 G4Step* G4VParticleChange::UpdateStepInfo(G4Step* pStep)
00203 {
00204
00205 pStep->SetStepLength( theTrueStepLength );
00206 pStep->AddTotalEnergyDeposit( theLocalEnergyDeposit );
00207 pStep->AddNonIonizingEnergyDeposit( theNonIonizingEnergyDeposit );
00208 pStep->SetControlFlag( theSteppingControlFlag );
00209
00210 if (theFirstStepInVolume) {pStep->SetFirstStepFlag();}
00211 else {pStep->ClearFirstStepFlag();}
00212 if (theLastStepInVolume) {pStep->SetLastStepFlag();}
00213 else {pStep->ClearLastStepFlag();}
00214
00215 return pStep;
00216 }
00217
00218
00219 G4Step* G4VParticleChange::UpdateStepForAtRest(G4Step* Step)
00220 {
00221 if (isParentWeightProposed ){
00222 Step->GetPostStepPoint()->SetWeight( theParentWeight );
00223 }
00224 return UpdateStepInfo(Step);
00225 }
00226
00227
00228 G4Step* G4VParticleChange::UpdateStepForAlongStep(G4Step* Step)
00229 {
00230 if (isParentWeightProposed ){
00231 G4double initialWeight = Step->GetPreStepPoint()->GetWeight();
00232 G4double currentWeight = Step->GetPostStepPoint()->GetWeight();
00233 G4double finalWeight = (theParentWeight/initialWeight)*currentWeight;
00234 Step->GetPostStepPoint()->SetWeight( finalWeight );
00235 }
00236 return UpdateStepInfo(Step);
00237 }
00238
00239 G4Step* G4VParticleChange::UpdateStepForPostStep(G4Step* Step)
00240 {
00241 if (isParentWeightProposed ){
00242 Step->GetPostStepPoint()->SetWeight( theParentWeight );
00243 }
00244 return UpdateStepInfo(Step);
00245 }
00246
00247
00248
00249
00250
00251
00252 void G4VParticleChange::DumpInfo() const
00253 {
00254
00255
00256 G4int olprc = G4cout.precision(3);
00257 G4cout << " -----------------------------------------------"
00258 << G4endl;
00259 G4cout << " G4ParticleChange Information " << std::setw(20) << G4endl;
00260 G4cout << " -----------------------------------------------"
00261 << G4endl;
00262
00263 G4cout << " # of 2ndaries : "
00264 << std::setw(20) << theNumberOfSecondaries
00265 << G4endl;
00266
00267 if (theNumberOfSecondaries >0) {
00268 G4cout << " Pointer to 2ndaries : "
00269 << std::setw(20) << GetSecondary(0)
00270 << G4endl;
00271 G4cout << " (Showed only 1st one)"
00272 << G4endl;
00273 }
00274 G4cout << " -----------------------------------------------"
00275 << G4endl;
00276
00277 G4cout << " Energy Deposit (MeV): "
00278 << std::setw(20) << theLocalEnergyDeposit/MeV
00279 << G4endl;
00280
00281 G4cout << " Non-ionizing Energy Deposit (MeV): "
00282 << std::setw(20) << theNonIonizingEnergyDeposit/MeV
00283 << G4endl;
00284
00285 G4cout << " Track Status : "
00286 << std::setw(20);
00287 if( theStatusChange == fAlive ){
00288 G4cout << " Alive";
00289 } else if( theStatusChange == fStopButAlive ){
00290 G4cout << " StopButAlive";
00291 } else if( theStatusChange == fStopAndKill ){
00292 G4cout << " StopAndKill";
00293 } else if( theStatusChange == fKillTrackAndSecondaries ){
00294 G4cout << " KillTrackAndSecondaries";
00295 } else if( theStatusChange == fSuspend ){
00296 G4cout << " Suspend";
00297 } else if( theStatusChange == fPostponeToNextEvent ){
00298 G4cout << " PostponeToNextEvent";
00299 }
00300 G4cout << G4endl;
00301 G4cout << " True Path Length (mm) : "
00302 << std::setw(20) << theTrueStepLength/mm
00303 << G4endl;
00304 G4cout << " Stepping Control : "
00305 << std::setw(20) << theSteppingControlFlag
00306 << G4endl;
00307 if (theFirstStepInVolume) {
00308 G4cout << " First Step In the voulme : " << G4endl;
00309 }
00310 if (theLastStepInVolume) {
00311 G4cout << " Last Step In the voulme : " << G4endl;
00312 }
00313 G4cout.precision(olprc);
00314 }
00315
00316 G4bool G4VParticleChange::CheckIt(const G4Track&
00317 #ifdef G4VERBOSE
00318 aTrack
00319 #endif
00320 )
00321 {
00322
00323 G4bool exitWithError = false;
00324 G4double accuracy;
00325 static G4int nError = 0;
00326 #ifdef G4VERBOSE
00327 const G4int maxError = 30;
00328 #endif
00329
00330
00331 G4bool itsOKforEnergy = true;
00332 accuracy = -1.0*theLocalEnergyDeposit/MeV;
00333 if (accuracy > accuracyForWarning) {
00334 itsOKforEnergy = false;
00335 nError += 1;
00336 exitWithError = (accuracy > accuracyForException);
00337 #ifdef G4VERBOSE
00338 if (nError < maxError) {
00339 G4cout << " G4VParticleChange::CheckIt : ";
00340 G4cout << "the energy deposit is negative !!"
00341 << " Difference: " << accuracy << "[MeV] " <<G4endl;
00342 G4cout << aTrack.GetDefinition()->GetParticleName()
00343 << " E=" << aTrack.GetKineticEnergy()/MeV
00344 << " pos=" << aTrack.GetPosition().x()/m
00345 << ", " << aTrack.GetPosition().y()/m
00346 << ", " << aTrack.GetPosition().z()/m
00347 <<G4endl;
00348 }
00349 #endif
00350 }
00351
00352
00353 G4bool itsOKforStepLength = true;
00354 accuracy = -1.0*theTrueStepLength/mm;
00355 if (accuracy > accuracyForWarning) {
00356 itsOKforStepLength = false;
00357 nError += 1;
00358 exitWithError = (accuracy > accuracyForException);
00359 #ifdef G4VERBOSE
00360 if (nError < maxError) {
00361 G4cout << " G4VParticleChange::CheckIt : ";
00362 G4cout << "the true step length is negative !!"
00363 << " Difference: " << accuracy << "[MeV] " <<G4endl;
00364 G4cout << aTrack.GetDefinition()->GetParticleName()
00365 << " E=" << aTrack.GetKineticEnergy()/MeV
00366 << " pos=" << aTrack.GetPosition().x()/m
00367 << ", " << aTrack.GetPosition().y()/m
00368 << ", " << aTrack.GetPosition().z()/m
00369 <<G4endl;
00370 }
00371 #endif
00372
00373 }
00374 #ifdef G4VERBOSE
00375 if (!itsOKforStepLength || !itsOKforEnergy ){
00376
00377 DumpInfo();
00378 }
00379 #endif
00380
00381
00382 if (exitWithError) {
00383 G4Exception("G4VParticleChange::CheckIt",
00384 "TRACK001", EventMustBeAborted,
00385 "Step length and/or energy deposit was illegal");
00386 }
00387
00388
00389 if ( !itsOKforStepLength ) {
00390 theTrueStepLength = (1.e-12)*mm;
00391 }
00392 if ( !itsOKforEnergy ) {
00393 theLocalEnergyDeposit = 0.0;
00394 }
00395 return (itsOKforStepLength && itsOKforEnergy );
00396 }
00397
00398 G4bool G4VParticleChange::CheckSecondary(G4Track& aTrack)
00399 {
00400 G4bool exitWithError = false;
00401 G4double accuracy;
00402 static G4int nError = 0;
00403 #ifdef G4VERBOSE
00404 const G4int maxError = 30;
00405 #endif
00406
00407
00408 G4bool itsOKforMomentum = true;
00409 if (aTrack.GetKineticEnergy()>0.){
00410 accuracy = std::fabs((aTrack.GetMomentumDirection()).mag2()-1.0);
00411 if (accuracy > accuracyForWarning) {
00412 itsOKforMomentum = false;
00413 nError += 1;
00414 exitWithError = exitWithError || (accuracy > accuracyForException);
00415 #ifdef G4VERBOSE
00416 if (nError < maxError) {
00417 G4cout << " G4VParticleChange::CheckSecondary : ";
00418 G4cout << "the Momentum direction is not unit vector !! "
00419 << " Difference: " << accuracy << G4endl;
00420 G4cout << aTrack.GetDefinition()->GetParticleName()
00421 << " E=" << aTrack.GetKineticEnergy()/MeV
00422 << " pos=" << aTrack.GetPosition().x()/m
00423 << ", " << aTrack.GetPosition().y()/m
00424 << ", " << aTrack.GetPosition().z()/m
00425 <<G4endl;
00426 }
00427 #endif
00428 }
00429 }
00430
00431
00432 G4bool itsOKforEnergy = true;
00433 accuracy = -1.0*(aTrack.GetKineticEnergy())/MeV;
00434 if (accuracy > accuracyForWarning) {
00435 itsOKforEnergy = false;
00436 nError += 1;
00437 exitWithError = exitWithError || (accuracy > accuracyForException);
00438 #ifdef G4VERBOSE
00439 if (nError < maxError) {
00440 G4cout << " G4VParticleChange::CheckSecondary : ";
00441 G4cout << "the kinetic energy is negative !!"
00442 << " Difference: " << accuracy << "[MeV] " <<G4endl;
00443 G4cout << " G4VParticleChange::CheckSecondary : ";
00444 G4cout << "the global time of secondary is earlier than the parent !!"
00445 << " Difference: " << accuracy << "[ns] " <<G4endl;
00446 G4cout << aTrack.GetDefinition()->GetParticleName()
00447 << " E=" << aTrack.GetKineticEnergy()/MeV
00448 << " pos=" << aTrack.GetPosition().x()/m
00449 << ", " << aTrack.GetPosition().y()/m
00450 << ", " << aTrack.GetPosition().z()/m
00451 <<G4endl;
00452 }
00453 #endif
00454 }
00455
00456 G4bool itsOKforTiming = true;
00457
00458 accuracy = (theParentGlobalTime - aTrack.GetGlobalTime())/ns;
00459 if (accuracy > accuracyForWarning){
00460 itsOKforTiming = false;
00461 nError += 1;
00462 exitWithError = (accuracy > accuracyForException);
00463 #ifdef G4VERBOSE
00464 if (nError < maxError) {
00465 G4cout << " G4VParticleChange::CheckSecondary : ";
00466 G4cout << "the global time of secondary goes back comapared to the parent !!"
00467 << " Difference: " << accuracy << "[ns] " <<G4endl;
00468 G4cout << aTrack.GetDefinition()->GetParticleName()
00469 << " E=" << aTrack.GetKineticEnergy()/MeV
00470 << " pos=" << aTrack.GetPosition().x()/m
00471 << ", " << aTrack.GetPosition().y()/m
00472 << ", " << aTrack.GetPosition().z()/m
00473 << " time=" << aTrack.GetGlobalTime()/ns
00474 << " parent time=" << theParentGlobalTime/ns
00475 <<G4endl;
00476 }
00477 #endif
00478 }
00479
00480
00481 if (exitWithError) {
00482 G4Exception("G4VParticleChange::CheckSecondary",
00483 "TRACK001", EventMustBeAborted,
00484 "Secondary with illegal energy/momentum ");
00485 }
00486
00487 G4bool itsOK = itsOKforMomentum && itsOKforEnergy && itsOKforTiming;
00488
00489
00490 if (!itsOKforMomentum) {
00491 G4double vmag = (aTrack.GetMomentumDirection()).mag();
00492 aTrack.SetMomentumDirection((1./vmag)*aTrack.GetMomentumDirection());
00493 }
00494 if (!itsOKforEnergy) {
00495 aTrack.SetKineticEnergy(0.0);
00496 }
00497
00498 return itsOK;
00499 }
00500
00501
00502 G4double G4VParticleChange::GetAccuracyForWarning() const
00503 {
00504 return accuracyForWarning;
00505 }
00506
00507 G4double G4VParticleChange::GetAccuracyForException() const
00508 {
00509 return accuracyForException;
00510 }
00511
00512
00514
00516 void G4VParticleChange::SetParentWeightByProcess(G4bool )
00517 {
00518 }
00519
00520
00521 G4bool G4VParticleChange::IsParentWeightSetByProcess() const
00522 {
00523 return true;
00524 }