diff --git a/converter/fluka/cards/beam_card.py b/converter/fluka/cards/beam_card.py new file mode 100644 index 00000000..509892b2 --- /dev/null +++ b/converter/fluka/cards/beam_card.py @@ -0,0 +1,79 @@ +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 + # in Fluka input file positive number is interpreted as momentum + # and negative as energy + # we store energy in FlukaBeam object as positive number, + # so we need to multiply it by -1 + # we also divide it by 1000 to convert from MeV to GeV + momentum_or_energy = format_float(-self.data.energy_MeV / 1000, 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 = [momentum_or_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") + + # z_negative is True if beam direction is negative in respect to z axis + 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 {-momentum_or_energy} 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 += f"* 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 diff --git a/converter/fluka/cards/card.py b/converter/fluka/cards/card.py index bd153b04..3f101277 100644 --- a/converter/fluka/cards/card.py +++ b/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.""" diff --git a/converter/fluka/helper_parsers/beam_parser.py b/converter/fluka/helper_parsers/beam_parser.py new file mode 100644 index 00000000..21b91d2c --- /dev/null +++ b/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_MeV: float = 150. + 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(beam_json: dict) -> float: + """Convert energy from MeV/nucl to MeV.""" + energy = beam_json["energy"] + 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_MeV = convert_energy(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 diff --git a/converter/fluka/input.py b/converter/fluka/input.py index 1082a85e..6ec9a9a3 100644 --- a/converter/fluka/input.py +++ b/converter/fluka/input.py @@ -1,7 +1,9 @@ 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 @@ -9,7 +11,7 @@ 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) ) diff --git a/converter/fluka/parser.py b/converter/fluka/parser.py index 9caa20f3..9838e03f 100644 --- a/converter/fluka/parser.py +++ b/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: """ diff --git a/tests/fluka/conftest.py b/tests/fluka/conftest.py index 22b8c5cf..a3768817 100644 --- a/tests/fluka/conftest.py +++ b/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) \ No newline at end of file diff --git a/tests/fluka/expected_fluka_output/fl_sim.py b/tests/fluka/expected_fluka_output/fl_sim.py index 385e95c7..8eb72044 100644 --- a/tests/fluka/expected_fluka_output/fl_sim.py +++ b/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 diff --git a/tests/fluka/test_fluka_beam.py b/tests/fluka/test_fluka_beam.py new file mode 100644 index 00000000..b93b87ce --- /dev/null +++ b/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_MeV == 70 + 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"] = 6 + beam_json["particle"]["z"] = 12 + + fluka_beam = parse_beam(beam_json) + + assert fluka_beam.energy_MeV == 420 + assert fluka_beam.particle_name == "HEAVYION" + assert fluka_beam.heavy_ion_a == 6 + assert fluka_beam.heavy_ion_z == 12 \ No newline at end of file