diff --git a/src/celeritas/global/alongstep/AlongStep.hh b/src/celeritas/global/alongstep/AlongStep.hh index 5cd2b80404..85cbd27d3d 100644 --- a/src/celeritas/global/alongstep/AlongStep.hh +++ b/src/celeritas/global/alongstep/AlongStep.hh @@ -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 diff --git a/src/celeritas/phys/ParticleData.hh b/src/celeritas/phys/ParticleData.hh index 507c98e8b6..3fd22b737b 100644 --- a/src/celeritas/phys/ParticleData.hh +++ b/src/celeritas/phys/ParticleData.hh @@ -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) //! Value of decay_constant for a stable particle static CELER_CONSTEXPR_FUNCTION real_type stable_decay_constant() diff --git a/src/celeritas/phys/ParticleParams.cc b/src/celeritas/phys/ParticleParams.cc index 9d83e4842a..de4d01a91c 100644 --- a/src/celeritas/phys/ParticleParams.cc +++ b/src/celeritas/phys/ParticleParams.cc @@ -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)); } diff --git a/src/celeritas/phys/ParticleTrackView.hh b/src/celeritas/phys/ParticleTrackView.hh index 252f48a05c..e97afa724f 100644 --- a/src/celeritas/phys/ParticleTrackView.hh +++ b/src/celeritas/phys/ParticleTrackView.hh @@ -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; @@ -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. diff --git a/src/celeritas/phys/ParticleView.hh b/src/celeritas/phys/ParticleView.hh index 812619a610..f3414e3a29 100644 --- a/src/celeritas/phys/ParticleView.hh +++ b/src/celeritas/phys/ParticleView.hh @@ -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_; @@ -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 diff --git a/test/celeritas/phys/Particle.test.cc b/test/celeritas/phys/Particle.test.cc index c356e2be06..f9c4154c69 100644 --- a/test/celeritas/phys/Particle.test.cc +++ b/test/celeritas/phys/Particle.test.cc @@ -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(std::move(defs)); } @@ -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))); @@ -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)); } } @@ -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()); @@ -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( @@ -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()); @@ -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()); }