Skip to content
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.

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

Updated Geant4 user actions #40252

Merged
merged 2 commits into from Dec 7, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
21 changes: 12 additions & 9 deletions SimG4Core/Application/interface/StackingAction.h
Expand Up @@ -15,6 +15,7 @@
class NewTrackAction;
class TrackingAction;
class CMSSteppingVerbose;
class G4VProcess;

class StackingAction : public G4UserStackingAction {
public:
Expand Down Expand Up @@ -46,7 +47,8 @@ class StackingAction : public G4UserStackingAction {
bool savePDandCinTracker, savePDandCinCalo;
bool savePDandCinMuon, saveFirstSecondary;
bool savePDandCinAll;
bool killInCalo, killInCaloEfH;
bool killInCalo{false};
bool killInCaloEfH{false};
bool killHeavy, trackNeutrino, killDeltaRay;
bool killExtra;
bool killGamma;
Expand All @@ -71,16 +73,17 @@ class StackingAction : public G4UserStackingAction {
G4VSolid* worldSolid;
const TrackingAction* trackAction;
const CMSSteppingVerbose* steppingVerbose;
const G4VProcess* m_Compton{nullptr};
NewTrackAction* newTA;
TrackInformationExtractor extractor;

// Russian roulette regions
const G4Region* regionEcal;
const G4Region* regionHcal;
const G4Region* regionMuonIron;
const G4Region* regionPreShower;
const G4Region* regionCastor;
const G4Region* regionWorld;
const G4Region* regionEcal{nullptr};
const G4Region* regionHcal{nullptr};
const G4Region* regionMuonIron{nullptr};
const G4Region* regionPreShower{nullptr};
const G4Region* regionCastor{nullptr};
const G4Region* regionWorld{nullptr};

// Russian roulette energy limits
double gRusRoEnerLim;
Expand All @@ -100,8 +103,8 @@ class StackingAction : public G4UserStackingAction {
double gRusRoWorld;
double nRusRoWorld;
// flags
bool gRRactive;
bool nRRactive;
bool gRRactive{false};
bool nRRactive{false};
};

#endif
46 changes: 33 additions & 13 deletions SimG4Core/Application/src/StackingAction.cc
Expand Up @@ -15,6 +15,7 @@
#include "G4VSolid.hh"
#include "G4TransportationManager.hh"
#include "G4GammaGeneralProcess.hh"
#include "G4LossTableManager.hh"

StackingAction::StackingAction(const TrackingAction* trka, const edm::ParameterSet& p, const CMSSteppingVerbose* sv)
: trackAction(trka), steppingVerbose(sv) {
Expand All @@ -38,16 +39,6 @@ StackingAction::StackingAction(const TrackingAction* trka, const edm::ParameterS
savePDandCinCalo = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInCalo", false);
savePDandCinMuon = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInMuon", false);
saveFirstSecondary = p.getUntrackedParameter<bool>("SaveFirstLevelSecondary", false);
killInCalo = false;
killInCaloEfH = false;

// Russian Roulette
regionEcal = nullptr;
regionHcal = nullptr;
regionMuonIron = nullptr;
regionPreShower = nullptr;
regionCastor = nullptr;
regionWorld = nullptr;

gRusRoEnerLim = p.getParameter<double>("RusRoGammaEnergyLimit") * CLHEP::MeV;
nRusRoEnerLim = p.getParameter<double>("RusRoNeutronEnergyLimit") * CLHEP::MeV;
Expand All @@ -66,8 +57,6 @@ StackingAction::StackingAction(const TrackingAction* trka, const edm::ParameterS
nRusRoCastor = p.getParameter<double>("RusRoCastorNeutron");
nRusRoWorld = p.getParameter<double>("RusRoWorldNeutron");

gRRactive = false;
nRRactive = false;
if (gRusRoEnerLim > 0.0 && (gRusRoEcal < 1.0 || gRusRoHcal < 1.0 || gRusRoMuonIron < 1.0 || gRusRoPreShower < 1.0 ||
gRusRoCastor < 1.0 || gRusRoWorld < 1.0)) {
gRRactive = true;
Expand Down Expand Up @@ -177,6 +166,17 @@ G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack(const G4Track* aTrac
auto track = const_cast<G4Track*>(aTrack);
const G4VProcess* creatorProc = aTrack->GetCreatorProcess();

if (creatorProc == nullptr && aTrack->GetParentID() != 0) {
edm::LogWarning("StackingAction::ClassifyNewTrack")
<< " TrackID=" << aTrack->GetTrackID() << " ParentID=" << aTrack->GetParentID() << " "
<< aTrack->GetDefinition()->GetParticleName() << " Ekin(MeV)=" << aTrack->GetKineticEnergy();
}
if (aTrack->GetKineticEnergy() < 0.0) {
edm::LogWarning("StackingAction::ClassifyNewTrack")
<< " TrackID=" << aTrack->GetTrackID() << " ParentID=" << aTrack->GetParentID() << " "
<< aTrack->GetDefinition()->GetParticleName() << " Ekin(MeV)=" << aTrack->GetKineticEnergy() << " creator "
<< creatorProc->GetProcessName();
}
// primary
if (creatorProc == nullptr || aTrack->GetParentID() == 0) {
if (!trackNeutrino && (abspdg == 12 || abspdg == 14 || abspdg == 16 || abspdg == 18)) {
Expand Down Expand Up @@ -214,14 +214,34 @@ G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack(const G4Track* aTrac
// potentially good for tracking
const double ke = aTrack->GetKineticEnergy();
G4int subType = (nullptr != creatorProc) ? creatorProc->GetProcessSubType() : 0;
// VI: this part of code is needed for Geant4 10.7 only
if (subType == 16) {
auto ptr = dynamic_cast<const G4GammaGeneralProcess*>(creatorProc);
if (nullptr != ptr) {
creatorProc = ptr->GetSelectedProcess();
subType = (nullptr != creatorProc) ? creatorProc->GetProcessSubType() : 0;
if (nullptr == creatorProc) {
if (nullptr == m_Compton) {
auto vp = G4LossTableManager::Instance()->GetEmProcessVector();
for (auto& p : vp) {
if (fComptonScattering == p->GetProcessSubType()) {
m_Compton = p;
break;
}
}
}
creatorProc = m_Compton;
}
subType = creatorProc->GetProcessSubType();
track->SetCreatorProcess(creatorProc);
}
if (creatorProc == nullptr) {
edm::LogWarning("StackingAction::ClassifyNewTrack")
<< " SubType=16 and no creatorProc; TrackID=" << aTrack->GetTrackID()
<< " ParentID=" << aTrack->GetParentID() << " " << aTrack->GetDefinition()->GetParticleName()
<< " Ekin(MeV)=" << ke << " SubType=" << subType;
}
}
// VI - end
LogDebug("SimG4CoreApplication") << "##StackingAction:Classify Track " << aTrack->GetTrackID() << " Parent "
<< aTrack->GetParentID() << " " << aTrack->GetDefinition()->GetParticleName()
<< " Ekin(MeV)=" << ke / CLHEP::MeV << " subType=" << subType << " ";
Expand Down
6 changes: 3 additions & 3 deletions SimG4Core/Application/src/SteppingAction.cc
Expand Up @@ -99,11 +99,11 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) {
TrackStatus tstat = (theTrack->GetTrackStatus() == fAlive) ? sAlive : sKilledByProcess;

if (theTrack->GetKineticEnergy() < 0.0) {
if (nWarnings < 5) {
if (nWarnings < 2) {
++nWarnings;
edm::LogWarning("SimG4CoreApplication")
<< "Track #" << theTrack->GetTrackID() << " " << theTrack->GetDefinition()->GetParticleName()
<< " Ekin(MeV)= " << theTrack->GetKineticEnergy() / MeV;
<< "SteppingAction::UserSteppingAction: Track #" << theTrack->GetTrackID() << " "
<< theTrack->GetDefinition()->GetParticleName() << " Ekin(MeV)= " << theTrack->GetKineticEnergy() / MeV;
}
theTrack->SetKineticEnergy(0.0);
}
Expand Down