Skip to content

Factor out nuclear interactions: available now for external users. #124

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mnovak42 merged 4 commits into from
Apr 3, 2025
Merged

Factor out nuclear interactions: available now for external users. #124

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d272006
Making the gamma-nuclear interaction available for external usage as …
mnovak42 Apr 3, 2025
0bac0b0
Make the pointers to the g4 nuclear interaction processes available (…
mnovak42 Apr 3, 2025
1f0b5ed
Make the method self contained, i.e. everything that needs to be don…
mnovak42 Apr 3, 2025
477973f
Include also lepton-nuclear interactions into a single method.
mnovak42 Apr 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 17 additions & 1 deletion G4HepEm/G4HepEm/include/G4HepEmTrackingManager.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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);

Expand Down
114 changes: 61 additions & 53 deletions G4HepEm/G4HepEm/src/G4HepEmTrackingManager.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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 ===
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -1303,17 +1312,16 @@ 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;
if (numSecondaries == 0) {
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();
Expand Down