155 fluka particle source

converter/common.py

@@ -1,5 +1,5 @@
 from pathlib import Path
-from math import log10, ceil, isclose
+from math import log10, ceil, isclose, sqrt, acos, degrees, atan2
 
 
 class Parser:
@@ -90,3 +90,22 @@ def format_float(number: float, n: int) -> float:
         return 0.
 
     return result
+
+
+def cartesian2spherical(vector: tuple[float, float, float]) -> tuple[float, float, float]:
+    """
+    Transform cartesian coordinates to spherical coordinates.
+
+    :param vector: cartesian coordinates
+    :return: spherical coordinates
+    """
+    x, y, z = vector
+    r = sqrt(x**2 + y**2 + z**2)
+    # acos returns the angle in radians between 0 and pi
+    theta = degrees(acos(z / r))
+    # atan2 returns the angle in radians between -pi and pi
+    phi = degrees(atan2(y, x))
+    # lets ensure the angle in degrees is always between 0 and 360, as SHIELD-HIT12A requires
+    if phi < 0.:
+        phi += 360.
+    return theta, phi, r

converter/fluka/cards/beam_card.py

@@ -0,0 +1,72 @@
+from dataclasses import dataclass, field
+from converter.common import format_float
+from converter.fluka.cards.card import Card
+from converter.fluka.helper_parsers.beam_parser import BeamShape, FlukaBeam
+
+
+@dataclass
+class BeamCard:
+    """Class representing description of beam in Fluka input"""
+
+    data: FlukaBeam = field(default_factory=lambda: FlukaBeam())  # skipcq: PYL-W0108
+
+    def __str__(self) -> str:
+        """Return the card as a string."""
+        beam_card = Card(tag="BEAM")
+        if self.data.shape == BeamShape.GAUSSIAN:
+            shape_what = -1
+            x_y_multiplier = -1
+        elif self.data.shape == BeamShape.CIRCULAR:
+            shape_what = -1
+            x_y_multiplier = 1
+        else:
+            shape_what = 0
+            x_y_multiplier = 1
+        energy = format_float(self.data.energy*-1, 10)
+        shape_x = format_float(self.data.shape_x*x_y_multiplier, 10)
+        shape_y = format_float(self.data.shape_y*x_y_multiplier, 10)
+        if self.data.shape == BeamShape.CIRCULAR:
+            # swap x and y if beam is circular
+            # as circular beam is defined maximum and minimum radius in that order
+            # and in radius is provided in y
+            shape_x, shape_y = shape_y, shape_x
+
+        beam_card.what = [energy, 0, 0,
+                          shape_x, shape_y, shape_what]
+        beam_card.sdum = self.data.particle_name
+        if self.data.particle_name == "HEAVYION":
+            hi_card = Card(tag="HI-PROPE")
+            hi_card.what = [self.data.heavy_ion_a, self.data.heavy_ion_z, 0, 0, 0, 0]
+
+        beamposition_card = Card(tag="BEAMPOS")
+        if self.data.z_negative:
+            z_sdum = "NEGATIVE"
+        else:
+            z_sdum = ""
+
+        pos_x = format_float(self.data.beam_pos[0], 10)
+        pos_y = format_float(self.data.beam_pos[1], 10)
+        pos_z = format_float(self.data.beam_pos[2], 10)
+        dir_x = format_float(self.data.beam_dir[0], 10)
+        dir_y = format_float(self.data.beam_dir[1], 10)
+        beamposition_card.what = [pos_x, pos_y, pos_z,
+                                  dir_x, dir_y, 0]
+        beamposition_card.sdum = z_sdum
+
+        result = f"* {self.data.particle_name} beam of energy {energy*-1} GeV\n"
+        if self.data.shape == BeamShape.CIRCULAR:
+            result += f"* {self.data.shape} shape with max radius={shape_x} cm, min radius={shape_y} cm\n"
+        else:
+            result += f"* {self.data.shape} shape with x={shape_x} cm, y={shape_y} cm\n"
+        result += beam_card.__str__() + "\n"
+        if self.data.particle_name == "HEAVYION":
+            result += "* heavy ion properties: a={self.data.heavy_ion_a}, z={self.data.heavy_ion_z}\n"
+            result += hi_card.__str__() + "\n"
+        result += (f"* beam position: ({pos_x}, {pos_y}, {pos_z}) cm\n"
+                   f"* beam direction cosines in respect to x: {dir_x}, y: {dir_y}\n")
+        if self.data.z_negative:
+            result += "* beam direction is negative in respect to z axis\n"
+        else:
+            result += "* beam direction is positive in respect to z axis\n"
+        result += beamposition_card.__str__()
+        return result

