From 1429953ba4c37d41df541064f27e38f6b7fcd2e3 Mon Sep 17 00:00:00 2001
From: Johannes Schumann <johannes.schumann@fau.de>
Date: Thu, 25 Nov 2021 13:22:33 +0000
Subject: [PATCH] Proposal 7
---
km3buu/config.py | 3 +-
km3buu/propagation.py | 78 +++++++++++++++++++-------------
km3buu/tests/test_propagation.py | 48 +++++++++++++-------
requirements.txt | 2 +-
4 files changed, 80 insertions(+), 51 deletions(-)
diff --git a/km3buu/config.py b/km3buu/config.py
index 25279ce..b9c7a0c 100644
--- a/km3buu/config.py
+++ b/km3buu/config.py
@@ -15,6 +15,7 @@ from . import IMAGE_NAME
from .environment import build_image
import mendeleev
import xml.etree.ElementTree as ElementTree
+import tempfile
__author__ = "Johannes Schumann"
__copyright__ = "Copyright 2020, Johannes Schumann and the KM3NeT collaboration."
@@ -101,7 +102,7 @@ class Config(object):
@property
def proposal_itp_tables(self):
- default_path = abspath(join(dirname(__file__), "../.tables"))
+ default_path = tempfile.gettempdir()
return self.get(PROPOSAL_SECTION, PROPOSAL_ITP_PATH_KEY, default_path)
@proposal_itp_tables.setter
diff --git a/km3buu/propagation.py b/km3buu/propagation.py
index 6a0564c..61b59a9 100644
--- a/km3buu/propagation.py
+++ b/km3buu/propagation.py
@@ -18,6 +18,7 @@ import proposal as pp
from particle import Particle
import awkward as ak
from collections import defaultdict
+import pathlib
from .config import Config
@@ -36,23 +37,34 @@ PROPOSAL_LEPTON_DEFINITIONS = {
-16: pp.particle.NuTauBarDef
}
+PROPOSAL_TARGET = pp.medium.AntaresWater()
-def _setup_propagator(max_distance, lepton_definition):
- sector = pp.SectorDefinition()
- sector.medium = pp.medium.AntaresWater(1.0)
- geometry = pp.geometry.Sphere(pp.Vector3D(), max_distance, 0)
- sector.geometry = geometry
- sector.scattering_model = pp.scattering.ScatteringModel.Highland
- cfg = Config()
- interpolation = pp.InterpolationDef()
- interpolation.path_to_tables = cfg.proposal_itp_tables
- interpolation.path_to_tables_readonly = cfg.proposal_itp_tables
+def _setup_propagator(max_distance, particle_definition):
+ itp_table_path = Config().proposal_itp_tables
+ pathlib.Path(itp_table_path).mkdir(parents=True, exist_ok=True)
+ pp.InterpolationSettings.tables_path = itp_table_path
+ crosssection = pp.crosssection.make_std_crosssection(
+ particle_def=particle_definition,
+ target=PROPOSAL_TARGET,
+ interpolate=True,
+ cuts=pp.EnergyCutSettings(np.inf, 0.05, False))
+ collection = pp.PropagationUtilityCollection()
+ collection.displacement = pp.make_displacement(crosssection, True)
+ collection.interaction = pp.make_interaction(crosssection, True)
+ collection.decay = pp.make_decay(crosssection, particle_definition, True)
+ collection.time = pp.make_time(crosssection, particle_definition, True)
- return pp.Propagator(sector_defs=[sector],
- particle_def=lepton_definition,
- detector=geometry,
- interpolation_def=interpolation)
+ utility = pp.PropagationUtility(collection=collection)
+
+ geometry = pp.geometry.Sphere(pp.Cartesian3D(), max_distance)
+ density_distr = pp.density_distribution.density_homogeneous(
+ PROPOSAL_TARGET.mass_density)
+
+ propagator = pp.Propagator(particle_definition,
+ [(geometry, utility, density_distr)])
+
+ return propagator
def propagate_lepton(lepout_data, pdgid):
@@ -76,28 +88,30 @@ def propagate_lepton(lepout_data, pdgid):
lepout_data.lepOut_E) / lepton_info.mass #[cm]
lepton_def = PROPOSAL_LEPTON_DEFINITIONS[pdgid]()
- lepton = pp.particle.DynamicData(lepton_def.particle_type)
- propagated_data = defaultdict(list)
propagator = _setup_propagator(10 * prop_range, lepton_def)
+ propagated_data = defaultdict(list)
+
for entry in lepout_data:
- lepton.energy = entry.lepOut_E * 1e3
- lepton.position = pp.Vector3D(0, 0, 0)
- lepton.direction = pp.Vector3D(entry.lepOut_Px, entry.lepOut_Py,
- entry.lepOut_Pz)
- lepton.direction.normalize()
- secondaries = propagator.propagate(lepton, 5 * prop_range)
-
- E = np.array(secondaries.energy) / 1e3
- pdgid = [p.type for p in secondaries.particles]
- Px = [p.direction.x * p.momentum / 1e3 for p in secondaries.particles]
- Py = [p.direction.y * p.momentum / 1e3 for p in secondaries.particles]
- Pz = [p.direction.z * p.momentum / 1e3 for p in secondaries.particles]
- x = [p.position.x / 100 for p in secondaries.particles]
- y = [p.position.y / 100 for p in secondaries.particles]
- z = [p.position.z / 100 for p in secondaries.particles]
- dt = [p.time for p in secondaries.particles]
+ init_state = pp.particle.ParticleState()
+ init_state.energy = entry.lepOut_E * 1e3
+ init_state.position = pp.Cartesian3D(0, 0, 0)
+ init_state.direction = pp.Cartesian3D(entry.lepOut_Px, entry.lepOut_Py,
+ entry.lepOut_Pz)
+ init_state.direction.normalize()
+ track = propagator.propagate(init_state, 5 * prop_range)
+ secondaries = track.decay_products()
+
+ E = [s.energy / 1e3 for s in secondaries]
+ pdgid = [p.type for p in secondaries]
+ Px = [p.direction.x * p.momentum / 1e3 for p in secondaries]
+ Py = [p.direction.y * p.momentum / 1e3 for p in secondaries]
+ Pz = [p.direction.z * p.momentum / 1e3 for p in secondaries]
+ x = [p.position.x / 100 for p in secondaries]
+ y = [p.position.y / 100 for p in secondaries]
+ z = [p.position.z / 100 for p in secondaries]
+ dt = [p.time for p in secondaries]
propagated_data["E"].append(E)
propagated_data["Px"].append(Px)
diff --git a/km3buu/tests/test_propagation.py b/km3buu/tests/test_propagation.py
index 6b013f0..27d1704 100644
--- a/km3buu/tests/test_propagation.py
+++ b/km3buu/tests/test_propagation.py
@@ -14,42 +14,53 @@ import csv
import unittest
import numpy as np
import uproot
+import pytest
from os.path import abspath, join, dirname
from thepipe.logger import get_logger
-from km3net_testdata import data_path
import proposal as pp
-from km3buu.output import GiBUUOutput
+import awkward as ak
from km3buu.propagation import propagate_lepton
-TESTDATA_DIR = data_path("gibuu")
-
pp.RandomGenerator.get().set_seed(1234)
+@pytest.mark.skip(reason="CI boost lib problem")
class TestTauPropagation(unittest.TestCase):
def setUp(self):
- log = get_logger("ctrl.py")
- log.setLevel("INFO")
- self.gibuu_output = GiBUUOutput(TESTDATA_DIR)
- fname = join(TESTDATA_DIR, self.gibuu_output.root_pert_files[0])
- fobj = uproot.open(fname)
- data = fobj["RootTuple"].arrays()
+ data = ak.Array({
+ "lepOut_E": [
+ 17.45830624434573, 3.1180989952362594, 21.270059768902005,
+ 5.262659790136034, 23.52185741888274
+ ],
+ "lepOut_Px": [
+ -0.42224402086330426, -1.0232258668453014, -0.5801431899058521,
+ -0.9038349288874724, 0.9022573877437422
+ ],
+ "lepOut_Py": [
+ 0.3644190693190108, -0.24542303987320932, 0.24499631087268617,
+ -1.1060562370375715, -3.982173292871768
+ ],
+ "lepOut_Pz": [
+ 17.35867612031871, 2.336148261778657, 21.186342871282157,
+ 4.743161507744377, 23.096499191566885
+ ]
+ })
self.sec = propagate_lepton(data, 15)
def test_secondary_momenta(self):
np.testing.assert_array_almost_equal(np.array(self.sec[0].E),
- [0.535, 1.316, 0.331],
+ [2.182, 13.348, 1.928],
decimal=3)
np.testing.assert_array_almost_equal(np.array(self.sec[0].Px),
- [-0.467, 0.321, -0.246],
+ [0.295, -0.48, -0.237],
decimal=3)
np.testing.assert_array_almost_equal(np.array(self.sec[0].Py),
- [0.127, -0.822, 0.218],
+ [-0.375, 0.784, -0.044],
decimal=3)
np.testing.assert_array_almost_equal(np.array(self.sec[0].Pz),
- [0.179, 0.967, -0.041],
+ [2.129, 13.316, 1.913],
decimal=3)
def test_secondary_types(self):
@@ -57,9 +68,12 @@ class TestTauPropagation(unittest.TestCase):
[13, 16, -14])
def test_secondary_positions(self):
- np.testing.assert_array_almost_equal(np.array(self.sec[0].x), [0, 0],
+ np.testing.assert_array_almost_equal(np.array(self.sec[0].x),
+ [-1.4e-05, -1.4e-05, -1.4e-05],
decimal=1)
- np.testing.assert_array_almost_equal(np.array(self.sec[0].y), [0, 0],
+ np.testing.assert_array_almost_equal(np.array(self.sec[0].y),
+ [1.2e-05, 1.2e-05, 1.2e-05],
decimal=1)
- np.testing.assert_array_almost_equal(np.array(self.sec[0].z), [0, 0],
+ np.testing.assert_array_almost_equal(np.array(self.sec[0].z),
+ [0., 0., 0.],
decimal=1)
diff --git a/requirements.txt b/requirements.txt
index d9a0795..73cffa8 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,4 +10,4 @@ uproot>=4.0.0
awkward>=1.0.0
pandas
mendeleev
-proposal==6.1.6
+proposal>=7.1.1
--
GitLab