Add MuBremsstrahlungInteractor #257

DoaaDeeb · 2021-06-21T19:14:19Z

This commit includes the following new files:

src/physics/em/BetheBlochModel.cc
src/physics/em/BetheBlochModel.hh
src/physics/em/detail/BetheBloch.cu
src/physics/em/detail/BetheBloch.hh
src/physics/em/detail/BetheBlochInteractor.hh
src/physics/em/detail/BetheBlochInteractor.i.hh
test/physics/em/BetheBloch.test.cc

and it updates the following files:

src/CMakeLists.txt
test/CMakeLists.txt

This commit includes the following new files: src/physics/em/BetheBlochModel.cc src/physics/em/BetheBlochModel.hh src/physics/em/detail/BetheBloch.cu src/physics/em/detail/BetheBloch.hh src/physics/em/detail/BetheBlochInteractor.hh src/physics/em/detail/BetheBlochInteractor.i.hh test/physics/em/BetheBloch.test.cc and it updates the following files: src/CMakeLists.txt test/CMakeLists.txt

src/physics/em/detail/BetheBloch.cu

src/physics/em/detail/BetheBlochInteractor.i.hh

src/physics/em/BetheBlochModel.hh

src/physics/em/detail/BetheBlochInteractor.hh

amandalund

Great job @DoaaDeeb on our first muon model! From the implementation and comments it looks like this is for bremsstrahlung, not ionization as the name BetheBloch suggests, correct? So we should probably rename to something like MuBremsstrahlung. Here's some initial feedback:

src/physics/em/detail/BetheBlochInteractor.i.hh

rename the classes to MuBremsstrahlung add Euler's number to src/base/Constants.hh use clamp_to_nonneg instead of using if statement pass an element to the interactor constructor instead of material use electron mass in MeV/c^2 fix some calculations

DoaaDeeb · 2021-06-24T01:44:06Z

Great job @DoaaDeeb on our first muon model! From the implementation and comments it looks like this is for bremsstrahlung, not ionization as the name BetheBloch suggests, correct? So we should probably rename to something like MuBremsstrahlung. Here's some initial feedback:

Thank you @amandalund for your feedback. I renamed the classes to MuBremsstrahlung as you suggested.

amandalund

Thanks @DoaaDeeb, looking good! One more round of mostly minor comments from my end.

src/physics/em/detail/MuBremsstrahlung.hh

src/physics/em/detail/MuBremsstrahlungInteractor.hh

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh

src/physics/em/detail/MuBremsstrahlungInteractor.hh

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh

test/physics/em/MuBremsstrahlung.test.cc

src/physics/em/MuBremsstrahlungModel.cc

src/physics/em/detail/MuBremsstrahlungInteractor.hh

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh

stognini · 2021-06-25T16:49:34Z

Great work @DoaaDeeb. I have two general comments that can be applied to the overall interactor -- I can help you find missing bits or discuss them further on Slack if you would like to:

There are many instances where a real_type should be set as a const real_type. So I suggest double checking the variables and set them as const whenever it makes sense.
The variable naming is very close to Geant4 ones, and this helps comparing, but at the same time it breaks our current code standards. So my suggestion is to not feel afraid of using longer variable naming for the sake of clarity. For example, using atomic_number instead of Z, inc_total_energy instead of E, inc_mass or muon_mass instead of mass, and so on. I believe code naming standards and clarity are better than using G4's nomenclature. Again, if you would like to discuss or need help on deciding variable naming, feel free to ping me on slack.

sethrj · 2021-06-26T12:32:50Z

