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.

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

Fix energy deposition for killed looping positrons #708

Merged
merged 2 commits into from
Mar 30, 2023
Merged

Fix energy deposition for killed looping positrons #708

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
44773c8
Add a flag for antiparticle
amandalund Mar 29, 2023
c7d6f80
Fix energy deposition for killed looping positrons
amandalund Mar 29, 2023
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
8 changes: 7 additions & 1 deletion src/celeritas/global/alongstep/AlongStep.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,8 +85,14 @@ inline CELER_FUNCTION void along_step(MH&& msc,
{
// If the track is looping (or if it's a stuck track that waa
// flagged as looping), deposit the energy locally.
auto deposited = particle.energy().value();
if (particle.is_antiparticle())
{
// Energy conservation for killed positrons
deposited += 2 * particle.mass().value();
}
track.make_physics_step_view().deposit_energy(
particle.energy());
ParticleTrackView::Energy{deposited});
particle.subtract_energy(particle.energy());

// Mark that this track was abandoned while looping
Expand Down
1 change: 1 addition & 0 deletions src/celeritas/phys/ParticleData.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ struct ParticleRecord
units::MevMass mass; //!< Rest mass [MeV / c^2]
units::ElementaryCharge charge; //!< Charge in units of [e]
real_type decay_constant; //!< Decay constant [1/s]
bool is_antiparticle; //!< Antiparticle (negative PDG number)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably a basic/dumb physics question: antiprotons and antielectrons obviously annihilate with regular matter. Do unstable antiparticles decay entirely to stable antiparticles? If not then it's not ideal to dump twice the rest mass into local energy deposition in all cases...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since only stable particles can be killed we should be safe either way.. but maybe since we really want to check if it's a positron we should just store the positron ID somewhere?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, I think this is certainly good enough for the time being since the only antiparticle we have is positron 🙃. The code in the loop is clear enough that when we add other antiparticles we'll know to revisit it.

🙃 + 🙂 → 💥💥


//! Value of decay_constant for a stable particle
static CELER_CONSTEXPR_FUNCTION real_type stable_decay_constant()
Expand Down
1 change: 1 addition & 0 deletions src/celeritas/phys/ParticleParams.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,7 @@ ParticleParams::ParticleParams(Input const& input)
host_def.mass = particle.mass;
host_def.charge = particle.charge;
host_def.decay_constant = particle.decay_constant;
host_def.is_antiparticle = particle.pdg_code.get() < 0;
particles.push_back(std::move(host_def));
}

Expand Down
12 changes: 12 additions & 0 deletions src/celeritas/phys/ParticleTrackView.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,9 @@ class ParticleTrackView
// Decay constant [1/s]
CELER_FORCEINLINE_FUNCTION real_type decay_constant() const;

// Whether it is an antiparticle
CELER_FORCEINLINE_FUNCTION bool is_antiparticle() const;

// Whether the particle is stable
CELER_FORCEINLINE_FUNCTION bool is_stable() const;

Expand Down Expand Up @@ -232,6 +235,15 @@ CELER_FUNCTION real_type ParticleTrackView::decay_constant() const
return this->particle_view().decay_constant();
}

//---------------------------------------------------------------------------//
/*!
* Whether it is an antiparticle.
*/
CELER_FUNCTION bool ParticleTrackView::is_antiparticle() const
{
return this->particle_view().is_antiparticle();
}

//---------------------------------------------------------------------------//
/*!
* Whether particle is stable.
Expand Down
12 changes: 12 additions & 0 deletions src/celeritas/phys/ParticleView.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,9 @@ class ParticleView
// Decay constant [1/s]
CELER_FORCEINLINE_FUNCTION real_type decay_constant() const;

// Whether it is an antiparticle
CELER_FORCEINLINE_FUNCTION bool is_antiparticle() const;

private:
ParticleParamsRef const& params_;
const ParticleId particle_;
Expand Down Expand Up @@ -96,5 +99,14 @@ CELER_FUNCTION real_type ParticleView::decay_constant() const
return params_.particles[particle_].decay_constant;
}

//---------------------------------------------------------------------------//
/*!
* Whether it is an antiparticle.
*/
CELER_FUNCTION bool ParticleView::is_antiparticle() const
{
return params_.particles[particle_].is_antiparticle;
}