converter/fluka/cards/card.py

@@ -1,7 +1,7 @@
 from dataclasses import dataclass, field
 
 
-@dataclass
+@dataclass(frozen=False)
 class Card:
     """Class representing one line (card) in Fluka input."""
 

converter/fluka/helper_parsers/beam_parser.py

@@ -0,0 +1,124 @@
+from math import cos, atan, pi
+from dataclasses import dataclass
+from enum import Enum
+
+
+class BeamShape(Enum):
+    """Enum representing beam shape"""
+
+    GAUSSIAN = 1
+    SQUARE = 2
+    CIRCULAR = 3
+
+    def __str__(self):
+        if self == BeamShape.GAUSSIAN:
+            return "gaussian"
+        if self == BeamShape.SQUARE:
+            return "flat square"
+        if self == BeamShape.CIRCULAR:
+            return "flat circular"
+        return ""
+
+
+@dataclass(frozen=False)
+class FlukaBeam:
+    """Class representing beam config in a FLUKA input file."""
+
+    energy: float = 150.  # [GeV]
+    beam_pos: tuple[float, float, float] = (0, 0, 0)  # [cm]
+    beam_dir: tuple[float, float] = (0, 0)  # cosines respective to x and y axes
+    z_negative: bool = False
+    shape: BeamShape = BeamShape.GAUSSIAN
+    shape_x: float = 0
+    shape_y: float = 0
+    particle_name: str = "PROTON"
+    heavy_ion_a: int = 1
+    heavy_ion_z: int = 1
+
+
+particle_dict = {
+    1: {"name": "NEUTRON", "a": 1},
+    2: {"name": "PROTON", "a": 1},
+    3: {"name": "PION-", "a": 1},
+    4: {"name": "PION+", "a": 1},
+    5: {"name": "PIZERO", "a": 1},
+    6: {"name": "ANEUTRON", "a": 1},
+    7: {"name": "APROTON", "a": 1},
+    8: {"name": "KAON-", "a": 1},
+    9: {"name": "KAON+", "a": 1},
+    10: {"name": "KAONZERO", "a": 1},
+    11: {"name": "KAONLONG", "a": 1},
+    12: {"name": "PHOTON", "a": 1},
+    15: {"name": "MUON-", "a": 1},
+    16: {"name": "MUON+", "a": 1},
+    21: {"name": "DEUTERON", "a": 2},
+    22: {"name": "TRITON", "a": 3},
+    23: {"name": "3-HELIUM", "a": 3},
+    24: {"name": "4-HELIUM", "a": 4},
+    25: {"name": "HEAVYION", "a": 1}
+}
+
+
+def convert_energy_to_gev(beam_json: dict) -> float:
+    """Convert energy from MeV/nucl to GeV."""
+    energy = beam_json["energy"] / 1000  # convert to GeV
+    particle = particle_dict[beam_json["particle"]["id"]]
+    if particle["name"] == "HEAVYION":
+        return energy * beam_json["particle"]["a"]
+    return energy * particle["a"]
+
+
+def parse_particle_name(particle_json: dict):
+    """Parse particle ID to FLUKA particle name."""
+    particle_id = particle_json["id"]
+    if particle_id in particle_dict:
+        particle = particle_dict[particle_id]
+        return particle["name"]
+    raise ValueError("Particle ID not supported by FLUKA")
+
+
+def parse_shape_params(shape_params_json: dict) -> tuple[BeamShape, float, float]:
+    """Parse shape params from JSON to FLUKA shape params."""
+    shape = shape_params_json["type"]
+    if shape == "Flat circular":
+        return BeamShape.CIRCULAR, shape_params_json["x"], shape_params_json["y"]
+    if shape == "Flat square":
+        return BeamShape.SQUARE, shape_params_json["x"], shape_params_json["y"]
+    if shape == "Gaussian":
+        return BeamShape.GAUSSIAN, shape_params_json["x"], shape_params_json["y"]
+    raise ValueError("Shape type not supported by FLUKA")
+
+
+def cartesian_to_spherical(coords: tuple[float, float, float]):
+    """
+    Convert cartesian coordinates to spherical coordinates
+    and return cosines of angles respective to x and y axes
+    """
+    x, y, z = coords
+    theta = pi/2
+    if x != 0:
+        theta = atan(z/x)
+    phi = pi/2
+    if y != 0:
+        phi = atan((x**2 + z**2)**0.5 / y)
+    return cos(theta), cos(phi)
+
+
+def parse_beam(beam_json: dict) -> FlukaBeam:
+    """Parse beam from JSON to FLUKA beam."""
+    fluka_beam = FlukaBeam()
+    fluka_beam.energy = convert_energy_to_gev(beam_json)
+    fluka_beam.particle_name = parse_particle_name(beam_json["particle"])
+    if fluka_beam.particle_name == "HEAVYION":
+        fluka_beam.heavy_ion_a = beam_json["particle"]["a"]
+        fluka_beam.heavy_ion_z = beam_json["particle"]["z"]
+    fluka_beam.beam_pos = tuple(beam_json["position"])
+    shape, shape_x, shape_y = parse_shape_params(beam_json["sigma"])
+    fluka_beam.shape = shape
+    fluka_beam.shape_x = shape_x
+    fluka_beam.shape_y = shape_y
+    theta, phi = cartesian_to_spherical(beam_json["direction"])
+    fluka_beam.beam_dir = (theta, phi)
+    if beam_json["direction"][2] < 0:
+        fluka_beam.z_negative = True
+    return fluka_beam