Oh noooo, the CI is having this error (which isn't your fault):

nvlink error   : Size doesn't match for '_ZN9celeritas9SecondaryC1Ev$172' in '/var/jenkins/workspace/Celeritas_PR-257/build/src/libceleritas_static.a:MuBremsstrahlung.cu.o', first specified in '/var/jenkins/workspace/Celeritas_PR-257/build/src/libceleritas_static.a:BetheHeitler.cu.o'

same as #118

amandalund · 2021-06-27T20:06:01Z

Oh noooo, the CI is having this error (which isn't your fault):

nvlink error   : Size doesn't match for '_ZN9celeritas9SecondaryC1Ev$172' in '/var/jenkins/workspace/Celeritas_PR-257/build/src/libceleritas_static.a:MuBremsstrahlung.cu.o', first specified in '/var/jenkins/workspace/Celeritas_PR-257/build/src/libceleritas_static.a:BetheHeitler.cu.o'

same as #118

Backing out the changes, I stop getting the nvlink error when the interactor takes MaterialView as an argument instead of ElementView. You could try this workaround, which at least on emmet is working for me:

diff --git a/src/physics/em/detail/MuBremsstrahlung.cu b/src/physics/em/detail/MuBremsstrahlung.cu
index 0b31052..e1750c0 100644
--- a/src/physics/em/detail/MuBremsstrahlung.cu
+++ b/src/physics/em/detail/MuBremsstrahlung.cu
@@ -59,13 +59,14 @@ __global__ void mu_bremsstrahlung_interact_kernel(
     if (physics.model_id() != mb.model_id)
         return;
 
-    ElementView element
-        = material_view.element_view(celeritas::ElementComponentId{0});
+    // TODO: sample an element. For now assume one element per material
+    const ElementComponentId elcomp_id{0};
     MuBremsstrahlungInteractor interact(mb,
                                         particle,
                                         model.states.direction[tid],
                                         allocate_secondaries,
-                                        element);
+                                        material_view,
+                                        elcomp_id);
 
     RngEngine rng(model.states.rng, tid);
     model.states.interactions[tid] = interact(rng);
diff --git a/src/physics/em/detail/MuBremsstrahlungInteractor.hh b/src/physics/em/detail/MuBremsstrahlungInteractor.hh
index 335899d..3868399 100644
--- a/src/physics/em/detail/MuBremsstrahlungInteractor.hh
+++ b/src/physics/em/detail/MuBremsstrahlungInteractor.hh
@@ -15,6 +15,7 @@
 #include "physics/base/Secondary.hh"
 #include "physics/base/Units.hh"
 #include "physics/material/ElementView.hh"
+#include "physics/material/MaterialView.hh"
 #include "MuBremsstrahlung.hh"
 
 namespace celeritas
@@ -38,7 +39,8 @@ class MuBremsstrahlungInteractor
         const ParticleTrackView&                  particle,
         const Real3&                              inc_direction,
         StackAllocator<Secondary>&                allocate,
-        ElementView&                              element);
+        const MaterialView&                       material,
+        ElementComponentId                        elcomp_id);
 
     // Sample an interaction with the given RNG
     template<class Engine>
@@ -74,7 +76,7 @@ class MuBremsstrahlungInteractor
     // Allocate space for one or more secondary particles
     StackAllocator<Secondary>& allocate_;
     // Element properties
-    const ElementView& element_;
+    const ElementView element_;
     // Incident muon mass
     const units::MevMass inc_mass_;
 };
diff --git a/src/physics/em/detail/MuBremsstrahlungInteractor.i.hh b/src/physics/em/detail/MuBremsstrahlungInteractor.i.hh
index 15cdb8a..901aab6 100644
--- a/src/physics/em/detail/MuBremsstrahlungInteractor.i.hh
+++ b/src/physics/em/detail/MuBremsstrahlungInteractor.i.hh
@@ -23,12 +23,13 @@ CELER_FUNCTION MuBremsstrahlungInteractor::MuBremsstrahlungInteractor(
     const ParticleTrackView&                  particle,
     const Real3&                              inc_direction,
     StackAllocator<Secondary>&                allocate,
-    ElementView&                              element)
+    const MaterialView&                       material,
+    ElementComponentId                        elcomp_id)
     : shared_(shared)
     , inc_energy_(particle.energy().value())
     , inc_direction_(inc_direction)
     , allocate_(allocate)
