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 #include "G4ITStepProcessor.hh"
00037 #include "G4LossTableManager.hh"
00038 #include "G4EnergyLossTables.hh"
00039 #include "G4ProductionCuts.hh"
00040 #include "G4ProductionCutsTable.hh"
00041 #include "G4VITProcess.hh"
00042 #include "G4TrackingInformation.hh"
00043 #include "G4IT.hh"
00044 #include "G4ITTrackingManager.hh"
00045 #include "G4ITTransportation.hh"
00046
00047 #include "G4ITNavigator.hh"
00048
00049
00050 void G4ITStepProcessor::DealWithSecondaries(G4int& counter)
00051 {
00052
00053 G4Track* tempSecondaryTrack;
00054
00055 for(G4int DSecLoop=0 ;
00056 DSecLoop<fpParticleChange->GetNumberOfSecondaries() ;
00057 DSecLoop++)
00058 {
00059 tempSecondaryTrack = fpParticleChange->GetSecondary(DSecLoop);
00060
00061 if(tempSecondaryTrack->GetDefinition()->GetApplyCutsFlag())
00062 {
00063 ApplyProductionCut(tempSecondaryTrack);
00064 }
00065
00066
00067 tempSecondaryTrack->SetParentID( fpTrack->GetTrackID() );
00068
00069
00070 tempSecondaryTrack->SetCreatorProcess( fpCurrentProcess );
00071
00072
00073
00074 if(tempSecondaryTrack->GetKineticEnergy() <= DBL_MIN)
00075 {
00076 G4ProcessManager* pm = tempSecondaryTrack->GetDefinition()->GetProcessManager();
00077 if (pm->GetAtRestProcessVector()->entries()>0){
00078 tempSecondaryTrack->SetTrackStatus( fStopButAlive );
00079 fpSecondary->push_back( tempSecondaryTrack );
00080 fN2ndariesAtRestDoIt++;
00081 } else {
00082 delete tempSecondaryTrack;
00083 }
00084 }
00085 else
00086 {
00087 fpSecondary->push_back( tempSecondaryTrack );
00088 counter++;
00089 }
00090 }
00091 }
00092
00093
00094
00095 void G4ITStepProcessor::Stepping(G4Track* track, const double & timeStep)
00096 {
00097 CleanProcessor();
00098 if(track == 0) return ;
00099 fTimeStep = timeStep ;
00100 SetTrack(track);
00101 DoStepping();
00102 }
00103
00104
00105
00106
00107
00108 void G4ITStepProcessor::DoStepping()
00109 {
00110 SetupMembers() ;
00111
00112 if(!fpProcessInfo)
00113 {
00114 G4ExceptionDescription exceptionDescription ;
00115 exceptionDescription << "No process info found for particle :"
00116 << fpTrack->GetDefinition()->GetParticleName();
00117 G4Exception("G4ITStepProcessor::DoStepping","ITStepProcessor0012",
00118 FatalErrorInArgument,exceptionDescription);
00119 return ;
00120 }
00121 else if(fpTrack->GetTrackStatus() == fStopAndKill )
00122 {
00123 fpState->fStepStatus = fUndefined;
00124 return ;
00125 }
00126
00127 if(fpProcessInfo->MAXofPostStepLoops == 0
00128 && fpProcessInfo->MAXofAlongStepLoops == 0
00129 && fpProcessInfo->MAXofAtRestLoops == 0)
00130 {
00131 fpTrack -> SetTrackStatus(fStopAndKill) ;
00132 fpState->fStepStatus = fUndefined;
00133 return ;
00134 }
00135
00136
00137
00138 else
00139 {
00140 fpNavigator->SetNavigatorState(fpITrack->GetTrackingInfo()->GetNavigatorState());
00141 fpNavigator->ResetHierarchyAndLocate( fpTrack->GetPosition(),
00142 fpTrack->GetMomentumDirection(),
00143 *((G4TouchableHistory*)fpTrack->GetTouchableHandle()()) );
00144 fpNavigator->SetNavigatorState(fpITrack->GetTrackingInfo()->GetNavigatorState());
00145
00146
00147
00148 if( fpTrack->GetTrackStatus() == fStopButAlive )
00149 {
00150 if( fpProcessInfo->MAXofAtRestLoops>0 &&
00151 fpProcessInfo->fpAtRestDoItVector != 0)
00152 {
00153
00154
00155
00156 InvokeAtRestDoItProcs();
00157 fpState->fStepStatus = fAtRestDoItProc;
00158 fpStep->GetPostStepPoint()->SetStepStatus( fpState->fStepStatus );
00159
00160 }
00161
00162 fpTrack->SetTrackStatus( fStopAndKill );
00163 }
00164 else
00165 {
00166 if(fpITrack == 0)
00167 {
00168 G4ExceptionDescription exceptionDescription ;
00169 exceptionDescription
00170 << " !!! TrackID : "<< fpTrack->GetTrackID() << G4endl
00171 << " !!! Track status : "<< fpTrack->GetTrackStatus() << G4endl
00172 << " !!! Particle Name : "<< fpTrack -> GetDefinition() -> GetParticleName() << G4endl
00173 << "No G4ITStepProcessor::fpITrack found" << G4endl;
00174
00175 G4Exception("G4ITStepProcessor::DoStepping","ITStepProcessor0013",
00176 FatalErrorInArgument,exceptionDescription);
00177 return ;
00178 }
00179
00180 if(fpITrack->GetTrackingInfo()->IsLeadingStep() == false)
00181 {
00182
00183
00184
00185
00186 fpState->fStepStatus = fPostStepDoItProc;
00187 fpStep->GetPostStepPoint()
00188 ->SetProcessDefinedStep(fpTransportation);
00189 FindTransportationStep();
00190 }
00191
00192
00193
00194 fpTrack->SetStepLength( fpState->fPhysicalStep );
00195 fpStep->SetStepLength( fpState->fPhysicalStep );
00196
00197 G4double GeomStepLength = fpState->fPhysicalStep;
00198
00199
00200 fpStep->GetPostStepPoint()->SetStepStatus( fpState->fStepStatus );
00201
00202
00203 InvokeAlongStepDoItProcs();
00204
00205
00206
00207
00208
00209 fpState->endpointSafOrigin= fpPostStepPoint->GetPosition();
00210
00211 fpState->endpointSafety= std::max( fpState->proposedSafety - GeomStepLength, kCarTolerance);
00212
00213 fpStep->GetPostStepPoint()->SetSafety( fpState->endpointSafety );
00214
00215 if(GetIT(fpTrack)->GetTrackingInfo()->IsLeadingStep())
00216 {
00217
00218 InvokePostStepDoItProcs();
00219 }
00220 else
00221 {
00222
00223 InvokeTransportationProc();
00224 }
00225 }
00226
00227 fpITrack->GetTrackingInfo()->SetNavigatorState(fpNavigator->GetNavigatorState());
00228 fpNavigator->SetNavigatorState(0);
00229 }
00230
00231
00232
00233
00234
00235 fpTrack->AddTrackLength(fpStep->GetStepLength());
00236 fpTrack->IncrementCurrentStepNumber();
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263 fpTrackingManager->AppendStep(fpTrack,fpStep);
00264
00265
00266
00267 }
00268
00269
00270
00271
00272
00273
00274
00275 void G4ITStepProcessor::InvokeAtRestDoItProcs()
00276 {
00277 fpStep->SetStepLength( 0. );
00278 fpTrack->SetStepLength( 0. );
00279
00280 G4SelectedAtRestDoItVector& selectedAtRestDoItVector = fpState->fSelectedAtRestDoItVector;
00281
00282
00283 for(size_t np=0; np < fpProcessInfo->MAXofAtRestLoops; np++)
00284 {
00285
00286
00287
00288
00289 if( selectedAtRestDoItVector[fpProcessInfo->MAXofAtRestLoops-np-1] != InActivated)
00290 {
00291 fpCurrentProcess = (G4VITProcess*) (*fpProcessInfo->fpAtRestDoItVector)[np];
00292
00293 fpCurrentProcess->SetProcessState(
00294 fpTrackingInfo->GetProcessState(fpCurrentProcess->GetProcessID()));
00295 fpParticleChange
00296 = fpCurrentProcess->AtRestDoIt( *fpTrack, *fpStep);
00297 fpCurrentProcess->SetProcessState(0);
00298
00299
00300 fpStep->GetPostStepPoint()
00301 ->SetProcessDefinedStep(fpCurrentProcess);
00302
00303
00304 fpParticleChange->UpdateStepForAtRest(fpStep);
00305
00306
00307 DealWithSecondaries(fN2ndariesAtRestDoIt) ;
00308
00309
00310 fpParticleChange->Clear();
00311
00312 }
00313 }
00314 fpStep->UpdateTrack();
00315
00316 fpTrack->SetTrackStatus( fStopAndKill );
00317 }
00318
00319
00320
00321
00322
00323
00324
00325 void G4ITStepProcessor::InvokeAlongStepDoItProcs()
00326 {
00327
00328
00329
00330 if(fpState->fStepStatus == fExclusivelyForcedProc)
00331 {
00332 return;
00333 }
00334
00335
00336 for( size_t ci=0 ; ci<fpProcessInfo->MAXofAlongStepLoops ; ci++ )
00337 {
00338 fpCurrentProcess = (G4VITProcess*) (*fpProcessInfo->fpAlongStepDoItVector)[ci];
00339 if (fpCurrentProcess== 0) continue;
00340
00341
00342 fpCurrentProcess->SetProcessState(fpTrackingInfo->GetProcessState(fpCurrentProcess->GetProcessID()));
00343 fpParticleChange
00344 = fpCurrentProcess->AlongStepDoIt( *fpTrack, *fpStep );
00345 fpCurrentProcess->SetProcessState(0);
00346
00347
00348 fpParticleChange->UpdateStepForAlongStep(fpStep);
00349
00350
00351 DealWithSecondaries(fN2ndariesAlongStepDoIt) ;
00352
00353
00354
00355 fpTrack->SetTrackStatus( fpParticleChange->GetTrackStatus() );
00356
00357
00358 fpParticleChange->Clear();
00359 }
00360
00361 fpStep->UpdateTrack();
00362
00363 G4TrackStatus fNewStatus = fpTrack->GetTrackStatus();
00364
00365 if ( fNewStatus == fAlive && fpTrack->GetKineticEnergy() <= DBL_MIN )
00366 {
00367
00368 if(fpProcessInfo->MAXofAtRestLoops>0) fNewStatus = fStopButAlive;
00369 else fNewStatus = fStopAndKill;
00370 fpTrack->SetTrackStatus( fNewStatus );
00371 }
00372
00373 }
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 void G4ITStepProcessor::InvokePostStepDoItProcs()
00384 {
00385
00386 G4SelectedPostStepDoItVector& selectedPostStepDoItVector = fpState->fSelectedPostStepDoItVector;
00387 G4StepStatus& stepStatus = fpState->fStepStatus;
00388
00389
00390 for(size_t np=0; np < fpProcessInfo->MAXofPostStepLoops; np++)
00391 {
00392
00393
00394
00395
00396 G4int Cond = selectedPostStepDoItVector[fpProcessInfo->MAXofPostStepLoops-np-1];
00397 if(Cond != InActivated)
00398 {
00399 if( ((Cond == NotForced) && (stepStatus == fPostStepDoItProc)) ||
00400 ((Cond == Forced) && (stepStatus != fExclusivelyForcedProc)) ||
00401
00402 ((Cond == ExclusivelyForced) && (stepStatus == fExclusivelyForcedProc)) ||
00403 ((Cond == StronglyForced) )
00404 )
00405 {
00406
00407 InvokePSDIP(np);
00408 }
00409 }
00410
00411
00412
00413 if(fpTrack->GetTrackStatus() == fStopAndKill)
00414 {
00415 for(size_t np1=np+1; np1 < fpProcessInfo->MAXofPostStepLoops; np1++)
00416 {
00417 G4int Cond2 = selectedPostStepDoItVector[fpProcessInfo->MAXofPostStepLoops-np1-1];
00418 if (Cond2 == StronglyForced)
00419 {
00420 InvokePSDIP(np1);
00421 }
00422 }
00423 break;
00424 }
00425 }
00426 }
00427
00428
00429
00430 void G4ITStepProcessor::InvokePSDIP(size_t np)
00431 {
00432 fpCurrentProcess = (G4VITProcess*) (*fpProcessInfo->fpPostStepDoItVector)[np];
00433
00434 fpCurrentProcess->SetProcessState(fpTrackingInfo->GetProcessState(fpCurrentProcess->GetProcessID()));
00435 fpParticleChange
00436 = fpCurrentProcess->PostStepDoIt( *fpTrack, *fpStep);
00437 fpCurrentProcess->SetProcessState(0);
00438
00439
00440 fpParticleChange->UpdateStepForPostStep(fpStep);
00441
00442
00443 fpStep->UpdateTrack();
00444
00445
00446 fpStep->GetPostStepPoint()->SetSafety( CalculateSafety() );
00447
00448
00449 DealWithSecondaries(fN2ndariesPostStepDoIt) ;
00450
00451
00452 fpTrack->SetTrackStatus( fpParticleChange->GetTrackStatus() );
00453
00454
00455 fpParticleChange->Clear();
00456 }
00457
00458
00459
00460
00461
00462
00463
00464 void G4ITStepProcessor::FindTransportationStep()
00465 {
00466 double physicalStep(0.) ;
00467
00468 fpTransportation = fpProcessInfo->fpTransportation;
00469
00470
00471 if(!fpTrack)
00472 {
00473 G4ExceptionDescription exceptionDescription ;
00474 exceptionDescription
00475 << "No G4ITStepProcessor::fpTrack found";
00476 G4Exception("G4ITStepProcessor::FindTransportationStep","ITStepProcessor0013",
00477 FatalErrorInArgument,exceptionDescription);
00478 return;
00479
00480 }
00481 if(!fpITrack)
00482 {
00483 G4ExceptionDescription exceptionDescription ;
00484 exceptionDescription
00485 << "No G4ITStepProcessor::fITrack" ;
00486 G4Exception("G4ITStepProcessor::FindTransportationStep","ITStepProcessor0014",
00487 FatalErrorInArgument,exceptionDescription);
00488 return;
00489 }
00490 if(!(fpITrack->GetTrack()))
00491 {
00492 G4ExceptionDescription exceptionDescription ;
00493 exceptionDescription
00494 << "No G4ITStepProcessor::fITrack->GetTrack()" ;
00495 G4Exception("G4ITStepProcessor::FindTransportationStep","ITStepProcessor0015",
00496 FatalErrorInArgument,exceptionDescription);
00497 return;
00498 }
00499
00500 if(fpTransportation)
00501 {
00502 fpTransportation->SetProcessState(fpTrackingInfo->GetProcessState(fpTransportation->GetProcessID()));
00503 fpTransportation->ComputeStep(*fpTrack, *fpStep, fTimeStep, physicalStep) ;
00504 fpTransportation->SetProcessState(0);
00505 }
00506
00507 if(physicalStep >= DBL_MAX)
00508 {
00509 fpTrack -> SetTrackStatus(fStopAndKill) ;
00510 return ;
00511 }
00512
00513 fpState->fPhysicalStep = physicalStep ;
00514 }
00515
00516
00517
00518 void G4ITStepProcessor::InvokeTransportationProc()
00519 {
00520 size_t _MAXofPostStepLoops = fpProcessInfo->MAXofPostStepLoops;
00521 G4SelectedPostStepDoItVector& selectedPostStepDoItVector = fpState->fSelectedPostStepDoItVector;
00522 G4StepStatus& stepStatus = fpState->fStepStatus;
00523
00524
00525 for(size_t np=0; np < _MAXofPostStepLoops; np++)
00526 {
00527
00528
00529
00530
00531 G4int Cond = selectedPostStepDoItVector[_MAXofPostStepLoops-np-1];
00532 if(Cond != InActivated)
00533 {
00534 if(
00535 ((Cond == Forced) && (stepStatus != fExclusivelyForcedProc)) ||
00536
00537 ((Cond == ExclusivelyForced) && (stepStatus == fExclusivelyForcedProc)) ||
00538 ((Cond == StronglyForced) )
00539 )
00540 {
00541
00542 InvokePSDIP(np);
00543 }
00544 }
00545
00546
00547
00548 if(fpTrack->GetTrackStatus() == fStopAndKill)
00549 {
00550 for(size_t np1=np+1; np1 < _MAXofPostStepLoops; np1++)
00551 {
00552 G4int Cond2 = selectedPostStepDoItVector[_MAXofPostStepLoops-np1-1];
00553 if (Cond2 == StronglyForced)
00554 {
00555 InvokePSDIP(np1);
00556 }
00557 }
00558 break;
00559 }
00560 }
00561 }
00562
00563
00564
00565
00566
00567
00568
00569
00570 void G4ITStepProcessor::ApplyProductionCut(G4Track* aSecondary)
00571 {
00572 G4bool tBelowCutEnergyAndSafety = false;
00573 G4int tPtclIdx
00574 = G4ProductionCuts::GetIndex(aSecondary->GetDefinition());
00575 if (tPtclIdx<0)
00576 {
00577 return;
00578 }
00579 G4ProductionCutsTable* tCutsTbl
00580 = G4ProductionCutsTable::GetProductionCutsTable();
00581 G4int tCoupleIdx
00582 = tCutsTbl->GetCoupleIndex(fpPreStepPoint->GetMaterialCutsCouple());
00583 G4double tProdThreshold
00584 = (*(tCutsTbl->GetEnergyCutsVector(tPtclIdx)))[tCoupleIdx];
00585 if( aSecondary->GetKineticEnergy()<tProdThreshold )
00586 {
00587 tBelowCutEnergyAndSafety = true;
00588 if(std::abs(aSecondary->GetDynamicParticle()->GetCharge()) > DBL_MIN)
00589 {
00590 G4double currentRange
00591 = G4LossTableManager::Instance()->GetRange(aSecondary->GetDefinition(),
00592 aSecondary->GetKineticEnergy(),
00593 fpPreStepPoint->GetMaterialCutsCouple());
00594 tBelowCutEnergyAndSafety = (currentRange < CalculateSafety() );
00595 }
00596 }
00597
00598 if( tBelowCutEnergyAndSafety )
00599 {
00600 if( !(aSecondary->IsGoodForTracking()) )
00601 {
00602
00603 fpStep->AddTotalEnergyDeposit(
00604 aSecondary->GetKineticEnergy() );
00605 aSecondary->SetKineticEnergy(0.0);
00606 }
00607 }
00608 }