converter/fluka/input.py

@@ -1,15 +1,17 @@
 from dataclasses import dataclass, field
+from converter.fluka.cards.beam_card import BeamCard
 from converter.fluka.cards.card import Card
 from converter.fluka.cards.figure_card import FiguresCard
 from converter.fluka.cards.region_card import RegionsCard
+from converter.fluka.helper_parsers.beam_parser import FlukaBeam
 from converter.solid_figures import SolidFigure
 
 
 @dataclass
 class Input:
     """Class mapping of the Fluka input file."""
 
-    energy_GeV: float = 0.07  # GeV FLUKA specific
+    beam: FlukaBeam = field(default_factory=lambda: FlukaBeam())  # skipcq: PYL-W0108
     number_of_particles: int = 10000
 
     figures: list[SolidFigure] = field(default_factory=lambda: [])
@@ -19,10 +21,7 @@ class Input:
 proton beam simulation
 * default physics settings for hadron therapy
 DEFAULTS                                                              HADROTHE
-* beam source
 {BEAM}
-* beam source position
-BEAMPOS          0.0       0.0    -100.0
 * geometry description starts here
 GEOBEGIN                                                              COMBNAME
     0    0
@@ -68,8 +67,8 @@ class Input:
     def __str__(self):
         """Return fluka input file as string"""
         return self.template.format(
-            BEAM=Card(tag="BEAM", what=[str(-self.energy_GeV)], sdum="PROTON"),
             START=Card(tag="START", what=[str(self.number_of_particles)]),
+            BEAM=BeamCard(data=self.beam),
             FIGURES=FiguresCard(data=self.figures),
             REGIONS=RegionsCard(data=self.regions)
         )

converter/fluka/parser.py

@@ -1,4 +1,5 @@
 from converter.common import Parser
