diff --git a/G4HepEm/G4HepEm/include/G4HepEmTrackingManager.hh b/G4HepEm/G4HepEm/include/G4HepEmTrackingManager.hh index 1904478..ecec2f9 100644 --- a/G4HepEm/G4HepEm/include/G4HepEmTrackingManager.hh +++ b/G4HepEm/G4HepEm/include/G4HepEmTrackingManager.hh @@ -38,6 +38,22 @@ public: // Control verbosity (0/1) (propagated to the G4HepEmRuManager) void SetVerbose(G4int verbose); + + // Returns the pointer to the Geant4 gamma-nuclear process (if any) + G4VProcess* GetGammaNuclearProcess() { return fGNucProcess; } + // Returns the pointer to the Geant4 electron-nuclear process (if any) + G4VProcess* GetElectronNuclearProcess() { return fENucProcess; } + // Returns the pointer to the Geant4 positron-nuclear process (if any) + G4VProcess* GetPositronNuclearProcess() { return fPNucProcess; } + + // Invokes the G4 electron/positron/gamma-nuclear process (if any and depending + // on the `particleID`={0/1/2} --> electron/positron/gamma-nuclear) updates the + // input step and track to the post interaction state, stacks the secondaries to + // the track vector of the step and returns the energy deposited in the interaction. + // NOTE: the step is assumed to be the one from the track (need non const. values) + double PerformNuclear(G4Track* aG4Track, G4Step* theG4Step, int particleID, bool isApplyCuts); + + // ATLAS XTR RELATED: // Set the names of the ATLAS specific transition radiation process and // radiator region (only for ATLAS and only if different than init.ed below) @@ -62,7 +78,7 @@ private: // Stacks secondaries created by Geant4 physics (if any) and returns with the // energy deposit while stacking due to applying secondary production cuts double StackG4Secondaries(G4VParticleChange* particleChange, - G4Track* aG4PrimaryTrack, + G4Track* aG4PrimaryTrack, G4Step* theStep, const G4VProcess* aG4CreatorProcess, int aG4IMC, bool isApplyCuts); diff --git a/G4HepEm/G4HepEm/src/G4HepEmTrackingManager.cc b/G4HepEm/G4HepEm/src/G4HepEmTrackingManager.cc index 98fda2e..4caf98a 100644 --- a/G4HepEm/G4HepEm/src/G4HepEmTrackingManager.cc +++ b/G4HepEm/G4HepEm/src/G4HepEmTrackingManager.cc @@ -332,12 +332,6 @@ bool G4HepEmTrackingManager::TrackElectron(G4Track *aTrack) { fFastSimProc->StartTracking(aTrack); } - // Invoke the electron/positron-nuclear process start tracking interace (if any) - G4VProcess* theNucProcess = isElectron ? fENucProcess : fPNucProcess; - if (theNucProcess != nullptr) { - theNucProcess->StartTracking(aTrack); - } - // === StartTracking === while(aTrack->GetTrackStatus() == fAlive) @@ -642,29 +636,17 @@ bool G4HepEmTrackingManager::TrackElectron(G4Track *aTrack) { } } else { // Electron/positron-nuclear: --> use Geant4 for the interaction: + // set the process pointer that is used for setting the process that defined this step proc = isElectron ? fElectronNoProcessVector[4] : fElectronNoProcessVector[5]; + // clear the corresponding number of interaction left thePrimaryTrack->SetNumIALeft(-1.0, iDProc); - // check if delta interaction happens - // Invoke the electron/positron-nuclear interaction using the Geant4 process - G4VParticleChange* particleChangeNuc = nullptr; + // check if there is nuclear interaction process and not delta interaction happens + G4VProcess* theNucProcess = isElectron ? fENucProcess : fPNucProcess; if (theNucProcess != nullptr && !G4HepEmElectronManager::CheckDelta(theHepEmData, thePrimaryTrack, theTLData->GetRNGEngine()->flat())) { - // call to set some fields of the process like material, energy etc.. used in its DoIt - G4ForceCondition forceCondition; - theNucProcess->PostStepGetPhysicalInteractionLength(*aTrack, 0.0, &forceCondition); - // - postStepPoint.SetStepStatus(fPostStepDoItProc); - particleChangeNuc = theNucProcess->PostStepDoIt(*aTrack, step); - // update the track and stack according to the result of the interaction - particleChangeNuc->UpdateStepForPostStep(&step); - step.UpdateTrack(); - aTrack->SetTrackStatus(particleChangeNuc->GetTrackStatus()); - // need to add secondaries to the secondary vector of the current track - // NOTE: as we use Geant4, we should care only those changes that are - // not included in the above update step and track, i.e. the energy - // deposited due to applying the cut when stacking the secondaries - thePrimaryTrack->AddEnergyDeposit( StackG4Secondaries(particleChangeNuc, aTrack, proc, g4IMC, isApplyCuts) ); - // done: clear the particle change - particleChangeNuc->Clear(); + // Invoke the electron/positron-nuclear interaction using the Geant4 process + // (step is updated and secondaries are stacked to the vector of the step) + int particleID = isElectron ? 0 : 1; + thePrimaryTrack->AddEnergyDeposit( PerformNuclear(aTrack, &step, particleID, isApplyCuts) ); // update the primary track kinetic energy and direction thePrimaryTrack->SetEKin(aTrack->GetKineticEnergy()); const G4ThreeVector& dir = aTrack->GetMomentumDirection(); @@ -887,10 +869,6 @@ bool G4HepEmTrackingManager::TrackGamma(G4Track *aTrack) { fFastSimProc->StartTracking(aTrack); } - if (fGNucProcess != nullptr) { - fGNucProcess->StartTracking(aTrack); - } - // Reset some Woodcock tracking related flags. G4bool isWDTOn = false; // === StartTracking === @@ -1128,26 +1106,11 @@ bool G4HepEmTrackingManager::TrackGamma(G4Track *aTrack) { // NOTE: it's destructive i.e. stopps and kills the gammma when the // interaction happens. thePrimaryTrack->SetEnergyDeposit(0.0); - // Invoke the gamma-nuclear interaction using the Geant4 process - G4VParticleChange* particleChangeGNuc = nullptr; if (fGNucProcess != nullptr) { - // call to set some fields of the process like material, energy etc... - G4ForceCondition forceCondition; - fGNucProcess->PostStepGetPhysicalInteractionLength(*aTrack, 0.0, &forceCondition); - // - postStepPoint.SetStepStatus(fPostStepDoItProc); - particleChangeGNuc = fGNucProcess->PostStepDoIt(*aTrack, step); - // update the track and stack according to the result of the interaction - particleChangeGNuc->UpdateStepForPostStep(&step); - step.UpdateTrack(); - aTrack->SetTrackStatus(particleChangeGNuc->GetTrackStatus()); - // need to add secondaries to the secondary vector of the current track - // NOTE: as we use Geant4, we should care only those changes that are - // not included in the above update step and track, i.e. the energy - // deposited due to applying the cut when stacking the secondaries - edep = StackG4Secondaries(particleChangeGNuc, aTrack, fGammaNoProcessVector[iDProc], g4IMC, isApplyCuts); - // done: clear the particle change - particleChangeGNuc->Clear(); + // Invoke the gamma-nuclear interaction using the Geant4 process + // (step is updated and secondaries are stacked to the vector of the step) + const int partileID = 2; + edep += PerformNuclear(aTrack, &step, partileID, isApplyCuts); } } // Set process defined setp and add edep to the step @@ -1220,6 +1183,52 @@ void G4HepEmTrackingManager::HandOverOneTrack(G4Track *aTrack) { } +double G4HepEmTrackingManager::PerformNuclear(G4Track* aG4Track, G4Step* theG4Step, int particleID, bool isApplyCuts) { + G4VProcess* theNuclearProcess = nullptr; + // could use the above g4process here but stay consistent + G4VProcess* theCreatorProcess = nullptr; + if (particleID == 2) { // gamma + theNuclearProcess = fGNucProcess; + theCreatorProcess = fGammaNoProcessVector[3]; + } else if (particleID == 0) { // e- + theNuclearProcess = fENucProcess; + theCreatorProcess = fElectronNoProcessVector[4]; + } else if (particleID == 1) { // e+ + theNuclearProcess = fPNucProcess; + theCreatorProcess = fElectronNoProcessVector[5]; + } + if (theNuclearProcess == nullptr) { + return 0.0; + } + double edep = 0.0; + G4VParticleChange* particleChangeGNuc = nullptr; + // calling `StartTracking` that sets the particle and dynamic partile fields of the `G4HadronicProcess` + theNuclearProcess->StartTracking(aG4Track); + // call to set some fields of the process like material, energy etc... + G4ForceCondition forceCondition; + theNuclearProcess->PostStepGetPhysicalInteractionLength(*aG4Track, 0.0, &forceCondition); + // perform the interaction + aG4Track->GetStep()->GetPostStepPoint()->SetStepStatus(fPostStepDoItProc); + particleChangeGNuc = theNuclearProcess->PostStepDoIt(*aG4Track, *theG4Step); + // update the track and stack according to the result of the interaction + particleChangeGNuc->UpdateStepForPostStep(theG4Step); + theG4Step->UpdateTrack(); + aG4Track->SetTrackStatus(particleChangeGNuc->GetTrackStatus()); + // need to add secondaries to the secondary vector of the current track + // NOTE: as we use Geant4, we should care only those changes that are + // not included in the above update step and track, i.e. the energy + // deposited due to applying the cut when stacking the secondaries + // + // the g4material-cuts couple index + const int g4IMC = aG4Track->GetTouchable()->GetVolume()->GetLogicalVolume()->GetMaterialCutsCouple()->GetIndex(); + edep = StackG4Secondaries(particleChangeGNuc, aG4Track, theG4Step, theCreatorProcess, g4IMC, isApplyCuts); + // done: clear the particle change + particleChangeGNuc->Clear(); + // return energy deposited in the interaction (or due to applying the cut) + return edep; +} + + // Helper that can be used to stack secondary e-/e+ and gamma i.e. everything // that HepEm physics can produce double G4HepEmTrackingManager::StackSecondaries(G4HepEmTLData* aTLData, G4Track* aG4PrimaryTrack, const G4VProcess* aG4CreatorProcess, int aG4IMC, bool isApplyCuts) { @@ -1303,7 +1312,7 @@ double G4HepEmTrackingManager::StackSecondaries(G4HepEmTLData* aTLData, G4Track* // Helper that can be used to stack secondary e-/e+ and gamma i.e. everything // that HepEm physics can produce -double G4HepEmTrackingManager::StackG4Secondaries(G4VParticleChange* particleChange, G4Track* aG4PrimaryTrack, const G4VProcess* aG4CreatorProcess, int aG4IMC, bool isApplyCuts) { +double G4HepEmTrackingManager::StackG4Secondaries(G4VParticleChange* particleChange, G4Track* aG4PrimaryTrack, G4Step* theStep, const G4VProcess* aG4CreatorProcess, int aG4IMC, bool isApplyCuts) { const int numSecondaries = particleChange->GetNumberOfSecondaries(); // return early if there are no secondaries created by the physics interaction double edep = 0.0; @@ -1311,9 +1320,8 @@ double G4HepEmTrackingManager::StackG4Secondaries(G4VParticleChange* particleCha return edep; } - G4Step& step = *fStep; - G4TrackVector& secondaries = *step.GetfSecondary(); - G4StepPoint& postStepPoint = *step.GetPostStepPoint(); + G4TrackVector& secondaries = *theStep->GetfSecondary(); + G4StepPoint& postStepPoint = *theStep->GetPostStepPoint(); const G4ThreeVector& theG4PostStepPointPosition = postStepPoint.GetPosition(); const G4double theG4PostStepGlobalTime = postStepPoint.GetGlobalTime();