-    , element_(element)
+    , element_(material.element_view(elcomp_id))
     , inc_mass_(particle.mass().value())
 {
     CELER_EXPECT(inc_energy_ >= this->min_incident_energy()
diff --git a/test/physics/em/MuBremsstrahlung.test.cc b/test/physics/em/MuBremsstrahlung.test.cc
index 3b25730..043fbbb 100644
--- a/test/physics/em/MuBremsstrahlung.test.cc
+++ b/test/physics/em/MuBremsstrahlung.test.cc
@@ -112,16 +112,15 @@ TEST_F(MuBremsstrahlungInteractorTest, basic)
     int num_samples = 4;
     this->resize_secondaries(num_samples);
 
-    celeritas::ElementView element(
-        this->material_track().material_view().element_view(
-            celeritas::ElementComponentId{0}));
+    auto material = this->material_track().material_view();
 
     // Create the interactor
     MuBremsstrahlungInteractor interact(pointers_,
                                         this->particle_track(),
                                         this->direction(),
                                         this->secondary_allocator(),
-                                        element);
+                                        material,
+                                        celeritas::ElementComponentId{0});
     RandomEngine&              rng_engine = this->rng();
 
     std::vector<double> energy;
@@ -185,16 +184,15 @@ TEST_F(MuBremsstrahlungInteractorTest, stress_test)
                 this->set_inc_direction(inc_dir);
                 this->resize_secondaries(num_samples);
 
-                celeritas::ElementView element(
-                    this->material_track().material_view().element_view(
-                        celeritas::ElementComponentId{0}));
+                auto material = this->material_track().material_view();
 
                 // Create interactor
                 MuBremsstrahlungInteractor interact(pointers_,
                                                     this->particle_track(),
                                                     this->direction(),
                                                     this->secondary_allocator(),
-                                                    element);
+                                                    material,
+                                                    celeritas::ElementComponentId{0});
 
                 for (unsigned int i = 0; i < num_samples; i++)
                 {

move min and max energy functions to MuBremsstrahlungInteractorPointers fix energy limits pass MaterialView and ElementComponentId to the interactor constructor rename variables to have descriptive names set variables as a const when suitable remove duplicate variables assign electron_mass in the test setup set sample_num to 10k in the stress test add comment separators

DoaaDeeb · 2021-06-28T20:10:00Z

Great work @DoaaDeeb. I have two general comments that can be applied to the overall interactor -- I can help you find missing bits or discuss them further on Slack if you would like to:

There are many instances where a real_type should be set as a const real_type. So I suggest double checking the variables and set them as const whenever it makes sense.

The variable naming is very close to Geant4 ones, and this helps comparing, but at the same time it breaks our current code standards. So my suggestion is to not feel afraid of using longer variable naming for the sake of clarity. For example, using atomic_number instead of Z, inc_total_energy instead of E, inc_mass or muon_mass instead of mass, and so on. I believe code naming standards and clarity are better than using G4's nomenclature. Again, if you would like to discuss or need help on deciding variable naming, feel free to ping me on slack.

Thank you, @stognini, for your help. I changed the variable names as you suggested.

amandalund · 2021-06-28T21:04:46Z

Excellent work @DoaaDeeb! ~~Looks like with the new changes we need another workaround for that nvlink error, so let's try this:~~

diff --git a/src/physics/em/detail/MuBremsstrahlungInteractor.i.hh b/src/physics/em/detail/MuBremsstrahlungInteractor.i.hh
index 1c9d268..7233d5f 100644
--- a/src/physics/em/detail/MuBremsstrahlungInteractor.i.hh
+++ b/src/physics/em/detail/MuBremsstrahlungInteractor.i.hh
@@ -8,7 +8,6 @@
 
 #include "base/ArrayUtils.hh"
 #include "base/Constants.hh"
-#include "base/Range.hh"
 #include "random/distributions/UniformRealDistribution.hh"
 
 namespace celeritas
@@ -81,7 +80,7 @@ CELER_FUNCTION Interaction MuBremsstrahlungInteractor::operator()(Engine& rng)
     const real_type tot_momentum = std::sqrt(inc_energy_.value() *
                                 (inc_energy_.value() + 2 * inc_mass_.value()));
     Real3 inc_direction;
-    for (int i : range(3))
+    for (int i = 0; i < 3; ++i)
     {
         inc_direction[i] = tot_momentum * inc_direction_[i]
                                         - epsilon * gamma_dir[i];

stognini · 2021-06-28T21:49:00Z

test/physics/em/MuBremsstrahlung.test.cc

+
+    for (auto particle : {pdg::mu_minus(), pdg::mu_plus()})
+    {
+        for (double inc_e : {1.5e4, 5e4, 10e4, 50e4, 100e4})


Does the value of avg_engine_samples becomes too high to be reasonable as it approaches the model limits (1e3 or 1e7 MeV)? It's important to know if we need to worry about having a never ending sampling loop during the interactor's do {...} while (...); in edge cases.

I tried values that approach the model limits, and it seems that avg_engine_samples are still reasonable.

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh

amandalund · 2021-06-29T05:36:47Z

Alright... so ignore that last workaround, and try removing member initialization from Secondary instead:

diff --git a/src/physics/base/Secondary.hh b/src/physics/base/Secondary.hh
index 69d1af4..1d2cf2c 100644
--- a/src/physics/base/Secondary.hh
+++ b/src/physics/base/Secondary.hh
@@ -21,9 +21,9 @@ namespace celeritas
  */
 struct Secondary
 {
-    ParticleId       particle_id{};           //!< New particle type
-    units::MevEnergy energy{zero_quantity()}; //!< New kinetic energy
-    Real3            direction;               //!< New direction
+    ParticleId       particle_id; //!< New particle type
+    units::MevEnergy energy;      //!< New kinetic energy
+    Real3            direction;   //!< New direction
 
     //! Whether the secondary survived cutoffs
     explicit CELER_FUNCTION operator bool() const

With this change I was able to build successfully in each of the cases we ran into that nvlink error previously.

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh

simplify the rejection loop in operator() remove member initialization from Secondary.hh remove tot-momentum and add inc_momentum_ as memeber variable instead

amandalund

Looks like you'll need to apply clang-format one more time, but otherwise looks good on my end. Great job @DoaaDeeb, and thanks for your patience with the back and forth as we worked around that nvlink bug!

DoaaDeeb · 2021-07-05T21:37:32Z

Looks like you'll need to apply clang-format one more time, but otherwise looks good on my end. Great job @DoaaDeeb, and thanks for your patience with the back and forth as we worked around that nvlink bug!

Thank you, @amandalund. I appreciate your help in fixing the bug.

sethrj

Thanks for the excellent work, @DoaaDeeb ! The only real issue for now is the skimpy documentation, though I understand most of the calculations are adopted from the Geant reference manuals or code. I don't think you need to write out every equation, but it would be good if you could reference specific equation numbers in the sampling functions.

src/physics/em/detail/MuBremsstrahlungInteractor.hh

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh

add more documentation to MuBremsstrahlungInteractor.hh store particle_ as a member variable remove inc_momentum_, inc_mass_, inc_energy_ and use particle_ methods use generate_canonical

DoaaDeeb · 2021-07-12T19:15:13Z

Thanks for the excellent work, @DoaaDeeb ! The only real issue for now is the skimpy documentation, though I understand most of the calculations are adopted from the Geant reference manuals or code. I don't think you need to write out every equation, but it would be good if you could reference specific equation numbers in the sampling functions.

Thanks @sethrj for your feedback. I noticed that the equations used in the code are not numbered in the reference manual. Should I write them out instead?

amandalund · 2021-07-12T19:55:20Z

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh

-        d_n_prime = d_n / std::pow(d_n, 1.0 / atomic_number);
-        b         = 183.0;
-        b1        = 1429.0;
+        d_n_prime = d_n / std::pow(d_n, 1 / atomic_number);


An integer division snuck in (1 / atomic_number will always be 0).

Also, I see the usage of both std::pow() and ipow<n>(). Aside from consistency, do we have to worry if they're used interchangeably?

There it is!! Thanks amanda, I was poring over the changeset for 15 minutes to find what changed before seeing your comment. Thanks!

You could also simplify the expression slightly while you're at it,

Suggested change

d_n_prime = d_n / std::pow(d_n, 1 / atomic_number);

d_n_prime = std::pow(d_n, 1 - real_type(1) / atomic_number);

Also, I see the usage of both std::pow() and ipow<n>(). Aside from consistency, do we have to worry if they're used interchangeably?

ipow should always be used for known integer n; otherwise you have to use std::pow (or some other runtime expression).

stognini · 2021-07-12T22:12:38Z

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh

+        b         = real_type(202.4);
+        b1        = 446;


I see this in a few other spots too, where some real_type, such as b, is explicitly casted again, while other numbers declared as real_type aren't, such as b1. Is there any reason to do so?

I think either way is technically fine, though it would be probably better to do one or the other. Interestingly in your example the first one could be a demotion (double -> float) while the second would be promotion (int -> float)...

sethrj · 2021-07-12T22:29:49Z

Thanks @sethrj for your feedback.

Thanks for addressing the comments! Looks good aside from the accidental integer division Amanda pointed out. (That's the danger of my advice of using integers and relying on floating point promotion: it's easier to get a completely wrong answer, as opposed to accidentally making your code slower and getting slightly different results.)

I noticed that the equations used in the code are not numbered in the reference manual. Should I write them out instead?

That's a pity, I'd say not to worry about it in that case, since you already referenced the section in the manual.

sethrj

Impressive work! Thanks for the final fix.

stognini reviewed Jun 21, 2021

View reviewed changes

src/physics/em/detail/BetheBloch.cu Outdated Show resolved Hide resolved

stognini reviewed Jun 21, 2021

View reviewed changes

src/physics/em/detail/BetheBloch.cu Outdated Show resolved Hide resolved

stognini reviewed Jun 21, 2021

View reviewed changes

src/physics/em/detail/BetheBlochInteractor.i.hh Outdated Show resolved Hide resolved

stognini reviewed Jun 21, 2021

View reviewed changes

src/physics/em/BetheBlochModel.hh Outdated Show resolved Hide resolved

stognini reviewed Jun 21, 2021

View reviewed changes

src/physics/em/detail/BetheBlochInteractor.hh Outdated Show resolved Hide resolved

amandalund requested changes Jun 22, 2021

View reviewed changes

DoaaDeeb changed the title ~~Add BetheBlochInteractor~~ Add MuBremsstrahlungInteractor Jun 24, 2021

amandalund requested changes Jun 24, 2021

View reviewed changes

stognini reviewed Jun 25, 2021

View reviewed changes

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh Show resolved Hide resolved

stognini reviewed Jun 25, 2021

View reviewed changes

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh Show resolved Hide resolved

stognini reviewed Jun 25, 2021

View reviewed changes

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh Outdated Show resolved Hide resolved

sethrj added the physics Particles, processes, and stepping algorithms label Jun 26, 2021

sethrj added this to the Y2 Goals milestone Jun 26, 2021

stognini reviewed Jun 28, 2021

View reviewed changes

amandalund reviewed Jun 29, 2021

View reviewed changes

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh Outdated Show resolved Hide resolved

amandalund reviewed Jun 29, 2021

View reviewed changes

src/physics/em/detail/MuBremsstrahlungInteractor.i.hh Outdated Show resolved Hide resolved

More updates:

85b7860

simplify the rejection loop in operator() remove member initialization from Secondary.hh remove tot-momentum and add inc_momentum_ as memeber variable instead

amandalund approved these changes Jul 2, 2021

View reviewed changes

apply clang-format

bbf6b62

sethrj requested changes Jul 7, 2021

View reviewed changes

Updates:

676bd40

add more documentation to MuBremsstrahlungInteractor.hh store particle_ as a member variable remove inc_momentum_, inc_mass_, inc_energy_ and use particle_ methods use generate_canonical

amandalund reviewed Jul 12, 2021

View reviewed changes

stognini reviewed Jul 12, 2021

View reviewed changes

fix integer division

2989193

sethrj approved these changes Jul 16, 2021

View reviewed changes

sethrj merged commit 5a33d32 into celeritas-project:master Jul 16, 2021

sethrj mentioned this pull request Jul 19, 2021

Add muon physics #180

Closed

11 tasks

sethrj modified the milestones: Y2Q4 Goals, Y2Q3 Goals, Y2Q2 Goals Sep 24, 2021

sethrj added the enhancement New feature or request label Feb 18, 2023

sethrj mentioned this pull request Jul 11, 2024

Fix muon bremsstrahlung energy limits and sampling #1319

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add MuBremsstrahlungInteractor #257

Add MuBremsstrahlungInteractor #257

DoaaDeeb commented Jun 21, 2021

amandalund left a comment

DoaaDeeb commented Jun 24, 2021 •

edited

Loading

amandalund left a comment

stognini commented Jun 25, 2021

sethrj commented Jun 26, 2021

amandalund commented Jun 27, 2021

DoaaDeeb commented Jun 28, 2021

amandalund commented Jun 28, 2021 •

edited

Loading

stognini Jun 28, 2021

DoaaDeeb Jul 1, 2021

amandalund commented Jun 29, 2021

amandalund left a comment

DoaaDeeb commented Jul 5, 2021

sethrj left a comment

DoaaDeeb commented Jul 12, 2021

amandalund Jul 12, 2021

stognini Jul 12, 2021

sethrj Jul 12, 2021

sethrj Jul 12, 2021

stognini Jul 12, 2021

sethrj Jul 12, 2021

sethrj commented Jul 12, 2021

sethrj left a comment

	d_n_prime = d_n / std::pow(d_n, 1 / atomic_number);
	d_n_prime = std::pow(d_n, 1 - real_type(1) / atomic_number);

Add MuBremsstrahlungInteractor #257

Add MuBremsstrahlungInteractor #257

Conversation

DoaaDeeb commented Jun 21, 2021

amandalund left a comment

Choose a reason for hiding this comment

DoaaDeeb commented Jun 24, 2021 • edited Loading

amandalund left a comment

Choose a reason for hiding this comment

stognini commented Jun 25, 2021

sethrj commented Jun 26, 2021

amandalund commented Jun 27, 2021

DoaaDeeb commented Jun 28, 2021

amandalund commented Jun 28, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

amandalund commented Jun 29, 2021

amandalund left a comment

Choose a reason for hiding this comment

DoaaDeeb commented Jul 5, 2021

sethrj left a comment

Choose a reason for hiding this comment

DoaaDeeb commented Jul 12, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sethrj commented Jul 12, 2021

sethrj left a comment

Choose a reason for hiding this comment

DoaaDeeb commented Jun 24, 2021 •

edited

Loading

amandalund commented Jun 28, 2021 •

edited

Loading