+from converter.fluka.helper_parsers.beam_parser import parse_beam
 from converter.fluka.helper_parsers.figure_parser import parse_figures
 from converter.fluka.helper_parsers.region_parser import parse_regions
 from converter.fluka.input import Input
@@ -18,13 +19,12 @@ def __init__(self) -> None:
 
     def parse_configs(self, json: dict) -> None:
         """Parse energy and number of particles from json."""
-        # Since energy in json is in MeV and FLUKA uses GeV, we need to convert it.
-        self.input.energy_GeV = float(json["beam"]["energy"]) * 1e-3
         self.input.number_of_particles = json["beam"]["numberOfParticles"]
 
         self.input.figures = parse_figures(json["figureManager"].get('figures'))
         self.input.regions, world_figures = parse_regions(json["zoneManager"], self.input.figures)
         self.input.figures.extend(world_figures)
+        self.input.beam = parse_beam(json["beam"])
 
     def get_configs_json(self) -> dict:
         """

tests/fluka/conftest.py

@@ -9,6 +9,6 @@ def project_fluka_path() -> Path:
 
 @pytest.fixture(scope='session')
 def project_fluka_json(project_fluka_path) -> dict:
-    """Dictionary with project data for SHIELD-HIT12A"""
+    """Dictionary with project data for Fluka"""
     with open(project_fluka_path, 'r') as file_handle:
         return json.load(file_handle)

tests/fluka/expected_fluka_output/fl_sim.py

@@ -2,10 +2,13 @@
 proton beam simulation
 * default physics settings for hadron therapy
 DEFAULTS                                                              HADROTHE
-* beam source
-BEAM           -0.07                                                  PROTON    
-* beam source position
-BEAMPOS          0.0       0.0    -100.0
+* PROTON beam of energy 0.07 GeV
+* flat circular shape with max radius=3.0 cm, min radius=0.0 cm
+BEAM           -0.07       0.0       0.0       3.0       0.0      -1.0PROTON    
+* beam position: (0.0, 0.0, -1.5) cm
+* beam direction cosines in respect to x: 0.0, y: 0.0
+* beam direction is positive in respect to z axis
+BEAMPOS          0.0       0.0      -1.5       0.0       0.0       0.0          
 * geometry description starts here
 GEOBEGIN                                                              COMBNAME
     0    0

tests/fluka/test_fluka_beam.py

@@ -0,0 +1,35 @@
+import copy
+from math import isclose
+from converter.fluka.helper_parsers.beam_parser import BeamShape, parse_beam
+
+
+def test_parse_fluka_beam(project_fluka_json):
+    """Test if Fluka beam is parsed correctly"""
+    beam_json = project_fluka_json["beam"]
+
+    fluka_beam = parse_beam(beam_json)
+
+    assert fluka_beam.energy == 0.07
+    assert fluka_beam.particle_name == "PROTON"
+    assert fluka_beam.shape == BeamShape.CIRCULAR
+    assert fluka_beam.shape_x == 0
+    assert fluka_beam.shape_y == 3
+    assert fluka_beam.z_negative is False
+    assert fluka_beam.beam_pos == (0, 0, -1.5)
+    assert isclose(fluka_beam.beam_dir[0], 0.0, abs_tol=1e-16)
+    assert isclose(fluka_beam.beam_dir[1], 0.0, abs_tol=1e-16)
+
+
+def test_parse_heavy_ions(project_fluka_json):
+    """Test if Fluka beam is parsed correctly"""
+    beam_json = copy.deepcopy(project_fluka_json["beam"])
+    beam_json["particle"]["id"] = 25
+    beam_json["particle"]["a"] = 2
+    beam_json["particle"]["z"] = 3
+
+    fluka_beam = parse_beam(beam_json)
+
+    assert fluka_beam.energy == 0.14
+    assert fluka_beam.particle_name == "HEAVYION"
+    assert fluka_beam.heavy_ion_a == 2
+    assert fluka_beam.heavy_ion_z == 3