diff --git a/opengate/SimulationEngine.py b/opengate/SimulationEngine.py index 1170a0ecf..7e3878e7b 100644 --- a/opengate/SimulationEngine.py +++ b/opengate/SimulationEngine.py @@ -4,7 +4,14 @@ import sys import os from .ExceptionHandler import * -from multiprocessing import Process, set_start_method, Queue +from multiprocessing import ( + Process, + set_start_method, + Manager, + Queue, + active_children, + cpu_count, +) from opengate_core import G4RunManagerFactory from .Decorators import requires_fatal from .helpers import fatal, warning @@ -131,12 +138,28 @@ def start(self): # (the "force" option is needed for notebooks) set_start_method("fork", force=True) # set_start_method("spawn") - q = Queue() + q = Manager().Queue() + # q = Queue() p = Process(target=self.init_and_start, args=(q,)) + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + print(f"Queue full: {q.full()}") + print("---start process---") p.start() + import time + + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + print(f"Queue full: {q.full()}") + while len(active_children()) >= cpu_count() + 4: + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + time.sleep(0.01) + print(q.full()) self.state = "started" - p.join() + p.join() # (timeout=10) # timeout might be needed self.state = "after" + print("AFTER") output = q.get() else: output = self.init_and_start(None) @@ -184,7 +207,15 @@ def init_and_start(self, queue): output.store_hook_log(self) output.current_random_seed = self.current_random_seed if queue is not None: + print("--- in process, before put ---") + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + print(f"Queue full: {queue.full()}") queue.put(output) + print("--- in process, after put ---") + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + print(f"Queue full: {queue.full()}") return None else: return output diff --git a/opengate/__init__.py b/opengate/__init__.py index 6bdb1e077..78efc3abc 100644 --- a/opengate/__init__.py +++ b/opengate/__init__.py @@ -10,6 +10,8 @@ from .helpers_tests import * from .helpers_tests_root import * from .helpers_transform import * +from .helpers_beamline import * +from .helpers_rt_plan import * # main mechanism for the 'elements': source, actor, volume from .UserInfo import * @@ -17,6 +19,7 @@ from .source.SourceBase import * from .actor.ActorBase import * from .geometry.VolumeBase import * +from .source.TreatmentPlanSource import * # main object from .Simulation import * diff --git a/opengate/helpers_beamline.py b/opengate/helpers_beamline.py new file mode 100644 index 000000000..02995878a --- /dev/null +++ b/opengate/helpers_beamline.py @@ -0,0 +1,69 @@ +# N.B: distances in mm, degrees in rad + + +class Rashi: + pass + + +class RangeMod: + pass + + +class BeamlineModel: + def __init__(self): + self.name = None + self.radiation_types = [] + self.rm = None # range modulator + self.rashi = None + # Nozzle entrance to Isocenter distance + self.distance_nozzle_iso = 0 # mm + # SMX (X bending magnet) to Isocenter distance + self.distance_stearmag_to_isocenter_x = float( + "inf" + ) # default infinity for parallel beams + # SMY (Y bending magnet) to Isocenter distance + self.distance_stearmag_to_isocenter_y = float("inf") + # polinomial coefficients + self.energy_mean_coeffs = [0] + self.energy_spread_coeffs = [0] + self.sigma_x_coeffs = [0] + self.theta_x_coeffs = [0] + self.epsilon_x_coeffs = [0] + self.sigma_y_coeffs = [0] + self.theta_y_coeffs = [0] + self.epsilon_y_coeffs = [0] + # convergence + self.conv_x = 0 + self.conv_y = 0 + + def _polynomial_map(self, base, coeff): + # coeff are given with decreasing degree (coeff[0]->max degree) + polyDegree = len(coeff) + exp = list(range(polyDegree)) + exp.reverse() + + return sum([c * (base ** (i)) for c, i in zip(coeff, exp)]) + + def get_energy(self, nominal_energy): + return self._polynomial_map(nominal_energy, self.energy_mean_coeffs) + + def get_sigma_energy(self, nominal_energy): + return self._polynomial_map(nominal_energy, self.energy_spread_coeffs) + + def get_sigma_x(self, energy): + return self._polynomial_map(energy, self.sigma_x_coeffs) + + def get_theta_x(self, energy): + return self._polynomial_map(energy, self.theta_x_coeffs) + + def get_epsilon_x(self, energy): + return self._polynomial_map(energy, self.epsilon_x_coeffs) + + def get_sigma_y(self, energy): + return self._polynomial_map(energy, self.sigma_y_coeffs) + + def get_theta_y(self, energy): + return self._polynomial_map(energy, self.theta_y_coeffs) + + def get_epsilon_y(self, energy): + return self._polynomial_map(energy, self.epsilon_y_coeffs) diff --git a/opengate/helpers_rt_plan.py b/opengate/helpers_rt_plan.py new file mode 100644 index 000000000..0d8ba9492 --- /dev/null +++ b/opengate/helpers_rt_plan.py @@ -0,0 +1,698 @@ +#!/usr/bin/env python3 +# ----------------------------------------------------------------------------- +# Copyright (C): MedAustron GmbH, ACMIT Gmbh and Medical University Vienna +# This software is distributed under the terms +# of the GNU Lesser General Public Licence (LGPL) +# See LICENSE for further details +# ----------------------------------------------------------------------------- + +import pydicom +import os +import re +import numpy as np + +# from utils.dose_info import dose_info +import logging + +logger = logging.getLogger(__name__) + + +def is_close(x, y, eps=1e-6): + sumabs = np.abs(x) + np.abs(y) + absdif = np.abs(x - y) + ok = (sumabs == 0) or (absdif < eps * 0.5 * sumabs) + return ok + + +class spot_info(object): + def __init__(self, xiec, yiec, w, e): + self.xiec = xiec + self.yiec = yiec + self.w = w + self.energy = e + self.particle_name = None + self.beamFraction = None + self.t0 = None + self.t1 = None + + def get_msw(self, t0, t1): + return self.w + + +class layer_info(object): + def __init__(self, ctrlpnt, j, cumsumchk=[], verbose=False, keep0=False): + self._cp = ctrlpnt + if verbose: + logger.debug("{}. control point with type {}".format(j, type(self._cp))) + for k in self._cp.keys(): + if pydicom.datadict.dictionary_has_tag(k): + kw = pydicom.datadict.keyword_for_tag(k) + else: + kw = "(UNKNOWN)" + logger.debug("k={} keyword={}".format(k, kw)) + nspot = int(self._cp.NumberOfScanSpotPositions) + # assert(self._cp.NominalBeamEnergyUnit == 'MEV') + if nspot == 1: + self.w = np.array([float(self._cp.ScanSpotMetersetWeights)]) + else: + self.w = np.array([float(w) for w in self._cp.ScanSpotMetersetWeights]) + assert nspot == len(self.w) + assert nspot * 2 == len(self._cp.ScanSpotPositionMap) + # self.cpindex = int(self._cp.ControlPointIndex) + # self.spotID = str(self._cp.ScanSpotTuneID) + cmsw = float(self._cp.CumulativeMetersetWeight) + if cumsumchk: + logger.debug( + "CumulativeMetersetWeight={0:14.8g} sum of previous spots={1:14.8g} diff={2:14.8g}".format( + cmsw, cumsumchk[0], cmsw - cumsumchk[0] + ) + ) + assert is_close(cmsw, cumsumchk[0]) + # self.npainting = int(self._cp.NumberOfPaintings) + xy = np.array([float(pos) for pos in self._cp.ScanSpotPositionMap]).reshape( + nspot, 2 + ) + self.x = np.array(xy[:, 0]) + self.y = np.array(xy[:, 1]) + if not keep0: + mask = self.w > 0.0 + self.w = self.w[mask] + self.x = self.x[mask] + self.y = self.y[mask] + + wsum = np.sum(self.w) + logger.debug( + "layer number {} has {} spots, of which {} with zero weight, cumsum={}, sum(w)={}".format( + j, len(self.w), np.sum(self.w <= 0), cmsw, wsum + ) + ) + + cumsumchk[0] += wsum + + @property + def energy(self): + # DICOM specifies energy per nucleon, Gate wants total kinetic energy + return float(self._cp.NominalBeamEnergy) + + @property + def tuneID(self): + return str(self._cp.ScanSpotTuneID) + + @property + def npainting(self): + return int(self._cp.NumberOfPaintings) + + @property + def mswtot(self): + return np.sum(self.w) + + @property + def nspots(self): + return len(self.w) + + @property + def weights(self): + return self.w + + @property + def spots(self): + e = self.energy + return [spot_info(x, y, w, e) for (x, y, w) in zip(self.x, self.y, self.w)] + + def get_spots(self, t0=None, t1=None): + e = self.energy + return [spot_info(x, y, w, e) for (x, y, w) in zip(self.x, self.y, self.w)] + + +class beam_info(object): + # def __init__(self,beam,rd,i,keep0=False): + def __init__(self, beam, i, override_number, keep0=False): + logger.debug("loading {}th beam".format(i)) + self._dcmbeam = beam # the DICOM beam object + self._warnings = list() # will hopefully stay empty + logger.debug("trying to access first control point") + self._icp0 = beam.IonControlPointSequence[0] # convenience: first control point + self._beam_number_is_fishy = ( + override_number # workaround for buggy TPSs, e.g. PDM + ) + self._index = i # the index in the beam sequence + self._layers = list() + mswchk = self.FinalCumulativeMetersetWeight + cumsumchk = [0.0] + logger.debug("going to read all layers") + for j, icp in enumerate(self._dcmbeam.IonControlPointSequence): + li = layer_info(icp, j, cumsumchk, False, keep0) + if 0.0 < li.mswtot or keep0: + self._layers.append(li) + logger.debug("survived reading all layers") + if not is_close(mswchk, cumsumchk[0]): + raise ValueError( + "final cumulative msw {} != sum of spot msw {}".format( + mswchk, cumsumchk[0] + ) + ) + logger.debug("survived cumulative MSW check") + + def GetAndClearWarnings(self): + # return and clear + w = self._warnings[:] + self._warnings = list() + return w + + @property + def FinalCumulativeMetersetWeight(self): + return float(self._dcmbeam.FinalCumulativeMetersetWeight) + + @property + def PatientSupportAngle(self): + return float(self._icp0.PatientSupportAngle) + + @property + def patient_angle(self): + return float(self._icp0.PatientSupportAngle) + + @property + def IsoCenter(self): + if "IsocenterPosition" in self._icp0: + if len(self._icp0.IsocenterPosition) == 3.0: + return [float(xyz) for xyz in self._icp0.IsocenterPosition] + else: + msg = "Got corrupted isocenter = '{}'; assuming [0,0,0] for now, keep fingers crossed.".format( + self._icp0.IsocenterPosition + ) + else: + msg = "No isocenter specified in treatment plan; assuming [0,0,0] for now, keep fingers crossed." + logger.error(msg) + self._warnings.append(msg) + # FIXME: what to do else? Cry? Throw segfaults? Drink bad coffee? + return [0.0, 0.0, 0.0] + + @property + def Name(self): + # TODO: the Name and name properties are identical, keep only one of them. + return str(self._dcmbeam.BeamName) + + @property + def Number(self): + # TODO: the Number and number properties are identical, keep only one of them. + nr = ( + str(self._index + 1) + if self._beam_number_is_fishy + else str(self._dcmbeam.BeamNumber) + ) + return nr + + @property + def name(self): + # TODO: the Name and name properties are identical, keep only one of them. + return str(self._dcmbeam.BeamName) + + @property + def number(self): + # TODO: the Number and number properties are identical, keep only one of them. + nr = ( + str(self._index + 1) + if self._beam_number_is_fishy + else str(self._dcmbeam.BeamNumber) + ) + return nr + + @property + def RadiationType(self): + radtype = str(self._dcmbeam.RadiationType) + radtypeOpengate = None + if radtype == "ION": + ionZ = str(self._dcmbeam.RadiationAtomicNumber) + ionA = str(self._dcmbeam.RadiationMassNumber) + ionQ = str(self._dcmbeam.RadiationChargeState) + radtype = "_".join(["ION", ionZ, ionA, ionQ]) + radtypeOpengate = f"ion {ionZ} {ionA}" + return [radtype, radtypeOpengate] + + @property + def gantry_angle(self): + return float(self._icp0.GantryAngle) + + @property + def TreatmentMachineName(self): + if "TreatmentMachineName" in self._dcmbeam: + return str(self._dcmbeam.TreatmentMachineName) + # RayStation 8b exports anonymized treatment plans without the treatment machine name! + if np.isclose(self.gantry_angle, 0.0): + # FIXME: should be solved in a way that works for any clinic, not just MedAustron + ducktape = str("IR2VBL") + elif np.isclose(self.gantry_angle, 90.0): + # FIXME: should be solved in a way that works for any clinic, not just MedAustron + ducktape = str("IR2HBL") + else: + raise ValueError( + "treatment machine name is missing and gantry angle {} does not enable a good guess".format( + self.gantry_angle + ) + ) + msg = "treatment machine name for beam name={} number={} missing in treatment plan, guessing '{}' from gantry angle {}".format( + self.name, self.number, ducktape, self.gantry_angle + ) + logger.error(msg) + self._warnings.append(msg) + return ducktape # ugly workaround! FIXME! + + @property + def SnoutID(self): + if "SnoutSequence" in self._dcmbeam: + return str(self._dcmbeam.SnoutSequence[0].SnoutID) + # FIXME: what to do else? + return str("NA") + + @property + def SnoutPosition(self): + if "SnoutPosition" in self._dcmbeam: + return float(self._icp0.SnoutPosition) + # FIXME: what to do else? + return str("NA") + + @property + def NumberOfRangeModulators(self): + return int(self._dcmbeam.NumberOfRangeModulators) + + @property + def RangeModulatorIDs(self): + if self.NumberOfRangeModulators > 0: + return [rm.RangeModulatorID for rm in self._dcmbeam.RangeModulatorSequence] + return list() + + @property + def NumberOfRangeShifters(self): + return int(self._dcmbeam.NumberOfRangeShifters) + + @property + def RangeShifterIDs(self): + if self.NumberOfRangeShifters > 0: + return [str(rs.RangeShifterID) for rs in self._dcmbeam.RangeShifterSequence] + return list() + + @property + def NumberOfEnergies(self): + return len( + set( + [icp.NominalBeamEnergy for icp in self._dcmbeam.IonControlPointSequence] + ) + ) + + @property + def nlayers(self): + return len(self._layers) + + @property + def layers(self): + return self._layers + + @property + def nspots(self): + return sum([l.nspots for l in self.layers]) + + @property + def mswtot(self): + return sum([l.mswtot for l in self._layers]) + + @property + def PrimaryDosimeterUnit(self): + return str(self._dcmbeam.PrimaryDosimeterUnit) + + +def sequence_check(obj, attr, nmin=1, nmax=0, name="object"): + logger.debug("checking that {} has attribute {}".format(name, attr)) + assert hasattr(obj, attr) + seq = getattr(obj, attr) + logger.debug( + "{} has length {}, will check if it >={} and <={}".format( + name, len(seq), nmin, nmax + ) + ) + assert len(seq) >= nmin + assert nmax == 0 or len(seq) <= nmax + + +class beamset_info(object): + """ + This class reads a DICOM 'RT Ion Plan Storage' file and collects related information such as TPS dose files. + It does NOT (yet) try to read a reffered structure set and/or CT images. + This acts as a wrapper (all DICOM access on the plan file happens here). This has a few advantages over direct + DICOM access in the other modules: + * we can deal with different "DICOM dialects" here; some TPSs may store their plans in different ways. + * if 'private tags' need to be taken into account then we can also do that here. + * We can make a similar class, with the same attributes, for a treatment plan stored in a different format, e.g. for research, commissioning or QA purposes. + + Then the rest of the code can work with that in the same way. + """ + + patient_attrs = ["Patient ID", "Patient Name", "Patient Birth Date", "Patient Sex"] + plan_req_attrs = [ + "RT Plan Label", + "SOP Instance UID", + "Referring Physician Name", + "Plan Intent", + ] + plan_opt_attrs = ["Operators Name", "Reviewer Name", "Review Date", "Review Time"] + plan_attrs = plan_req_attrs + plan_opt_attrs + bs_attrs = [ + "Number Of Beams", + "RT Plan Label", + "Prescription Dose", + "Target ROI Name", + "Radiation Type", + "Treatment Machine(s)", + ] + + def __init__(self, rpfp): + self._warnings = list() # will hopefully stay empty + self._beam_numbers_corrupt = False # e.g. PDM does not define beam numbers + self._rp = pydicom.read_file(rpfp) + self._rpfp = rpfp + logger.debug("beamset: survived reading DICOM file {}".format(rpfp)) + self._rpdir = os.path.dirname(rpfp) + self._rpuid = str(self._rp.SOPInstanceUID) + self._dose_roiname = ( + None # stays None for CT-less plans, e.g. commissioning plans + ) + self._dose_roinumber = ( + None # stays None for CT-less plans, e.g. commissioning plans + ) + logger.debug("beamset: going to do some checks") + self._chkrp() + logger.debug("beamset: survived check, loading beams") + self._beams = [ + beam_info(b, i, self._beam_numbers_corrupt) + for i, b in enumerate(self._rp.IonBeamSequence) + ] + logger.debug("beamset: DONE") + + def GetAndClearWarnings(self): + # return a copy + for b in self._beams: + # bwarnings = b.GetAndClearWarnings() + for w in b.GetAndClearWarnings(): + if w not in self._warnings: + self._warnings.append(w) + allw = self._warnings[:] + self._warnings = list() + return allw + + def __getitem__(self, k): + if type(k) == int: + if k >= 0 and k < len(self._beams): + return self._beams[k] + else: + raise IndexError( + "attempt to get nonexisting beam with index {}".format(k) + ) + for b in self._beams: + if str(k) == b.name or str(k) == b.number: + return b + raise KeyError("attempt to get nonexisting beam with key {}".format(k)) + + def _chkrp(self): + if "SOPClassUID" not in self._rp: + raise IOError("bad DICOM file {},\nmissing SOPClassUID".format(self._rpfp)) + sop_class_name = pydicom.uid.UID_dictionary[self._rp.SOPClassUID][0] + if sop_class_name != "RT Ion Plan Storage": + raise IOError( + "bad plan file {},\nwrong SOPClassUID: {}='{}',\nexpecting an 'RT Ion Plan Storage' file instead.".format( + self._rpfp, self._rp.SOPClassUID, sop_class_name + ) + ) + missing_attrs = list() + for a in ["IonBeamSequence"] + self.plan_req_attrs + self.patient_attrs: + b = a.replace(" ", "") + if not hasattr(self._rp, b): + missing_attrs.append(b) + if missing_attrs: + raise IOError( + "bad plan file {},\nmissing keys: {}".format( + self._rpfp, ", ".join(missing_attrs) + ) + ) + # self._get_rds() + if hasattr(self._rp, "DoseReferenceSequence"): + sequence_check(self._rp, "DoseReferenceSequence", 1, 1) + if hasattr(self._rp.DoseReferenceSequence[0], "ReferencedROINumber"): + self._dose_roinumber = int( + self._rp.DoseReferenceSequence[0].ReferencedROINumber + ) + if self._dose_roinumber is None: + logger.info( + "no target ROI specified (probably because of missing DoseReferenceSequence)" + ) + sequence_check(self._rp, "IonBeamSequence", 1, 0) + sequence_check(self._rp, "FractionGroupSequence", 1, 1) + frac0 = self._rp.FractionGroupSequence[0] + sequence_check( + frac0, + "ReferencedBeamSequence", + len(self._rp.IonBeamSequence), + len(self._rp.IonBeamSequence), + ) + number_set = set() + for dcmbeam in self._rp.IonBeamSequence: + nr = int(dcmbeam.BeamNumber) + if nr < 0: + self._beam_numbers_corrupt = True + logger.error( + "CORRUPT INPUT: found a negative beam number {}.".format(nr) + ) + if nr in number_set: + self._beam_numbers_corrupt = True + logger.error( + "CORRUPT INPUT: found same beam number {} for multiple beams.".format( + nr + ) + ) + number_set.add(nr) + if self._beam_numbers_corrupt: + msg = "Beam numbers are corrupt! Will override them with simple enumeration, keep fingers crossed." + logger.error(msg) + self._warnings.append(msg) + logger.debug("checked planfile, looks like all attributes are available...") + + @property + def mswtot(self): + return sum([b.mswtot for b in self._beams]) + + @property + def name(self): + # It looks like RTPlanLabel is for the beamset, + # and RTPlanName is for the entire plan including possibly several beamsets. + # the RTPlanName is not exported anymore in RayStation 8b, so let's forget about the plan name + # some anonymization methods are butchering all useful labeling information, even the label and name of an RT plan. + return str(self._rp.get("RTPlanLabel", "anonymized")) + + @property + def fields(self): + # GATE synomym for 'beams' + return self._beams + + @property + def beams(self): + return self._beams + + @property + def uid(self): + return self._rpuid + + @property + def beam_angles(self): + return [str(bm.gantry_angle) for bm in self._beams] + + @property + def beam_names(self): + return [str(bm.Name) for bm in self._beams] + + @property + def beam_numbers(self): + return [str(bm.Number) for bm in self._beams] + + @property + def Nfractions(self): + # FIXME: some evil DICOM plan files have "NumberOfFractionsPlanned" equal to zero. Force this to be one, or is this zero somehow meaningful & useful? + nfrac = int(self._rp.FractionGroupSequence[0].NumberOfFractionsPlanned) + if nfrac > 0: + return nfrac + logger.error("Got Nfractions={} ???! Using nfrac=1 instead.".format(nfrac)) + return 1 + + @property + def target_ROI_name(self): + return self._dose_roiname + + @target_ROI_name.setter + def target_ROI_name(self, roiname): + self._dose_roiname = roiname + + @property + def target_ROI_number(self): + return self._dose_roinumber + + @property + def prescription_dose(self): + if hasattr(self._rp, "DoseReferenceSequence"): + if hasattr(self._rp.DoseReferenceSequence[0], "TargetPrescriptionDose"): + return float(self._rp.DoseReferenceSequence[0].TargetPrescriptionDose) + return "NA" + + @property + def plan_label(self): + return str(self._rp.get("RTPlanLabel", "anonymized")) + + @property + def sanitized_plan_label(self): + badchars = re.compile("[^a-zA-Z0-9_]") + return re.sub(badchars, "_", self.plan_label) + + @property + def patient_info(self): + return dict( + [ + (a, str(getattr(self._rp, a.replace(" ", "")))) + for a in self.patient_attrs + ] + ) + + @property + def plan_info(self): + reqs = dict( + [ + (a, str(getattr(self._rp, a.replace(" ", "")))) + for a in self.plan_req_attrs + ] + ) + opts = dict( + [ + ( + a, + "NA" + if not hasattr(self._rp, a.replace(" ", "")) + else str(getattr(self._rp, a.replace(" ", ""))), + ) + for a in self.plan_opt_attrs + ] + ) + reqs.update(opts) + return reqs + + @property + def bs_info(self): + info = dict( + [ + ("Number Of Beams", str(len(self._beams))), + ("RT Plan Label", self.plan_label), + ("Prescription Dose", str(self.prescription_dose)), + ("Target ROI Name", str(self.target_ROI_name)), + ( + "Radiation Type", + ", ".join( + set([str(beam.RadiationType[0]) for beam in self._beams]) + ), + ), + ( + "Radiation Type Opengate", + ", ".join( + set([str(beam.RadiationType[1]) for beam in self._beams]) + ), + ), + ( + "Treatment Machine(s)", + ", ".join( + set([str(beam.TreatmentMachineName) for beam in self._beams]) + ), + ), + ] + ) + dirty = self.plan_label + sanitized = self.sanitized_plan_label + if dirty != sanitized: + info["Sanitized RT Plan Label"] = sanitized + return info + + def __repr__(self): + s = "\nPLAN\n\t" + "\n\t".join( + ["{0:30s}: {1}".format(a, self.plan_info[a]) for a in self.plan_attrs] + ) + s += "\nPATIENT\n\t" + "\n\t".join( + ["{0:30s}: {1}".format(a, self.patient_info[a]) for a in self.patient_attrs] + ) + s += "\nBEAMSET\n\t" + "\n\t".join( + ["{0:30s}: {1}".format(a, self.bs_info[a]) for a in self.bs_attrs] + ) + return s + + +def spots_info_from_txt(txtFile, ionType): + # initialize empty variables + nFields = 0 + ntot = [] + energies = [] + nSpots = [] + spots = [] + start_index = [] + G = 0 + + # read content + with open(txtFile, "r") as f: + lines = f.readlines() + + # get plan's info + for i, line in enumerate(lines): + if line.startswith("###GantryAngle"): + l = lines[i + 1].split("\n")[0] + G = int(l) + if line.startswith("##NumberOfFields"): + l = lines[i + 1].split("\n")[0] + nFields = int(l) + if line.startswith("###FinalCumulativeMeterSetWeight"): + l = lines[i + 1].split("\n")[0] + ntot.append(float(l)) + if line.startswith("####Energy"): + l = lines[i + 1].split("\n")[0] + energies.append(float(l)) + if line.startswith("####NbOfScannedSpots"): + l = lines[i + 1].split("\n")[0] + nSpots.append(int(l)) + if line.startswith("####X Y Weight"): + start_index.append(i + 1) + + np = sum(ntot) + for k in range(nFields): + # np = ntot[k] + for i in range(len(energies)): + e = energies[i] + print(f"ENERGY: {e}") + start = start_index[i] + end = start_index[i] + nSpots[i] + for j in range(start, end): + l = lines[j].split("\n")[0].split() + spot = spot_info(float(l[0]), float(l[1]), float(l[2]), e) + spot.beamFraction = float(l[2]) / np + spot.particle_name = ionType + spots.append(spot) + + return spots, np, energies, G + + +def get_spots_from_beamset(beamset): + rad_type = beamset.bs_info["Radiation Type Opengate"] + spots_array = [] + mswtot = beamset.mswtot + for beam in beamset.beams: + # mswtot = beam.mswtot + for energy_layer in beam.layers: + for spot in energy_layer.spots: + nPlannedSpot = spot.w + spot.beamFraction = ( + nPlannedSpot / mswtot + ) # nr particles planned for the spot/tot particles planned for the beam + spot.particle_name = rad_type + spots_array.append(spot) + return spots_array + + +# vim: set et softtabstop=4 sw=4 smartindent: diff --git a/opengate/helpers_tests.py b/opengate/helpers_tests.py index ed9bc540d..46ef0ceb9 100644 --- a/opengate/helpers_tests.py +++ b/opengate/helpers_tests.py @@ -209,7 +209,7 @@ def plot_img_y(ax, img, label): data = itk.GetArrayViewFromImage(img) y = np.sum(data, 2) y = np.sum(y, 0) - x = np.arange(len(y)) * img.GetSpacing()[2] + x = np.arange(len(y)) * img.GetSpacing()[1] ax.plot(x, y, label=label) ax.legend() @@ -219,7 +219,7 @@ def plot_img_x(ax, img, label): data = itk.GetArrayViewFromImage(img) y = np.sum(data, 1) y = np.sum(y, 0) - x = np.arange(len(y)) * img.GetSpacing()[2] + x = np.arange(len(y)) * img.GetSpacing()[0] ax.plot(x, y, label=label) ax.legend() @@ -975,10 +975,13 @@ def plot_gauss_fit(positionVec, dose, fit, show=False): return fig -def create_position_vector(length, spacing): +def create_position_vector(length, spacing, centered=True): # cretae position vector, with origin in the image plane's center width = length * spacing - positionVec = np.arange(0, width, spacing) - width / 2 + spacing / 2 + if centered: + positionVec = np.arange(0, width, spacing) - width / 2 + spacing / 2 + else: + positionVec = np.arange(0, width, spacing) return positionVec @@ -1007,14 +1010,31 @@ def read_mhd(filename): return data, spacing, shape -def create_2D_Edep_colorMap(filepath, show=False): +def plot2D(twodarray, label, show=False): + fig = plt.figure(figsize=(20, 20)) + ax = fig.add_subplot(111) + ax.set_title(label) + plt.imshow(twodarray) + ax.set_aspect("equal") + plt.colorbar(orientation="vertical") + if show: + plt.show() + return fig + + +def create_2D_Edep_colorMap(filepath, show=False, axis="z"): img = itk.imread(str(filepath)) data = itk.GetArrayViewFromImage(img) fig = plt.figure(figsize=(20, 20)) ax = fig.add_subplot(111) ax.set_title("colorMap") - plt.imshow(data[0, :, :]) + if axis == "z": + plt.imshow(data[0, :, :]) + elif axis == "x": + plt.imshow(data[:, :, 0]) + else: + plt.imshow(data[:, 0, :]) ax.set_aspect("equal") plt.colorbar(orientation="vertical") if show: @@ -1192,6 +1212,229 @@ def test_weights(expected_ratio, mhd_1, mhd_2, thresh=0.1): return is_ok +def test_tps_spot_size_positions(data, ref, spacing, thresh=0.1, abs_tol=0.3): + if not np.array_equal(data.size, ref.size): + print("Images do not have the same size") + return False + ok = True + # beam along x + size = data.shape + # get gaussian fit of Edep only around the i-th spot + param_y_out, _, _ = gate.extract_gauss_param_1D(data, size[1], spacing[1], axis=0) + param_y_ref, _, _ = gate.extract_gauss_param_1D(ref, size[1], spacing[1], axis=0) + + param_z_out, _, _ = gate.extract_gauss_param_1D(data, size[0], spacing[2], axis=1) + param_z_ref, _, _ = gate.extract_gauss_param_1D(ref, size[0], spacing[2], axis=1) + + # check positions + print("Check position of the spot") + print(f" opengate: ({param_y_out[1]:.2f},{param_z_out[1]:.2f})") + print(f" gate: ({param_y_ref[1]:.2f},{param_z_ref[1]:.2f})") + + diffmY = param_y_out[1] - param_y_ref[1] # / param_y_ref[1] + diffmZ = param_z_out[1] - param_z_ref[1] # / param_z_ref[1] + mean_diff = np.mean([diffmY, diffmZ]) + + if ( + (abs(diffmY) > 2 * abs_tol) + or (abs(diffmZ) > 2 * abs_tol) + or (abs(mean_diff) > abs_tol) + ): + print( + f"\033[91m Position error above threshold. DiffX={diffmY:.2f}, diffY={diffmZ:.2f}, threshold is 0.3mm \033[0m" + ) + ok = False + + # check sizes + print("Check size of the spot") + print(f" opengate: ({param_y_out[2]:.2f},{param_z_out[2]:.2f})") + print(f" gate: ({param_y_ref[2]:.2f},{param_z_ref[2]:.2f})") + + diffsY = (param_y_out[2] - param_y_ref[2]) / param_y_ref[2] + diffsZ = (param_z_out[2] - param_z_ref[2]) / param_z_ref[2] + + if (diffsY > thresh) or (diffsZ > thresh): + print("\033[91m Size error above threshold \033[0m") + ok = False + + return ok + + +def scale_dose(path, scaling, outpath): + img_mhd_in = itk.imread(path) + data = itk.GetArrayViewFromImage(img_mhd_in) + dose = data * scaling + spacing = img_mhd_in.GetSpacing() + img = gate.itk_image_view_from_array(dose) + img.SetSpacing(spacing) + itk.imwrite(img, outpath) + return outpath + + +def check_dose_grid_geometry(dose_mhd_path, dose_actor): + img = itk.imread(dose_mhd_path) + data = itk.GetArrayViewFromImage(img) + shape = data.shape + spacing = img.GetSpacing() + shape_ref = tuple(np.flip(dose_actor.size)) + spacing_ref = dose_actor.spacing + + ok = True + if shape != shape_ref: + print(f"{shape=} not the same as {shape_ref=}!") + ok = False + + if spacing != spacing_ref: + print(f"{spacing=} not the same as {spacing_ref=}!") + ok = False + + return ok + + +def arangeDx(dx, xV, includeUB=False, lb=[], ub=[]): + if not lb: + lb = np.amin(xV) + if not ub: + ub = np.amax(xV) + if includeUB: + x_int = np.arange(lb, ub + dx / 10, dx) + else: + x_int = np.arange(lb, ub, dx) + return x_int + + +def interpolate1Dprofile(xV, dV, dx=0.01, interpolMethod="cubic"): + f = scipy.interpolate.interp1d( + xV, dV, kind=interpolMethod, fill_value="extrapolate" + ) + xVfine = arangeDx(dx, xV, includeUB=True, lb=np.amin(xV), ub=np.amax(xV)) + dVfine = f(xVfine) + return xVfine, dVfine + + +def getRange(xV, dV, percentLevel=0.8): + dx = 0.01 + xVfine, dVfine = interpolate1Dprofile(xV, dV, dx, "cubic") + + indMaxFine = np.argmax(dVfine) + indR80 = np.argmax( + np.logical_and( + xVfine > xVfine[indMaxFine], dVfine <= percentLevel * dVfine[indMaxFine] + ) + ) + r80 = xVfine[indR80] + dAtR80 = dVfine[indR80] + + return (r80, dAtR80) + + +def get_range_from_image(volume, shape, spacing, axis="y"): + x1, d1 = get_1D_profile(volume, shape, spacing, axis=axis) + r, _ = getRange(x1, d1) + + return r + + +def compareRange( + volume1, + volume2, + shape1, + shape2, + spacing1, + spacing2, + axis1="y", + axis2="y", + thresh=2.0, +): + ok = True + x1, d1 = get_1D_profile(volume1, shape1, spacing1, axis=axis1) + x2, d2 = get_1D_profile(volume2, shape2, spacing2, axis=axis2) + + print("---RANGE80---") + r1, _ = getRange(x1, d1) + r2, _ = getRange(x2, d2) + print(r1) + print(r2) + diff = abs(r2 - r1) + + if diff > thresh: + print(f"\033[91mRange difference is {diff}mm, threshold is {thresh}mm \033[0m") + ok = False + + return ok + + +def get_1D_profile(data, shape, spacing, axis="z"): + if axis == "x": + d1 = np.sum(np.sum(data, 1), 0) + x1 = create_position_vector(shape[2], spacing[0], centered=False) + + if axis == "y": + d1 = np.sum(np.sum(data, 2), 0) + x1 = create_position_vector(shape[1], spacing[1], centered=False) + + if axis == "z": + d1 = np.sum(np.sum(data, 2), 1) + x1 = create_position_vector(shape[0], spacing[2], centered=False) + + return x1, d1 + + +def compare_dose_at_points( + pointsV, + dose1, + dose2, + shape1, + shape2, + spacing1, + spacing2, + axis1="z", + axis2="z", + rel_tol=0.03, +): + ok = True + s1 = 0 + s2 = 0 + x1, doseV1 = get_1D_profile(dose1, shape1, spacing1, axis=axis1) + x2, doseV2 = get_1D_profile(dose2, shape2, spacing2, axis=axis2) + # plt.plot(x1, doseV1) + # plt.plot(x2, doseV2) + # plt.show() + for p in pointsV: + # get dose at the position p [mm] + cp1 = min(x1, key=lambda x: abs(x - p)) + d1_p = doseV1[np.where(x1 == cp1)] + + cp2 = min(x2, key=lambda x: abs(x - p)) + d2_p = doseV2[np.where(x2 == cp2)] + + s1 += d1_p + s2 += d2_p + + print(abs(s1 - s2) / s2) + + # print(f"Dose difference at {p} mm is {diff_pc}%") + if abs(s1 - s2) / s2 > rel_tol: + print(f"\033[91mDose difference above threshold \033[0m") + ok = False + return ok + + +def assert_img_sum(img1, img2, sum_tolerance=5): + data1 = itk.GetArrayViewFromImage(img1).ravel() + data2 = itk.GetArrayViewFromImage(img2).ravel() + + s1 = np.sum(data1) + s2 = np.sum(data2) + if s1 == 0 and s2 == 0: + t = 0 + else: + t = np.fabs((s1 - s2) / s1) * 100 + b = t < sum_tolerance + print_test(b, f"Img sums {s1} vs {s2} : {t:.2f} % (tol {sum_tolerance:.2f} %)") + return b + + def check_diff(value1, value2, tolerance, txt): diff = np.fabs(value1 - value2) / value1 * 100 t = diff < tolerance diff --git a/opengate/source/TreatmentPlanSource.py b/opengate/source/TreatmentPlanSource.py new file mode 100644 index 000000000..6089fcc58 --- /dev/null +++ b/opengate/source/TreatmentPlanSource.py @@ -0,0 +1,135 @@ +import numpy as np +from scipy.spatial.transform import Rotation +import opengate as gate + + +class TreatmentPlanSource: + def __init__(self, name, sim): + self.name = name + # self.mother = None + self.rotation = Rotation.identity() + self.translation = [0, 0, 0] + self.spots = None + self.beamline_model = None + self.n_sim = 0 + self.sim = sim # simulation obj to which we want to add the tp source + + def __del__(self): + pass + + def set_particles_to_simulate(self, n_sim): + self.n_sim = n_sim + + def set_spots(self, spots): + self.spots = spots + + def set_spots_from_rtplan(self, rt_plan_path): + beamset = gate.beamset_info(rt_plan_path) + gantry_angle = beamset.beam_angles[0] + spots = gate.get_spots_from_beamset(beamset) + self.spots = spots + self.rotation = Rotation.from_euler("z", gantry_angle, degrees=True) + + def set_beamline_model(self, beamline): + self.beamline_model = beamline + + def initialize_tpsource(self): + # some alias + spots_array = self.spots + sim = self.sim + nSim = self.n_sim + beamline = self.beamline_model + self.d_nozzle_to_iso = beamline.distance_nozzle_iso + self.d_stearMag_to_iso_x = beamline.distance_stearmag_to_isocenter_x + self.d_stearMag_to_iso_y = beamline.distance_stearmag_to_isocenter_y + + # mapping factors between iso center plane and nozzle plane (due to steering magnets) + cal_proportion_factor = ( + lambda d_magnet_iso: 1 + if (d_magnet_iso == float("inf")) + else (d_magnet_iso - self.d_nozzle_to_iso) / d_magnet_iso + ) + self.proportion_factor_x = cal_proportion_factor(self.d_stearMag_to_iso_x) + self.proportion_factor_y = cal_proportion_factor(self.d_stearMag_to_iso_y) + + # initialize a pencil beam for each spot + for i, spot in enumerate(spots_array): + nspot = np.round( + spot.beamFraction * nSim + ) # simulate a fraction of the beam particles for this spot + if nspot == 0: + continue + + source = sim.add_source("PencilBeamSource", f"{self.name}_spot_{i}") + + # set energy + source.energy.type = "gauss" + source.energy.mono = beamline.get_energy(nominal_energy=spot.energy) + source.energy.sigma_gauss = beamline.get_sigma_energy( + nominal_energy=spot.energy + ) + + source.particle = spot.particle_name + source.position.type = "disc" # pos = Beam, shape = circle + sigma + + # # set mother + # if self.mother is not None: + # source.mother = self.mother + + # POSITION: + source.position.translation = self._get_pbs_position(spot) + + # ROTATION: + source.position.rotation = self._get_pbs_rotation(spot) + + # add weight + # source.weight = -1 + source.n = nspot + + # set optics parameters + source.direction.partPhSp_x = [ + beamline.get_sigma_x(spot.energy), + beamline.get_theta_x(spot.energy), + beamline.get_epsilon_x(spot.energy), + beamline.conv_x, + ] + source.direction.partPhSp_y = [ + beamline.get_sigma_y(spot.energy), + beamline.get_theta_y(spot.energy), + beamline.get_epsilon_y(spot.energy), + beamline.conv_y, + ] + + def _get_pbs_position(self, spot): + # (x,y) referr to isocenter plane. + # Need to be corrected to referr to nozzle plane + pos = [ + (spot.xiec) * self.proportion_factor_x, + (spot.yiec) * self.proportion_factor_y, + self.d_nozzle_to_iso, + ] + # Gantry angle = 0 -> source comes from +y and is positioned along negative side of y-axis + # https://opengate.readthedocs.io/en/latest/source_and_particle_management.html + + position = (self.rotation * Rotation.from_euler("x", np.pi / 2)).apply( + pos + ) + self.translation + + return position + + def _get_pbs_rotation(self, spot): + # by default the source points in direction z+. + # Need to account for SM direction deviation and rotation thoward isocenter (270 deg around x) + # then rotate of gantry angle + rotation = [0.0, 0.0, 0.0] + beta = np.arctan(spot.yiec / self.d_stearMag_to_iso_y) + alpha = np.arctan(spot.xiec / self.d_stearMag_to_iso_x) + rotation[0] = -np.pi / 2 + beta + rotation[2] = -alpha + + # apply gantry angle + spot_rotation = ( + self.rotation * Rotation.from_euler("xyz", rotation) + ).as_matrix() + + return spot_rotation diff --git a/opengate/tests/src/test059_TPSource_abs_dose.py b/opengate/tests/src/test059_TPSource_abs_dose.py new file mode 100644 index 000000000..89a319e23 --- /dev/null +++ b/opengate/tests/src/test059_TPSource_abs_dose.py @@ -0,0 +1,194 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +rad = gate.g4_units("rad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.translation = [1148 * mm, 0.0, 0.0] +box.rotation = Rotation.from_euler("y", -90, degrees=True).as_matrix() +box.material = "Vacuum" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# target +phantom = sim.add_volume("Box", "phantom") +phantom.size = [500 * mm, 500 * mm, 400 * mm] +phantom.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() +phantom.translation = [-200.0, 0.0, 0] +phantom.material = "G4_WATER" +phantom.color = [0, 0, 1, 1] + +# roos chamber +roos = sim.add_volume("Tubs", "roos") +roos.mother = phantom.name +roos.material = "G4_WATER" +roos.rmax = 7.8 +roos.rmin = 0 +roos.dz = 200 +roos.color = [1, 0, 1, 1] + + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" #'QGSP_BIC_HP_EMZ' #"FTFP_INCLXX_EMZ" + +sim.set_cut("world", "all", 1000 * km) + + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "abs_dose_roos.mhd" +dose.mother = roos.name +dose.size = [1, 1, 800] +dose.spacing = [15.6, 15.6, 0.5] +dose.hit_type = "random" +dose.gray = True + + +## ---------- DEFINE BEAMLINE MODEL -------------## +IR2HBL = gate.BeamlineModel() +IR2HBL.name = None +IR2HBL.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +IR2HBL.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_x = 6700.00 +# SMY to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_y = 7420.00 +# polinomial coefficients +IR2HBL.energy_mean_coeffs = [11.91893485094217, -9.539517997860457] +IR2HBL.energy_spread_coeffs = [0.0004790681841295621, 5.253257865904452] +IR2HBL.sigma_x_coeffs = [2.3335753978880014] +IR2HBL.theta_x_coeffs = [0.0002944903217664001] +IR2HBL.epsilon_x_coeffs = [0.0007872786903040108] +IR2HBL.sigma_y_coeffs = [1.9643343053823967] +IR2HBL.theta_y_coeffs = [0.0007911780133478402] +IR2HBL.epsilon_y_coeffs = [0.0024916149017600447] + +## --------START PENCIL BEAM SCANNING---------- ## +# NOTE: HBL means that the beam is coming from -x (90 degree rot around y) + +nSim = 200000 # 328935 # particles to simulate per beam + +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan4Gate-F5x5cm_E120MeVn.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("RT_plan", sim) +tps.set_beamline_model(IR2HBL) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", G, degrees=True) +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +## -------------END SCANNING------------- ## +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## +dose_path = gate.scale_dose( + str(dose.output).replace(".mhd", "_dose.mhd"), + ntot / nSim, + output_path / "threeDdoseWaternew.mhd", +) + +# ABSOLUTE DOSE + +# read output and ref +img_mhd_out = itk.imread(dose_path) +img_mhd_ref = itk.imread( + ref_path / "idc-PHANTOM-roos-F5x5cm_E120MeVn-PLAN-Physical.mhd" +) +data = itk.GetArrayViewFromImage(img_mhd_out) +data_ref = itk.GetArrayViewFromImage(img_mhd_ref) +shape = data.shape +spacing = img_mhd_out.GetSpacing() +spacing_ref = np.flip(img_mhd_ref.GetSpacing()) + +ok = gate.assert_img_sum( + img_mhd_out, + img_mhd_ref, +) + +points = 400 - np.linspace(10, 14, 9) +ok = ( + gate.compare_dose_at_points( + points, + data, + data_ref, + shape, + data_ref.shape, + spacing, + spacing_ref, + axis1="z", + axis2="x", + rel_tol=0.03, + ) + and ok +) + +# 1D +# fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +# gate.plot_img_axis(ax, img_mhd_out, "x profile", axis="z") +# gate.plot_img_axis(ax, img_mhd_ref, "x ref", axis="z") +# plt.show() + +# fig.savefig(output_path / "dose_profiles_water.png") + + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_gantry_rot.py b/opengate/tests/src/test059_TPSource_gantry_rot.py new file mode 100644 index 000000000..ec3d7f595 --- /dev/null +++ b/opengate/tests/src/test059_TPSource_gantry_rot.py @@ -0,0 +1,201 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +## ---------- DEFINE BEAMLINE MODEL -------------## +beamline = gate.BeamlineModel() +beamline.name = None +beamline.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +beamline.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +beamline.distance_stearmag_to_isocenter_x = 6700.00 +# SMY to Isocenter distance +beamline.distance_stearmag_to_isocenter_y = 7420.00 +# polinomial coefficients +beamline.energy_mean_coeffs = [11.91893485094217, -9.539517997860457] +beamline.energy_spread_coeffs = [0.0004790681841295621, 5.253257865904452] +beamline.sigma_x_coeffs = [2.3335753978880014] +beamline.theta_x_coeffs = [0.0002944903217664001] +beamline.epsilon_x_coeffs = [0.0007872786903040108] +beamline.sigma_y_coeffs = [1.9643343053823967] +beamline.theta_y_coeffs = [0.0007911780133478402] +beamline.epsilon_y_coeffs = [0.0024916149017600447] + +# NOTE: HBL means that the beam is coming from -x (90 degree rot around y) +nSim = 20000 # particles to simulate per beam +# rt_plan = ref_path / "RP1.2.752.243.1.1.20220406175810679.4500.52008_tagman.dcm" +# beamset = gate.beamset_info(rt_plan) +# G = float(beamset.beam_angles[0]) + + +## ---- VBL Nozzle --- +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.rotation = Rotation.from_euler("x", -90, degrees=True).as_matrix() +box.translation = [0.0, -1148 * mm, 0.0] +box.material = "Vacuum" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# Rashi +rashi = sim.add_volume("Box", "rashi") +rashi.mother = box.name +rashi.size = [500 * mm, 500 * mm, 5 * mm] +rashi.translation = [0.0, 0.0, 200 * mm] +rashi.material = "G4_LUCITE" +rashi.color = [1, 0, 1, 1] + +## ---- HBL Nozzle --- +box_rot = sim.add_volume("Box", "box_rot") +gate.copy_user_info(box, box_rot) +box_rot.rotation = Rotation.from_euler("y", -90, degrees=True).as_matrix() +box_rot.translation = [1148.0, 0.0, 1000.0] + +nozzle_rot = sim.add_volume("Box", "nozzle_rot") +gate.copy_user_info(nozzle, nozzle_rot) +nozzle_rot.mother = box_rot.name + +rashi_rot = sim.add_volume("Box", "rashi_rot") +gate.copy_user_info(rashi, rashi_rot) +rashi_rot.mother = box_rot.name + +# ----------------------------------- + +# target 1 VBL +phantom = sim.add_volume("Box", "phantom") +phantom.size = [324 * mm, 324 * mm, 324 * mm] +phantom.translation = [0 * mm, 0.0, 0.0] +phantom.material = "G4_WATER" +phantom.color = [0, 0, 1, 1] + +# target 2 HBL +phantom_rot = sim.add_volume("Box", "phantom_rot") +gate.copy_user_info(phantom, phantom_rot) +phantom_rot.rotation = Rotation.from_euler("z", 90, degrees=True).as_matrix() +phantom_rot.translation = [0.0, 0.0, 1000.0] + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "testTPSgantry.mhd" +dose.mother = phantom.name +dose.size = [162, 1620, 162] +dose.spacing = [2.0, 0.2, 2.0] +dose.hit_type = "random" +dose.gray = True + +dose_rot = sim.add_actor("DoseActor", "doseInXYZ_rot") +gate.copy_user_info(dose, dose_rot) +dose_rot.mother = phantom_rot.name +dose_rot.output = output_path / "testTPSganry_rot.mhd" + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" +sim.set_cut("world", "all", 1000 * km) + +# add TPSources +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan4Gate-1D_HBL_120.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("VBL", sim) +tps.set_beamline_model(beamline) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", 0, degrees=True) +tps.initialize_tpsource() + +tps_rot = gate.TreatmentPlanSource("HBL", sim) +tps_rot.set_beamline_model(beamline) +tps_rot.set_particles_to_simulate(nSim) +tps_rot.set_spots(spots) +tps_rot.rotation = Rotation.from_euler("z", G, degrees=True) +tps_rot.translation = [0.0, 0.0, 1000.0] +tps_rot.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## + +# ABSOLUTE DOSE +# ok = gate.assert_images( +# dose.output, +# dose_rot.output, +# stat, +# tolerance=50, +# ignore_value=0, +# ) +ok = True + +# read output and ref +img_mhd_out = itk.imread(dose_rot.output) +img_mhd_ref = itk.imread(dose.output) +data = itk.GetArrayViewFromImage(img_mhd_out) +data_ref = itk.GetArrayViewFromImage(img_mhd_ref) +spacing = img_mhd_out.GetSpacing() + +# Range 80 +ok = ( + gate.compareRange(data, data_ref, data.shape, data_ref.shape, spacing, spacing) + and ok +) + +# 1D plots +# fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +# gate.plot_img_axis(ax, img_mhd_out, "y profile", axis="y") +# gate.plot_img_axis(ax, img_mhd_ref, "y ref", axis="y") +# fig.savefig(output_path / "dose_profiles_water.png") +# plt.show() + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_optics.py b/opengate/tests/src/test059_TPSource_optics.py new file mode 100644 index 000000000..be2152c28 --- /dev/null +++ b/opengate/tests/src/test059_TPSource_optics.py @@ -0,0 +1,225 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +rad = gate.g4_units("rad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.translation = [1148 * mm, 0.0, 0.0] +box.rotation = Rotation.from_euler("y", -90, degrees=True).as_matrix() +box.material = "Vacuum" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# target +phantom = sim.add_volume("Box", "phantom") +phantom.size = [300 * mm, 310 * mm, 310 * mm] +phantom.translation = [150 * mm, 0.0, 0.0] +# phantom.rotation = Rotation.from_euler('x',-90,degrees=True).as_matrix() +phantom.material = "G4_AIR" +phantom.color = [0, 0, 1, 1] + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" # "Shielding_EMZ"#"FTFP_INCLXX_EMZ" + +# p.physics_list_name = "QGSP_BIC_EMZ" +sim.set_cut("world", "all", 1000 * km) + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "testTPSoptics.mhd" +dose.mother = phantom.name +dose.size = [30, 620, 620] +dose.spacing = [10.0, 0.5, 0.5] +dose.hit_type = "random" +dose.gray = True + +## ---------- DEFINE BEAMLINE MODEL -------------## +IR2HBL = gate.BeamlineModel() +IR2HBL.name = None +IR2HBL.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +IR2HBL.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_x = 6700.00 +# SMY to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_y = 7420.00 +# polinomial coefficients +IR2HBL.energy_mean_coeffs = [11.91893485094217, -9.539517997860457] +IR2HBL.energy_spread_coeffs = [0.0004790681841295621, 5.253257865904452] +IR2HBL.sigma_x_coeffs = [-0.00011142901344618727, 2.346946879501544] +IR2HBL.theta_x_coeffs = [-3.6368814874049214e-07, 0.0003381328996152591] +IR2HBL.epsilon_x_coeffs = [3.1292233857396716e-06, 0.0004117718840152502] +IR2HBL.sigma_y_coeffs = [-0.0004009682717802152, 2.0124504979960225] +IR2HBL.theta_y_coeffs = [-8.437400716390318e-07, 0.000892426821944524] +IR2HBL.epsilon_y_coeffs = [-8.757558864087579e-08, 0.00250212397239695] + +## --------START PENCIL BEAM SCANNING---------- ## +# NOTE: HBL means that the beam is coming from -x (90 degree rot around y) +# nSim = 328935 # particles to simulate per beam +nSim = 20000 +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan4Gate-gate_test51_TP_1_old.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("RT_plan", sim) +tps.set_beamline_model(IR2HBL) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", G, degrees=True) +# rt_plan = ref_path / "RP1.2.752.243.1.1.20230119115736709.2000.75541.dcm" +# tps.set_spots_from_rtplan(rt_plan) # no need to set rotation here +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +## -------------END SCANNING------------- ## +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## +dose_path = gate.scale_dose( + str(dose.output).replace(".mhd", "_dose.mhd"), + ntot / nSim, + output_path / "threeDdoseAirSpots.mhd", +) + +# SPOT POSITIONS COMPARISON +# read output and ref +img_mhd_out = itk.imread(dose_path) +img_mhd_ref = itk.imread( + ref_path / "idc-PHANTOM-air_box-gate_test51_TP_1-PLAN-Physical.mhd" +) +data = itk.GetArrayViewFromImage(img_mhd_out) +data_ref = itk.GetArrayViewFromImage(img_mhd_ref) +shape = data.shape +spacing = img_mhd_out.GetSpacing() + +# spot comparison (N.B x and z are inverted in np array!) +# spots in the plan file +yz = [ + 50, + -50, + -50, + -50, + -100, + 0, + 100, + 0, + 50, + 50, + 0, + 100, + -100, + -100, + 0, + -100, + 100, + -100, + 0, + 0, + -50, + 50, + -100, + 100, + 100, + 100, +] + +yzM = np.array(yz).reshape(int(len(yz) / 2), 2) +# convert from mm (wrt image center) to voxel +spot_y = [int(y / dose.spacing[1]) + int(dose.size[1] / 2) for y in yzM[:, 0]] +spot_z = [int(z / dose.spacing[1]) + int(dose.size[1] / 2) for z in yzM[:, 1]] + +thresh = 0.8 ## OSS: need step limiter + +# # 1D +# fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +# gate.plot_img_axis(ax, img_mhd_out, "z profile", axis="z") +# #gate.plot_img_axis(ax, img_mhd_out, "x profile", axis="x") +# gate.plot_img_axis(ax, img_mhd_out, "y profile", axis="y") + + +# # fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +# gate.plot_img_axis(ax, img_mhd_ref, "z ref", axis="z") +# #gate.plot_img_axis(ax, img_mhd_ref, "x ref", axis="x") +# gate.plot_img_axis(ax, img_mhd_ref, "y ref", axis="y") + +# fig.savefig(output_path / "dose_profiles_spots.png") + + +ok = True +for i in range(1, shape[2], shape[2] // 3): + # check i-th slab + print(f"Airslab nr. {i}") + # gate.plot2D(data[:, :, i], "2D Edep opengate", show=True) + # gate.plot2D(data_ref[:, :, i], "2D Edep gate", show=True) + for y, z in zip(spot_y, spot_z): + # i = 0 + print(f" ({y:.2f},{z:.2f})") + # 'cut' the slab around the spot expected in y,z + w = 30 # cut window's half size + d_out = data[z - w : z + w, y - w : y + w, i : i + 1] + d_ref = data_ref[z - w : z + w, y - w : y + w, i : i + 1] + ok = ( + gate.test_tps_spot_size_positions( + d_out, d_ref, spacing, thresh=thresh, abs_tol=1.0 + ) + and ok + ) + + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_optics_vbl.py b/opengate/tests/src/test059_TPSource_optics_vbl.py new file mode 100644 index 000000000..3a58ea53e --- /dev/null +++ b/opengate/tests/src/test059_TPSource_optics_vbl.py @@ -0,0 +1,211 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +rad = gate.g4_units("rad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.rotation = Rotation.from_euler("x", -90, degrees=True).as_matrix() +box.translation = [0.0, -1148 * mm, 0.0] +box.material = "Air" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# target +phantom = sim.add_volume("Box", "phantom") +phantom.size = [300 * mm, 310 * mm, 310 * mm] +phantom.rotation = Rotation.from_euler("z", -90, degrees=True).as_matrix() +phantom.translation = [0.0, -150 * mm, 0.0] +phantom.material = "G4_AIR" +phantom.color = [0, 0, 1, 1] + + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" +# p.physics_list_name = "QGSP_BIC_EMZ" +sim.set_cut("world", "all", 1000 * km) + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "TPS_optics_vbl.mhd" +dose.mother = phantom.name +dose.size = [30, 620, 620] +dose.spacing = [10.0, 0.5, 0.5] +dose.hit_type = "random" +dose.gray = True + +## ---------- DEFINE BEAMLINE MODEL -------------## +IR2VBL = gate.BeamlineModel() +IR2VBL.name = None +IR2VBL.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +IR2VBL.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +IR2VBL.distance_stearmag_to_isocenter_x = 6700.00 +# SMY to Isocenter distance +IR2VBL.distance_stearmag_to_isocenter_y = 7420.00 +# polinomial coefficients +IR2VBL.energy_mean_coeffs = [11.91893485094217, -9.539517997860457] +IR2VBL.energy_spread_coeffs = [0.0004790681841295621, 5.253257865904452] +IR2VBL.sigma_x_coeffs = [-0.00011142901344618727, 2.346946879501544] +IR2VBL.theta_x_coeffs = [-3.6368814874049214e-07, 0.0003381328996152591] +IR2VBL.epsilon_x_coeffs = [3.1292233857396716e-06, 0.0004117718840152502] +IR2VBL.sigma_y_coeffs = [-0.0004009682717802152, 2.0124504979960225] +IR2VBL.theta_y_coeffs = [-8.437400716390318e-07, 0.000892426821944524] +IR2VBL.epsilon_y_coeffs = [-8.757558864087579e-08, 0.00250212397239695] + + +## --------START PENCIL BEAM SCANNING---------- ## +# nSim = 328935 # particles to simulate per beam +nSim = 20000 +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan4Gate-gate_test51tps_v.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("RT_plan", sim) +tps.set_beamline_model(IR2VBL) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", G, degrees=True) +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +## -------------END SCANNING------------- ## +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## +dose_path = gate.scale_dose( + str(dose.output).replace(".mhd", "_dose.mhd"), + ntot / nSim, + output_path / "threeDdoseAirSpots_vbl.mhd", +) + +# SPOT POSITIONS COMPARISON +# read output and ref +img_mhd_out = itk.imread(dose_path) +img_mhd_ref = itk.imread( + ref_path / "idc-PHANTOM-air_box_vbl-gate_test51tps_v-PLAN-Physical.mhd" +) +data = itk.GetArrayViewFromImage(img_mhd_out) +data_ref = itk.GetArrayViewFromImage(img_mhd_ref) +shape = data.shape +spacing = img_mhd_out.GetSpacing() + + +# spot comparison (N.B x and z are inverted in np array!) +# spots in the plan file +yz = [ + 0, + 50, + 0, + 0, + -62.507, + -21.669, + -61.951, + -72.23, + 50, + -50, + 50, + 0, + 50, + 50, + -50, + 50, +] + +yzM = np.array(yz).reshape(int(len(yz) / 2), 2) +# convert from mm (wrt image center) to voxel +spot_y = [int(y / dose.spacing[1]) + int(dose.size[1] / 2) for y in yzM[:, 0]] +spot_z = [int(z / dose.spacing[1]) + int(dose.size[1] / 2) for z in yzM[:, 1]] + +thresh = 0.8 ## OSS: need step limiter + +# 1D +fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +gate.plot_img_axis(ax, img_mhd_out, "z profile", axis="z") +gate.plot_img_axis(ax, img_mhd_out, "x profile", axis="x") +gate.plot_img_axis(ax, img_mhd_out, "y profile", axis="y") + +# fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +gate.plot_img_axis(ax, img_mhd_ref, "z ref", axis="z") +gate.plot_img_axis(ax, img_mhd_ref, "x ref", axis="x") +gate.plot_img_axis(ax, img_mhd_ref, "y ref", axis="y") +fig.savefig(output_path / "dose_profiles_spots_vbl.png") + +ok = True +for i in range(1, shape[2], shape[2] // 3): + # check i-th slab + print(f"Airslab nr. {i}") + # gate.plot2D(data[:, :, i], "2D Edep opengate", show=True) + # gate.plot2D(data_ref[:, :, i], "2D Edep gate", show=True) + for y, z in zip(spot_y, spot_z): + # i = 0 + print(f" ({y:.2f},{z:.2f})") + # 'cut' the slab around the spot expected in y,z + w = 30 # cut window's half size + d_out = data[z - w : z + w, y - w : y + w, i : i + 1] + d_ref = data_ref[z - w : z + w, y - w : y + w, i : i + 1] + ok = ( + gate.test_tps_spot_size_positions( + d_out, d_ref, spacing, thresh=thresh, abs_tol=1.0 + ) + and ok + ) + + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_range_ref.py b/opengate/tests/src/test059_TPSource_range_ref.py new file mode 100644 index 000000000..ef4d641eb --- /dev/null +++ b/opengate/tests/src/test059_TPSource_range_ref.py @@ -0,0 +1,161 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +rad = gate.g4_units("rad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.translation = [1148 * mm, 0.0, 0.0] +box.rotation = Rotation.from_euler("y", -90, degrees=True).as_matrix() +box.material = "Vacuum" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# target +phantom = sim.add_volume("Box", "phantom") +phantom.size = [500 * mm, 500 * mm, 400 * mm] +phantom.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() +phantom.translation = [-200.0, 0.0, 0] +phantom.material = "G4_WATER" +phantom.color = [0, 0, 1, 1] + +# roos chamber +peak_finder = sim.add_volume("Tubs", "peak_finder") +peak_finder.mother = phantom.name +peak_finder.material = "G4_WATER" +peak_finder.rmax = 40.3 +peak_finder.rmin = 0 +peak_finder.dz = 200 +peak_finder.color = [1, 0, 1, 1] + + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" #'QGSP_BIC_HP_EMZ' #"FTFP_INCLXX_EMZ" +sim.set_cut("world", "all", 1000 * km) + + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "dose_peak_finder.mhd" +dose.mother = peak_finder.name +dose.size = [1, 1, 8000] +dose.spacing = [80.6, 80.6, 0.05] +dose.hit_type = "random" +dose.gray = True + + +## ---------- DEFINE BEAMLINE MODEL -------------## +IR2HBL = gate.BeamlineModel() +IR2HBL.name = None +IR2HBL.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +IR2HBL.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_x = 2000.00 +# SMY to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_y = 2000.00 +# polinomial coefficients +IR2HBL.energy_mean_coeffs = [12.0, -9.54] +IR2HBL.energy_spread_coeffs = [0.00048, 5.2532] +IR2HBL.sigma_x_coeffs = [2.33] +IR2HBL.theta_x_coeffs = [0.00029] +IR2HBL.epsilon_x_coeffs = [0.00078] +IR2HBL.sigma_y_coeffs = [1.96] +IR2HBL.theta_y_coeffs = [0.00079] +IR2HBL.epsilon_y_coeffs = [0.0024] + +## --------START PENCIL BEAM SCANNING---------- ## +# NOTE: HBL means that the beam is coming from -x (90 degree rot around y) +nSim = 20000 # 328935 # particles to simulate per beam +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "PlanCentralSpot_1440MeV.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("RT_plan", sim) +tps.set_beamline_model(IR2HBL) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", G, degrees=True) +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +## -------------END SCANNING------------- ## +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## +dose_path = str(dose.output).replace(".mhd", "_dose.mhd") + +# RANGE + +# read output and ref +img_mhd_out = itk.imread(dose_path) +data = itk.GetArrayViewFromImage(img_mhd_out) +shape = data.shape +spacing = img_mhd_out.GetSpacing() + +# Range 80 +range80_gate9_E120MeV = 367.06 +range_opengate = gate.get_range_from_image(data, data.shape, spacing, axis="z") + +thresh = 2.0 * mm +ok = True +if abs(range_opengate - range80_gate9_E120MeV) > thresh: + ok = False + + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_weights.py b/opengate/tests/src/test059_TPSource_weights.py new file mode 100644 index 000000000..6f655abcf --- /dev/null +++ b/opengate/tests/src/test059_TPSource_weights.py @@ -0,0 +1,206 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from scipy.spatial.transform import Rotation +import os + +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 123654789 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +mrad = gate.g4_units("mrad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +## FIRST DETECTOR ## +# box +# translation and rotation like in the Gate macro +box1 = sim.add_volume("Box", "box1") +box1.size = [10 * cm, 10 * cm, 105 * cm] +box1.translation = [0 * cm, 0 * cm, 52.5 * cm] +box1.material = "Vacuum" +box1.color = [0, 0, 1, 1] + +# phantoms +m = Rotation.identity().as_matrix() + +phantom = sim.add_volume("Box", "phantom_a_1") +phantom.mother = "box1" +phantom.size = [100 * mm, 100 * mm, 50 * mm] +phantom.translation = [0 * mm, 0 * mm, -500 * mm] +phantom.rotation = m +phantom.material = "G4_AIR" +phantom.color = [1, 0, 1, 1] + + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInYZ_1") +filename = "phantom_a_1.mhd" +dose.output = output_path / filename +dose.mother = "phantom_a_1" +dose.size = [250, 250, 1] +dose.spacing = [0.4, 0.4, 2] +dose.hit_type = "random" + + +## SECOND DETECTOR ## +# box +# translation and rotation like in the Gate macro +box2 = sim.add_volume("Box", "box2") +box2.size = [10 * cm, 10 * cm, 105 * cm] +box2.translation = [30 * cm, 0 * cm, 52.5 * cm] +box2.material = "Vacuum" +box2.color = [0, 0, 1, 1] + +# phantoms +m = Rotation.identity().as_matrix() + +phantom2 = sim.add_volume("Box", "phantom_a_2") +phantom2.mother = "box2" +phantom2.size = [100 * mm, 100 * mm, 50 * mm] +phantom2.translation = [0 * mm, 0 * mm, -500 * mm] +phantom2.rotation = m +phantom2.material = "G4_AIR" +phantom2.color = [1, 0, 1, 1] + + +# add dose actor +dose2 = sim.add_actor("DoseActor", "doseInYZ_2") +filename = "phantom_a_2.mhd" +dose2.output = output_path / filename +dose2.mother = "phantom_a_2" +dose2.size = [250, 250, 1] +dose2.spacing = [0.4, 0.4, 2] +dose2.hit_type = "random" + +## TPS SOURCE ## +# beamline model +beamline = gate.BeamlineModel() +beamline.name = None +beamline.radiation_types = "proton" + +# polinomial coefficients +beamline.energy_mean_coeffs = [1, 0] +beamline.energy_spread_coeffs = [0.4417036946562556] +beamline.sigma_x_coeffs = [2.3335754] +beamline.theta_x_coeffs = [2.3335754e-3] +beamline.epsilon_x_coeffs = [0.00078728e-3] +beamline.sigma_y_coeffs = [1.96433431] +beamline.theta_y_coeffs = [0.00079118e-3] +beamline.epsilon_y_coeffs = [0.00249161e-3] + +# tps +nSim = 60000 # particles to simulate per beam +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan2Spots.txt", "proton" +) +tps = gate.TreatmentPlanSource("test", sim) +tps.set_beamline_model(beamline) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("x", 90, degrees=True) +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" +sim.set_cut("world", "all", 1000 * km) +# sim.set_user_limits("phantom_a_2","max_step_size",1,['proton']) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +# start simulation +sim.run() +output = sim.output + +# print results at the end +stat = output.get_actor("Stats") +print(stat) + + +# ---------------------------------------------------------------------------------------------------------------- +# tests + +# energy deposition: we expect the edep from source two +# to be double the one of source one + +print("Compare tps Edep to single pb sources") +print(" --------------------------------------- ") +mhd_1 = "phantom_a_1.mhd" +mhd_2 = "phantom_a_2.mhd" +test = True + +# check first spot +test = ( + gate.assert_images( + output_path / mhd_1, + ref_path / mhd_1, + stat, + tolerance=70, + ignore_value=0, + ) + and test +) + +# check second spot +test = ( + gate.assert_images( + output_path / mhd_1, + ref_path / mhd_1, + stat, + tolerance=70, + ignore_value=0, + ) + and test +) +print(" --------------------------------------- ") +# fig1 = gate.create_2D_Edep_colorMap(output_path / mhd_1, show=True) +# fig2 = gate.create_2D_Edep_colorMap(output_path / mhd_2, show=True) + +print("Compare ratio of the two spots with expected ratio") + +# Total Edep +is_ok = ( + gate.test_weights( + 2, + output_path / mhd_1, + output_path / mhd_2, + thresh=0.2, + ) + and test +) + + +gate.test_ok(is_ok)