Merge pull request #50 from soleti/wion

Move from MEV2ELECTRONS to W_ION
DUNE · Dec 4, 2021 · 9914df5 · 9914df5
2 parents 9646d20 + fd2f446
commit 9914df5
Show file tree

Hide file tree

Showing 3 changed files with 12 additions and 15 deletions.
diff --git a/larndsim/consts/physics.py b/larndsim/consts/physics.py
@@ -13,10 +13,8 @@
 BIRKS_kb = 0.0486 # g/cm2/MeV Amoruso, et al NIM A 523 (2004) 275
 #: Electron charge in Coulomb
 E_CHARGE = 1.602e-19
-#: Average energy expended per ion pair in LAr in :math:`eV` from Phys. Rev. A 10, 1452
-W_ION = 23.6
-#: From MeV to electrons in LAr 
-MEV2ELECTRONS = 1/W_ION*1e6
+#: Average energy expended per ion pair in LAr in :math:`MeV` from Phys. Rev. A 10, 1452
+W_ION = 23.6e-6
 
 ## Quenching parameters
 BOX = 1

diff --git a/larndsim/quenching.py b/larndsim/quenching.py
@@ -7,7 +7,6 @@
 from numba import cuda
 
 from .consts import detector, physics, light
-from . import consts
 
 @cuda.jit
 def quench(tracks, mode):
@@ -41,5 +40,5 @@ def quench(tracks, mode):
         if isnan(recomb):
             raise RuntimeError("Invalid recombination value")
 
-        tracks[itrk]["n_electrons"] = recomb * dE * physics.MEV2ELECTRONS
+        tracks[itrk]["n_electrons"] = recomb * dE / physics.W_ION
         tracks[itrk]["n_photons"] = (dE/light.W_PH - tracks[itrk]["n_electrons"]) * light.SCINT_PRESCALE
diff --git a/tests/testQuenching.py b/tests/testQuenching.py
@@ -14,7 +14,7 @@ class TestQuenching:
     """
     #normal valid values
     tracks = np.zeros((100, 22))
-    tracks = np.core.records.fromarrays(tracks.transpose(), 
+    tracks = np.core.records.fromarrays(tracks.transpose(),
                                         names="eventID, dEdx, x_start, dE, t_start, z_end, trackID, x_end, y_end, n_electrons, n_photons, t, dx, pdgId, y, x, long_diff, z, z_start, y_start, tran_diff, t_end, pixel_plane",
                                         formats = "i8, f8, f8, f8, f8, f8, i8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, i8")
     tracks["dE"] = np.random.uniform(0.1, 100, 100)
@@ -23,14 +23,14 @@ class TestQuenching:
     #extreme valid values
     #a track with dEdx = 0, dE was set to 1 (any non-zero value) to test the recombination factor calculation
     track_zero = np.zeros((1,22))
-    track_zero = np.core.records.fromarrays(track_zero.transpose(), 
+    track_zero = np.core.records.fromarrays(track_zero.transpose(),
                                             names="eventID, dEdx, x_start, dE, t_start, z_end, trackID, x_end, y_end, n_electrons, n_photons, t, dx, pdgId, y, x, long_diff, z, z_start, y_start, tran_diff, t_end, pixel_plane",
                                            formats = "i8, f8, f8, f8, f8, f8, i8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, i8")
     track_zero["dE"] = 1
 
     #a track with extremely high dEdx
     track_inf = np.zeros((1,22))
-    track_inf = np.core.records.fromarrays(track_inf.transpose(), 
+    track_inf = np.core.records.fromarrays(track_inf.transpose(),
                                            names="eventID, dEdx, x_start, dE, t_start, z_end, trackID, x_end, y_end, n_electrons, n_photons, t, dx, pdgId, y, x, long_diff, z, z_start, y_start, tran_diff, t_end, pixel_plane",
                                            formats = "i8, f8, f8, f8, f8, f8, i8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, f8, i8")
     track_inf["dE"] = 1e10
@@ -49,7 +49,7 @@ def test_birksModel(self):
 
         recomb = physics.BIRKS_Ab / (1 + physics.BIRKS_kb * dedx / (detector.E_FIELD * detector.LAR_DENSITY))
 
-        assert nelectrons == pytest.approx(recomb * de * physics.MEV2ELECTRONS)
+        assert nelectrons == pytest.approx(recomb * de / physics.W_ION)
 
     def test_boxModel(self):
 
@@ -65,7 +65,7 @@ def test_boxModel(self):
         csi = physics.BOX_BETA * dedx / (detector.E_FIELD * detector.LAR_DENSITY)
         recomb = np.log(physics.BOX_ALPHA + csi)/csi
 
-        assert nelectrons == pytest.approx(recomb * de * physics.MEV2ELECTRONS)
+        assert nelectrons == pytest.approx(recomb * de / physics.W_ION)
 
     def test_birksModel_zero(self):
 
@@ -79,7 +79,7 @@ def test_birksModel_zero(self):
 
         recomb = physics.BIRKS_Ab
 
-        assert nelectrons == pytest.approx(recomb * de * physics.MEV2ELECTRONS)
+        assert nelectrons == pytest.approx(recomb * de / physics.W_ION)
 
     def test_boxModel_zero(self):
 
@@ -93,7 +93,7 @@ def test_boxModel_zero(self):
 
         recomb = 0.0
 
-        assert nelectrons == pytest.approx(recomb * de * physics.MEV2ELECTRONS)
+        assert nelectrons == pytest.approx(recomb * de / physics.W_ION)
 
     def test_birksModel_inf(self):
 
@@ -105,7 +105,7 @@ def test_birksModel_inf(self):
         quenching.quench[BPG,TPB](track_birks_inf, physics.BIRKS)
         nelectrons = track_birks_inf["n_electrons"]
 
-        recomb = nelectrons/(de * physics.MEV2ELECTRONS)
+        recomb = nelectrons/(de / physics.W_ION)
 
         assert recomb > 0 and recomb < 1e-6
 
@@ -119,6 +119,6 @@ def test_boxModel_inf(self):
         quenching.quench[BPG,TPB](tracks_box_inf, physics.BOX)
         nelectrons = tracks_box_inf["n_electrons"]
 
-        recomb = nelectrons/(de * physics.MEV2ELECTRONS)
+        recomb = nelectrons/(de / physics.W_ION)
 
         assert recomb > 0 and recomb < 1e-6