-
Notifications
You must be signed in to change notification settings - Fork 8
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
Add Particle filter stacking action option #244
base: main
Are you sure you want to change the base?
Changes from all commits
2296c90
8625b44
52053e9
2c5c34f
da1e6f1
a69af07
63f9f67
43a0e9e
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,83 @@ | ||
// Copyright (C) 2022 Luigi Pertoldi <[email protected]> | ||
// | ||
// This program is free software: you can redistribute it and/or modify it under | ||
// the terms of the GNU Lesser General Public License as published by the Free | ||
// Software Foundation, either version 3 of the License, or (at your option) any | ||
// later version. | ||
// | ||
// This program is distributed in the hope that it will be useful, but WITHOUT | ||
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS | ||
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more | ||
// details. | ||
// | ||
// You should have received a copy of the GNU Lesser General Public License | ||
// along with this program. If not, see <https://www.gnu.org/licenses/>. | ||
|
||
#ifndef _RMG_PARTICLE_FILTER_OUTPUT_SCHEME_HH_ | ||
#define _RMG_PARTICLE_FILTER_OUTPUT_SCHEME_HH_ | ||
|
||
#include <optional> | ||
#include <set> | ||
|
||
#include "G4GenericMessenger.hh" | ||
|
||
#include "RMGVOutputScheme.hh" | ||
|
||
/** @brief Filter-output scheme for particles. | ||
* | ||
* @details This optional output scheme filters out all particles specified via | ||
* their PDG code. Optionally, it can apply this filter to a specified volume | ||
* or ignore a specified volume. Properties need to be specified per macro. | ||
* | ||
* Does nothing if the output scheme is not enabled per macro before run initialization. | ||
* Also does nothing if no particle is specified. If no volume is specified the filter | ||
* is applied to all volumes. | ||
*/ | ||
class RMGParticleFilterOutputScheme : public RMGVOutputScheme { | ||
|
||
public: | ||
|
||
RMGParticleFilterOutputScheme(); | ||
|
||
/** @brief Wraps @c G4UserStackingAction::StackingActionClassify | ||
* @details This is used to classify all specified particles as @c fKill if | ||
* they are in the specified volumes (or no volume is specified). | ||
* | ||
* If the primary particle is filtered out here, the simulation crashes. To | ||
* avoid the crash, the particle will be simulated anyways and a warning message will | ||
* be shown. | ||
*/ | ||
std::optional<G4ClassificationOfNewTrack> StackingActionClassify(const G4Track*, int) override; | ||
|
||
/** @brief Add a particle, identified by its PDG code, to the list of particles to kill. */ | ||
inline void AddParticle(int pdg) { fParticles.insert(pdg); } | ||
|
||
/** @brief Add a physical volume, by name, to the volumes in which the filter will | ||
* not be applied. | ||
* @details This means that specified particles outside of the specified | ||
* keep-volumes will be filtered. It is therefore not possible to specify keep-volumes | ||
* and kill-volumes in the same geometry. | ||
*/ | ||
void AddKeepVolume(std::string name); | ||
|
||
/** @brief Add a physical volume, by name, to the volumes in which the filter will | ||
* be applied. | ||
* @details This means that specified particles outside of the specified | ||
* kill-volumes will not be affected by the filter. It is therefore not possible to | ||
* specify keep-volumes and kill-volumes in the same geometry. | ||
*/ | ||
void AddKillVolume(std::string name); | ||
|
||
private: | ||
|
||
std::unique_ptr<G4GenericMessenger> fMessenger; | ||
void DefineCommands(); | ||
|
||
std::set<int> fParticles; | ||
std::set<std::string> fKeepVolumes; | ||
std::set<std::string> fKillVolumes; | ||
}; | ||
|
||
#endif | ||
|
||
// vim: tabstop=2 shiftwidth=2 expandtab |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,104 @@ | ||
// Copyright (C) 2022 Luigi Pertoldi <[email protected]> | ||
// | ||
// This program is free software: you can redistribute it and/or modify it under | ||
// the terms of the GNU Lesser General Public License as published by the Free | ||
// Software Foundation, either version 3 of the License, or (at your option) any | ||
// later version. | ||
// | ||
// This program is distributed in the hope that it will be useful, but WITHOUT | ||
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS | ||
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more | ||
// details. | ||
// | ||
// You should have received a copy of the GNU Lesser General Public License | ||
// along with this program. If not, see <https://www.gnu.org/licenses/>. | ||
|
||
#include "RMGParticleFilterOutputScheme.hh" | ||
|
||
#include <set> | ||
|
||
#include "G4Event.hh" | ||
#include "G4EventManager.hh" | ||
|
||
#include "RMGLog.hh" | ||
|
||
RMGParticleFilterOutputScheme::RMGParticleFilterOutputScheme() { this->DefineCommands(); } | ||
|
||
void RMGParticleFilterOutputScheme::AddKeepVolume(std::string name) { | ||
if (!fKillVolumes.empty()) { | ||
RMGLog::OutDev(RMGLog::fatal, "Conflicting requests for kill/keep volume in ParticleFilter. " | ||
"Trying to assign keep-volume but a kill-volume already exists."); | ||
} | ||
|
||
fKeepVolumes.insert(name); | ||
} | ||
|
||
void RMGParticleFilterOutputScheme::AddKillVolume(std::string name) { | ||
if (!fKeepVolumes.empty()) { | ||
RMGLog::OutDev(RMGLog::fatal, "Conflicting requests for kill/keep volume in ParticleFilter. " | ||
"Trying to assign kill-volume but a keep-volume already exists."); | ||
} | ||
|
||
fKillVolumes.insert(name); | ||
} | ||
|
||
std::optional<G4ClassificationOfNewTrack> RMGParticleFilterOutputScheme:: | ||
StackingActionClassify(const G4Track* aTrack, int stage) { | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Do we want to only stack particles after stage 1, eg. I am thinking about stacking nuclei only after they decay? Probably for nuclei its a bit difficult since you might have additional decays. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I am fairly certain the stages are something that is inherently handled by the user (us) and not G4? So far we run one stage, meaning everything until the Particle-Stack is empty. Then we check if any OutputScheme requests a new stage. If one OutputScheme requests a new stage, all tracks that have been marked with |
||
|
||
const int pdg = aTrack->GetDefinition()->GetPDGEncoding(); | ||
// If the particle is not marked to kill, let it go | ||
if (fParticles.find(pdg) == fParticles.end()) return std::nullopt; | ||
|
||
const auto pv = aTrack->GetTouchableHandle()->GetVolume(); | ||
// I am not 100% sure the physical volume will exist at this point. Can remove after some tests | ||
if (!pv) { RMGLog::OutDev(RMGLog::debug, "Physical volume does not exist in ParticleFilter!"); } | ||
const auto pv_name = pv->GetName(); | ||
// If a kill volume is specified only kill if in the kill volume. | ||
if (!fKillVolumes.empty() && fKillVolumes.find(pv_name) == fKillVolumes.end()) | ||
return std::nullopt; | ||
// If a keep volume is specified only keep particle if in that volume. | ||
if (fKeepVolumes.find(pv_name) != fKeepVolumes.end()) return std::nullopt; | ||
|
||
// If this is the primary particle we can not kill it without crashing the simulation | ||
if (aTrack->GetTrackID() == 0) { | ||
RMGLog::OutDev(RMGLog::warning, | ||
"ParticleFilter: The primary particle was specified to kill. " | ||
"This would cause the simulation to crash... so simulating it anyways. "); | ||
return std::nullopt; | ||
} | ||
|
||
// We land here if | ||
// i) Particle is marked to kill. | ||
// ii) No Kill volume specified or the particle is in the kill volume. | ||
// iii) Particle is not in the keep volume. | ||
// iiii) Particle is not the primary particle. | ||
RMGLog::OutDev(RMGLog::debug, "Filtering out particle with PDG code ", pdg, | ||
" in RMGParticleFilterOutputScheme"); | ||
return fKill; | ||
} | ||
|
||
void RMGParticleFilterOutputScheme::DefineCommands() { | ||
|
||
fMessenger = std::make_unique<G4GenericMessenger>(this, "/RMG/Output/ParticleFilter/", | ||
"Commands for filtering particles out by PDG encoding."); | ||
|
||
fMessenger->DeclareMethod("AddParticle", &RMGParticleFilterOutputScheme::AddParticle) | ||
.SetGuidance("Add a particle to be filtered out by its PDG code. User is responsible for " | ||
"correct PDG code.") | ||
.SetParameterName(0, "PDGcode", false, false) | ||
.SetStates(G4State_Idle); | ||
|
||
fMessenger->DeclareMethod("AddKeepVolume", &RMGParticleFilterOutputScheme::AddKeepVolume) | ||
.SetGuidance("Add a physical volume by name in which all specified Particles will be kept. " | ||
"They will be killed everywhere else. Can NOT be mixed with KillVolumes.") | ||
.SetParameterName(0, "PhysicalVolumeName", false, false) | ||
.SetStates(G4State_Idle); | ||
|
||
fMessenger->DeclareMethod("AddKillVolume", &RMGParticleFilterOutputScheme::AddKillVolume) | ||
.SetGuidance("Add a physical volume by name in which all specified Particles will be killed. " | ||
"They will only be killed in this volume. Can NOT be mixed with KeepVolumes.") | ||
.SetParameterName(0, "PhysicalVolumeName", false, false) | ||
.SetStates(G4State_Idle); | ||
} | ||
|
||
// vim: tabstop=2 shiftwidth=2 expandtab |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,19 @@ | ||
# collect auxiliary files | ||
file( | ||
GLOB _aux | ||
RELATIVE ${PROJECT_SOURCE_DIR} | ||
macros/* gdml/*.gdml gdml/*.xml run-test-*.sh run-test-*.py) | ||
|
||
# copy them to the build area | ||
foreach(_file ${_aux}) | ||
configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY) | ||
endforeach() | ||
|
||
add_test(NAME particlefilter/basicfilter COMMAND ${PYTHONPATH} run-test-basicfilter.py) | ||
set_tests_properties(particlefilter/basicfilter PROPERTIES LABELS extra) | ||
|
||
add_test(NAME particlefilter/killvolume COMMAND ${PYTHONPATH} run-test-killvolume.py) | ||
set_tests_properties(particlefilter/killvolume PROPERTIES LABELS extra) | ||
|
||
add_test(NAME particlefilter/keepvolume COMMAND ${PYTHONPATH} run-test-keepvolume.py) | ||
set_tests_properties(particlefilter/keepvolume PROPERTIES LABELS extra) |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,45 @@ | ||
<?xml version="1.0" ?> | ||
<gdml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd"> | ||
<define/> | ||
<materials/> | ||
<solids> | ||
<box name="world" x="2" y="2" z="2" lunit="m"/> | ||
<box name="scint" x="0.5" y="1" z="1" lunit="m"/> | ||
<box name="det" x="0.1" y="0.5" z="0.5" lunit="m"/> | ||
</solids> | ||
<structure> | ||
<volume name="det"> | ||
<materialref ref="G4_Ge"/> | ||
<solidref ref="det"/> | ||
</volume> | ||
<volume name="scint1"> | ||
<materialref ref="G4_lAr"/> | ||
<solidref ref="scint"/> | ||
<physvol name="det1"> | ||
<volumeref ref="det"/> | ||
</physvol> | ||
</volume> | ||
<volume name="scint2"> | ||
<materialref ref="G4_lAr"/> | ||
<solidref ref="scint"/> | ||
<physvol name="det2"> | ||
<volumeref ref="det"/> | ||
</physvol> | ||
</volume> | ||
<volume name="world"> | ||
<materialref ref="G4_Galactic"/> | ||
<solidref ref="world"/> | ||
<physvol name="scint1"> | ||
<volumeref ref="scint1"/> | ||
<position name="scint1_pos" x="-255.000000000000000" y="0.000000000000000" z="0.000000000000000" unit="mm"/> | ||
</physvol> | ||
<physvol name="scint2"> | ||
<volumeref ref="scint2"/> | ||
<position name="scint2_pos" x="255.000000000000000" y="0.000000000000000" z="0.000000000000000" unit="mm"/> | ||
</physvol> | ||
</volume> | ||
</structure> | ||
<setup name="Default" version="1.0"> | ||
<world ref="world"/> | ||
</setup> | ||
</gdml> |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,21 @@ | ||
/RMG/Geometry/GDMLDisableOverlapCheck true | ||
|
||
/RMG/Output/ActivateOutputScheme Track | ||
/RMG/Output/ActivateOutputScheme ParticleFilter | ||
|
||
/run/initialize | ||
|
||
/RMG/Output/ParticleFilter/AddParticle 1000020040 | ||
|
||
/RMG/Generator/Confine UnConfined | ||
|
||
/RMG/Generator/Select GPS | ||
/gps/pos/type Volume | ||
/gps/pos/shape Sphere | ||
/gps/pos/radius 0.5 m | ||
/gps/particle ion | ||
/gps/ion 90 228 | ||
/gps/energy 0 keV | ||
/gps/ang/type iso | ||
|
||
/run/beamOn 10 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
you own this code ;)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
:D i have not really cared much about the copyright stuff, i have mostly thought about this as enveloping the entire remage thingy and not that particular file.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would prefer the actual original author being credited. @tdixon97 @ManuelHu @MoritzNeuberger Please also sweep through the remage source files and check this.