//---------------------------------------------------------------------------//
} // namespace celeritas
26 changes: 25 additions & 1 deletion test/celeritas/phys/Particle.test.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,11 @@ class ParticleTest : public Test
MevMass{939.565413},
zero,
1.0 / (879.4 * second)});
defs.push_back({"positron",
pdg::positron(),
MevMass{0.5109989461},
ElementaryCharge{1},
stable});

particle_params = std::make_shared<ParticleParams>(std::move(defs));
}
Expand All @@ -69,10 +74,12 @@ TEST_F(ParticleTest, params_accessors)
EXPECT_EQ(ParticleId(0), defs.find(PDGNumber(11)));
EXPECT_EQ(ParticleId(1), defs.find(PDGNumber(22)));
EXPECT_EQ(ParticleId(2), defs.find(PDGNumber(2112)));
EXPECT_EQ(ParticleId(3), defs.find(PDGNumber(-11)));

EXPECT_EQ(ParticleId(0), defs.find("electron"));
EXPECT_EQ(ParticleId(1), defs.find("gamma"));
EXPECT_EQ(ParticleId(2), defs.find("neutron"));
EXPECT_EQ(ParticleId(3), defs.find("positron"));

EXPECT_EQ("electron", defs.id_to_label(ParticleId(0)));
EXPECT_EQ(PDGNumber(11), defs.id_to_pdg(ParticleId(0)));
Expand All @@ -86,7 +93,7 @@ TEST_F(ParticleTest, output)
if (CELERITAS_USE_JSON)
{
EXPECT_EQ(
R"json([{"label":"electron","pdg":11},{"label":"gamma","pdg":22},{"label":"neutron","pdg":2112}])json",
R"json([{"label":"electron","pdg":11},{"label":"gamma","pdg":22},{"label":"neutron","pdg":2112},{"label":"positron","pdg":-11}])json",
to_string(out));
}
}
Expand Down Expand Up @@ -171,6 +178,7 @@ TEST_F(ParticleTestHost, electron)
EXPECT_DOUBLE_EQ(0.5109989461, particle.mass().value());
EXPECT_DOUBLE_EQ(-1., particle.charge().value());
EXPECT_DOUBLE_EQ(0.0, particle.decay_constant());
EXPECT_FALSE(particle.is_antiparticle());
EXPECT_TRUE(particle.is_stable());
EXPECT_SOFT_EQ(0.74453076757415848, particle.beta_sq());
EXPECT_SOFT_EQ(0.86286196322132447, particle.speed().value());
Expand All @@ -188,6 +196,20 @@ TEST_F(ParticleTestHost, electron)
EXPECT_DOUBLE_EQ(0.0, particle.energy().value());
}

TEST_F(ParticleTestHost, positron)
{
ParticleTrackView particle(
particle_params->host_ref(), state_ref, TrackSlotId(0));
particle = Initializer_t{ParticleId{3}, MevEnergy{1}};

EXPECT_DOUBLE_EQ(1, particle.energy().value());
EXPECT_DOUBLE_EQ(0.5109989461, particle.mass().value());
EXPECT_DOUBLE_EQ(1., particle.charge().value());
EXPECT_DOUBLE_EQ(0.0, particle.decay_constant());
EXPECT_TRUE(particle.is_antiparticle());
EXPECT_TRUE(particle.is_stable());
}

TEST_F(ParticleTestHost, gamma)
{
ParticleTrackView particle(
Expand All @@ -196,6 +218,7 @@ TEST_F(ParticleTestHost, gamma)

EXPECT_DOUBLE_EQ(0, particle.mass().value());
EXPECT_DOUBLE_EQ(10, particle.energy().value());
EXPECT_FALSE(particle.is_antiparticle());
EXPECT_TRUE(particle.is_stable());
EXPECT_DOUBLE_EQ(1.0, particle.beta_sq());
EXPECT_DOUBLE_EQ(1.0, particle.speed().value());
Expand All @@ -210,6 +233,7 @@ TEST_F(ParticleTestHost, neutron)

EXPECT_DOUBLE_EQ(20, particle.energy().value());
EXPECT_DOUBLE_EQ(1.0 / 879.4, particle.decay_constant());
EXPECT_FALSE(particle.is_antiparticle());
EXPECT_FALSE(particle.is_stable());
}

Expand Down