Skip to content

Commit ec11383

Browse files
cmsbuildcmsbuild
authored
Merge pull request cms-sw#33923 from civanch/bug_fix_stepping
Bug fix in Geant4 stepping and stacking actions
2 parents ad03f8b + 4f5e150 commit ec11383

File tree

4 files changed

+54
-44
lines changed

4 files changed

+54
-44
lines changed

‎SimG4Core/Application/interface/StackingAction.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ class StackingAction : public G4UserStackingAction {
3636

3737
bool rrApplicable(const G4Track*, const G4Track&) const;
3838

39-
bool isItOutOfTimeWindow(const G4Region*, const G4Track*) const;
39+
bool isItOutOfTimeWindow(const G4Region*, const double&) const;
4040

4141
bool isThisRegion(const G4Region*, std::vector<const G4Region*>&) const;
4242

‎SimG4Core/Application/interface/SteppingAction.h‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -42,10 +42,10 @@ class SteppingAction : public G4UserSteppingAction {
4242
bool initPointer();
4343

4444
inline bool isInsideDeadRegion(const G4Region* reg) const;
45-
inline bool isOutOfTimeWindow(const G4Track* theTrack, const G4Region* reg) const;
45+
inline bool isOutOfTimeWindow(const G4Region* reg, const double& time) const;
4646
inline bool isThisVolume(const G4VTouchable* touch, const G4VPhysicalVolume* pv) const;
4747

48-
bool isLowEnergy(const G4Step* aStep) const;
48+
bool isLowEnergy(const G4LogicalVolume*, const G4Track*) const;
4949
void PrintKilledTrack(const G4Track*, const TrackStatus&) const;
5050

5151
EventAction* eventAction_;
@@ -88,15 +88,15 @@ inline bool SteppingAction::isInsideDeadRegion(const G4Region* reg) const {
8888
return res;
8989
}
9090

91-
inline bool SteppingAction::isOutOfTimeWindow(const G4Track* theTrack, const G4Region* reg) const {
91+
inline bool SteppingAction::isOutOfTimeWindow(const G4Region* reg, const double& time) const {
9292
double tofM = maxTrackTime;
9393
for (unsigned int i = 0; i < numberTimes; ++i) {
9494
if (reg == maxTimeRegions[i]) {
9595
tofM = maxTrackTimes[i];
9696
break;
9797
}
9898
}
99-
return (theTrack->GetGlobalTime() > tofM);
99+
return (time > tofM);
100100
}
101101

102102
inline bool SteppingAction::isThisVolume(const G4VTouchable* touch, const G4VPhysicalVolume* pv) const {

‎SimG4Core/Application/src/StackingAction.cc‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -186,6 +186,7 @@ G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack(const G4Track* aTrac
186186
} else {
187187
// secondary
188188
const G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
189+
const double time = aTrack->GetGlobalTime();
189190

190191
// definetly killed tracks
191192
if (aTrack->GetTrackStatus() == fStopAndKill) {
@@ -195,14 +196,14 @@ G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack(const G4Track* aTrac
195196

196197
} else if (std::abs(aTrack->GetPosition().z()) >= maxZCentralCMS) {
197198
// very forward secondary
198-
if (aTrack->GetGlobalTime() > maxTrackTimeForward) {
199+
if (time > maxTrackTimeForward) {
199200
classification = fKill;
200201
} else {
201202
const G4Track* mother = trackAction->geant4Track();
202203
newTA->secondary(aTrack, *mother, 0);
203204
}
204205

205-
} else if (isItOutOfTimeWindow(reg, aTrack)) {
206+
} else if (isItOutOfTimeWindow(reg, time)) {
206207
// time window check
207208
classification = fKill;
208209

@@ -217,7 +218,7 @@ G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack(const G4Track* aTrac
217218
if (classification != fKill && ke <= limitEnergyForVacuum && isThisRegion(reg, lowdensRegions)) {
218219
classification = fKill;
219220

220-
} else {
221+
} else if (classification != fKill) {
221222
// very low-energy gamma
222223
if (pdg == 22 && killGamma && ke < kmaxGamma) {
223224
classification = fKill;
@@ -434,8 +435,8 @@ bool StackingAction::rrApplicable(const G4Track* aTrack, const G4Track& mother)
434435
const TrackInformation& motherInfo(extractor(mother));
435436

436437
// Check whether mother is gamma, e+, e-
437-
const int genID = std::abs(motherInfo.genParticlePID());
438-
return (22 != genID && 11 != genID);
438+
const int genID = motherInfo.genParticlePID();
439+
return (22 != genID && 11 != std::abs(genID));
439440
}
440441

441442
int StackingAction::isItFromPrimary(const G4Track& mother, int flagIn) const {
@@ -449,15 +450,15 @@ int StackingAction::isItFromPrimary(const G4Track& mother, int flagIn) const {
449450
return flag;
450451
}
451452

452-
bool StackingAction::isItOutOfTimeWindow(const G4Region* reg, const G4Track* aTrack) const {
453+
bool StackingAction::isItOutOfTimeWindow(const G4Region* reg, const double& t) const {
453454
double tofM = maxTrackTime;
454455
for (unsigned int i = 0; i < numberTimes; ++i) {
455456
if (reg == maxTimeRegions[i]) {
456457
tofM = maxTrackTimes[i];
457458
break;
458459
}
459460
}
460-
return (aTrack->GetGlobalTime() > tofM);
461+
return (t > tofM);
461462
}
462463

463464
void StackingAction::printRegions(const std::vector<const G4Region*>& reg, const std::string& word) const {

‎SimG4Core/Application/src/SteppingAction.cc‎

Lines changed: 41 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -103,7 +103,7 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) {
103103
const G4StepPoint* postStep = aStep->GetPostStepPoint();
104104

105105
// NaN energy deposit
106-
if (sAlive == tstat && edm::isNotFinite(aStep->GetTotalEnergyDeposit())) {
106+
if (edm::isNotFinite(aStep->GetTotalEnergyDeposit())) {
107107
tstat = sEnergyDepNaN;
108108
if (nWarnings < 5) {
109109
++nWarnings;
@@ -115,48 +115,60 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) {
115115
}
116116
}
117117

118+
// the track is killed by the process
119+
if (tstat == sKilledByProcess) {
120+
if (nullptr != steppingVerbose) {
121+
steppingVerbose->NextStep(aStep, fpSteppingManager, false);
122+
}
123+
return;
124+
}
125+
118126
if (sAlive == tstat && theTrack->GetCurrentStepNumber() > maxNumberOfSteps) {
119127
tstat = sNumberOfSteps;
120128
if (nWarnings < 5) {
121129
++nWarnings;
122130
edm::LogWarning("SimG4CoreApplication")
123131
<< "Track #" << theTrack->GetTrackID() << " " << theTrack->GetDefinition()->GetParticleName()
124132
<< " E(MeV)= " << preStep->GetKineticEnergy() / MeV << " Nstep= " << theTrack->GetCurrentStepNumber()
125-
<< " is killed due to limit on number of steps;PV: " << preStep->GetPhysicalVolume()->GetName() << " at "
133+
<< " is killed due to limit on number of steps;/n PV: " << preStep->GetPhysicalVolume()->GetName() << " at "
126134
<< theTrack->GetPosition() << " StepLen(mm)= " << aStep->GetStepLength();
127135
}
128136
}
137+
const double time = theTrack->GetGlobalTime();
129138

130139
// check Z-coordinate
131140
if (sAlive == tstat && std::abs(theTrack->GetPosition().z()) >= maxZCentralCMS) {
132-
tstat = (theTrack->GetGlobalTime() > maxTrackTimeForward) ? sOutOfTime : sVeryForward;
141+
tstat = (time > maxTrackTimeForward) ? sOutOfTime : sVeryForward;
133142
}
134143

135144
// check G4Region
136-
if (sAlive == tstat && postStep->GetPhysicalVolume() != nullptr) {
137-
const G4Region* theRegion = preStep->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
145+
if (sAlive == tstat) {
146+
// next logical volume and next region
147+
const G4LogicalVolume* lv = postStep->GetPhysicalVolume()->GetLogicalVolume();
148+
const G4Region* theRegion = lv->GetRegion();
138149

139150
// kill in dead regions
140-
if (isInsideDeadRegion(theRegion)) {
151+
if (isInsideDeadRegion(theRegion))
141152
tstat = sDeadRegion;
142-
}
143153

144154
// kill out of time
145-
if (sAlive == tstat && isOutOfTimeWindow(theTrack, theRegion)) {
146-
tstat = sOutOfTime;
155+
if (sAlive == tstat) {
156+
if (isOutOfTimeWindow(theRegion, time))
157+
tstat = sOutOfTime;
147158
}
148159

149160
// kill low-energy in volumes on demand
150-
if (sAlive == tstat && numberEkins > 0 && isLowEnergy(aStep)) {
151-
tstat = sLowEnergy;
161+
if (sAlive == tstat && numberEkins > 0) {
162+
if (isLowEnergy(lv, theTrack))
163+
tstat = sLowEnergy;
152164
}
153165

154166
// kill low-energy in vacuum
155-
const G4double kinEnergy = theTrack->GetKineticEnergy();
156-
if (sAlive == tstat && killBeamPipe && kinEnergy < theCriticalEnergyForVacuum &&
157-
theTrack->GetDefinition()->GetPDGCharge() != 0.0 && kinEnergy > 0.0 &&
158-
theTrack->GetNextVolume()->GetLogicalVolume()->GetMaterial()->GetDensity() <= theCriticalDensity) {
159-
tstat = sLowEnergyInVacuum;
167+
if (sAlive == tstat && killBeamPipe) {
168+
if (theTrack->GetKineticEnergy() < theCriticalEnergyForVacuum &&
169+
theTrack->GetDefinition()->GetPDGCharge() != 0.0 && lv->GetMaterial()->GetDensity() <= theCriticalDensity) {
170+
tstat = sLowEnergyInVacuum;
171+
}
160172
}
161173
}
162174
// check transition tracker/calo
@@ -185,11 +197,9 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) {
185197
}
186198
}
187199

188-
bool SteppingAction::isLowEnergy(const G4Step* aStep) const {
189-
const G4StepPoint* sp = aStep->GetPostStepPoint();
190-
const G4LogicalVolume* lv = sp->GetPhysicalVolume()->GetLogicalVolume();
191-
const double ekin = sp->GetKineticEnergy();
192-
int pCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
200+
bool SteppingAction::isLowEnergy(const G4LogicalVolume* lv, const G4Track* theTrack) const {
201+
const double ekin = theTrack->GetKineticEnergy();
202+
int pCode = theTrack->GetDefinition()->GetPDGEncoding();
193203

194204
for (auto& vol : ekinVolumes) {
195205
if (lv == vol) {
@@ -281,9 +291,6 @@ void SteppingAction::PrintKilledTrack(const G4Track* aTrack, const TrackStatus&
281291
std::string rname = "";
282292
std::string typ = " ";
283293
switch (tst) {
284-
case sKilledByProcess:
285-
typ = " G4Process ";
286-
break;
287294
case sDeadRegion:
288295
typ = " in dead region ";
289296
break;
@@ -309,14 +316,16 @@ void SteppingAction::PrintKilledTrack(const G4Track* aTrack, const TrackStatus&
309316
break;
310317
}
311318
G4VPhysicalVolume* pv = aTrack->GetNextVolume();
312-
if (pv) {
313-
vname = pv->GetLogicalVolume()->GetName();
314-
rname = pv->GetLogicalVolume()->GetRegion()->GetName();
315-
}
319+
vname = pv->GetLogicalVolume()->GetName();
320+
rname = pv->GetLogicalVolume()->GetRegion()->GetName();
316321

317-
edm::LogVerbatim("SimG4CoreApplication")
322+
const double ekin = aTrack->GetKineticEnergy();
323+
if (ekin < 2 * CLHEP::MeV) {
324+
return;
325+
}
326+
edm::LogWarning("SimG4CoreApplication")
318327
<< "Track #" << aTrack->GetTrackID() << " StepN= " << aTrack->GetCurrentStepNumber() << " "
319-
<< aTrack->GetDefinition()->GetParticleName() << " E(MeV)= " << aTrack->GetKineticEnergy() / MeV
320-
<< " T(ns)= " << aTrack->GetGlobalTime() / ns << " is killed due to " << typ << " inside LV: " << vname << " ("
321-
<< rname << ") at " << aTrack->GetPosition();
328+
<< aTrack->GetDefinition()->GetParticleName() << " E(MeV)=" << ekin / CLHEP::MeV
329+
<< " T(ns)=" << aTrack->GetGlobalTime() / CLHEP::ns << " is killed due to " << typ << "\n LV: " << vname << " ("
330+
<< rname << ") at " << aTrack->GetPosition() << " step(cm)=" << aTrack->GetStep()->GetStepLength() / CLHEP::cm;
322331
}

0 commit comments

Comments
 (0)