diff --git a/opengate/actors/coincidences.py b/opengate/actors/coincidences.py index a70dcd60b..49c3246a9 100644 --- a/opengate/actors/coincidences.py +++ b/opengate/actors/coincidences.py @@ -938,3 +938,144 @@ def ccmod_make_cones( ) # flatten index Do not save old index from data return data_cones + + +def kill_multiple_coinc(df: pd.DataFrame, group_col: str = "CoincID") -> pd.DataFrame: + """ + Remove coincidences that contain more than 2 singles. + """ + singles_per_coinc = df.groupby(group_col, sort=False).size() + valid_ids = singles_per_coinc[singles_per_coinc == 2].index + return df[df[group_col].isin(valid_ids)].copy() + + +def kill_same_volume_pairs( + df: pd.DataFrame, + group_col: str = "CoincID", + volume_col: str = "PreStepUniqueVolumeID", +) -> pd.DataFrame: + """ + Remove 2‑single coincidences where both singles occur in the same volume. + + Assumes: + - Each CoincID has exactly 2 singles (e.g. after kill_multiple_coinc). + """ + unique_volumes_per_coinc = df.groupby(group_col, sort=False)[volume_col].nunique() + good_ids = unique_volumes_per_coinc[unique_volumes_per_coinc >= 2].index + return df[df[group_col].isin(good_ids)].copy() + + +def _first_tree(f: uproot.ReadOnlyFile): + key = f.keys()[0] + name = key.split(";")[0] + return f[name] + + +def _get_tree(f: uproot.ReadOnlyFile, preferred: str | None): + if preferred and preferred in f: + return f[preferred] + return _first_tree(f) + + +def ccmod_merge_several_singles_root_into_one( + scatt_root_filename, + abs_root_filename, + output_root_filename, + scatt_tree_name, + abs_tree_name, + output_tree_name, + overwrite: bool = True, + save_branches: list[str] | None = None, +): + scatt_root_filename = Path(scatt_root_filename) + abs_root_filename = Path(abs_root_filename) + output_root_filename = Path(output_root_filename) + output_root_filename.parent.mkdir(parents=True, exist_ok=True) + + if output_root_filename.exists(): + if overwrite: + output_root_filename.unlink() + else: + return output_root_filename + + required_branches = ( + "EventID", + "GlobalTime", + "PreStepUniqueVolumeID", + "TotalEnergyDeposit", + "PostPosition_X", + "PostPosition_Y", + "PostPosition_Z", + ) + + with ( + uproot.open(scatt_root_filename) as f_sc, + uproot.open(abs_root_filename) as f_ab, + ): + t_sc = _get_tree(f_sc, scatt_tree_name) + t_ab = _get_tree(f_ab, abs_tree_name) + + branches_sc = set(t_sc.keys()) + branches_ab = set(t_ab.keys()) + + # check that all needed columns exist + miss_sc = set(required_branches) - branches_sc + miss_ab = set(required_branches) - branches_ab + if miss_sc: + raise ValueError(f"Scatterer singles missing branches: {sorted(miss_sc)}") + if miss_ab: + raise ValueError(f"Absorber singles missing branches: {sorted(miss_ab)}") + # check that both trees have the same branches + if branches_sc != branches_ab: + only_sc = branches_sc - branches_ab + only_ab = branches_ab - branches_sc + + if only_sc: + raise ValueError(f"Branches only in scatterer: {sorted(only_sc)}") + if only_ab: + raise ValueError(f"Branches only in absorber: {sorted(only_ab)}") + + # Load all branches from each tree into pandas + sc_df = t_sc.arrays(library="pd") + ab_df = t_ab.arrays(library="pd") + + merged = pd.concat([sc_df, ab_df], ignore_index=True) + # Optionally filter branches before writing + if save_branches is not None: + missing = set(save_branches) - set(merged.columns) + if missing: + raise ValueError( + f"Branches not found in merged data: {sorted(missing)}" + ) + + merged = merged[save_branches] + + # cc_coincidences_sorter expects time-ordered singles + merged = merged.sort_values( + ["EventID", "GlobalTime"], kind="mergesort" + ).reset_index(drop=True) + + # Write merged tree as a TTree + for col in merged.columns: + if merged[col].dtype == object: + merged[col] = pd.Categorical(merged[col]) + data = merged.to_dict(orient="list") + data = root_tree_get_branch_data(data) + types = root_tree_get_branch_types(data) + with uproot.recreate(output_root_filename) as f_out: + root_write_tree(f_out, output_tree_name, types, data) + + # Validate exact branch set + with uproot.open(output_root_filename) as f_chk: + if output_tree_name not in f_chk: + raise ValueError("Merged singles tree is empty; no data to write.") + out_keys = set(f_chk[output_tree_name].keys()) + + # Warn if "required branches" for sorting are not included in the merged file + missing_required = set(required_branches) - out_keys + if missing_required: + print( + "WARNING: The merged tree does not include all required branches for sorting.\n" + f"Missing required branches: {sorted(missing_required)}\n" + f"Required branches are: {sorted(required_branches)}" + ) diff --git a/opengate/bin/opengate_tests.py b/opengate/bin/opengate_tests.py index f34648f19..ab34e628c 100755 --- a/opengate/bin/opengate_tests.py +++ b/opengate/bin/opengate_tests.py @@ -64,6 +64,13 @@ default="", help="Only for developers: overwrite the used geant4 version str to pass the check, style: v11.4.0", ) +@click.option( + "--print_last_test", + "-l", + is_flag=True, + default=False, + help="Print the number of the last test and exit", +) def go( start_id, end_id, @@ -74,9 +81,24 @@ def go( run_previously_failed_jobs, num_processes, g4_version, + print_last_test, ): path_tests_src = return_tests_path() # returns the path to the tests/src dir + if print_last_test: + import re + + all_files = helpers.discover_all_tests(path_tests_src) + pattern = re.compile(r"^test([0-9]+)") + test_numbers = [] + for f in all_files: + match = pattern.match(os.path.basename(f)) + if match: + test_numbers.append(int(match.group(1))) + if test_numbers: + print(max(test_numbers)) + return + test_dir_path = path_tests_src.parent start = time.time() diff --git a/opengate/contrib/compton_camera/constants.yaml b/opengate/contrib/compton_camera/constants.yaml new file mode 100644 index 000000000..ab7c57642 --- /dev/null +++ b/opengate/contrib/compton_camera/constants.yaml @@ -0,0 +1,72 @@ +# SPDX-FileCopyrightText: 2024 Vincent Lequertier , Voichita Maxim +# +# SPDX-License-Identifier: MIT + +avogadro: 6.02214179e+23 +m_e: 510.99895000 # electron rest mass in keV +r_e: 2.8179403227e-13 # electron classical radius in cm + +# TODO: This depends on the material (Silicium currently). These values were +# computed by former students, Enrique or Yuemeng. Add new keys if we add other +# materials +doppler_broadening: + energy_140: + a1: 0.3773 + a2: 0.1443 + sigma_beta_1: 0.002090928993768552 # rad + sigma_beta_2: 0.01835770181680156 # rad + energy_245: + a1: 0.6291 + a2: 0.2335 + sigma_beta_1: 0.0014484966569412875 # rad + sigma_beta_2: 0.011749867543634142 # rad + energy_364: + a1: 0.8393 + a2: 0.2735 + sigma_beta_1: 0.0013006032426011388 # rad + sigma_beta_2: 0.009809405280466551 # rad + energy_511: + a1: 1.0000 + a2: 0.2980 + sigma_beta_1: 0.0012218929989592862 # rad + sigma_beta_2: 0.008513331113645892 # rad + energy_1275: + a1: 1.00000000 + a2: 0.21889611 + sigma_beta_1: 0.0010469077761151514 # rad + sigma_beta_2: 0.006182086979184491 # rad + energy_2170: + a1: 0.42537069 + a2: 0.33833700 + sigma_beta_1: 0.00040536484092147874 # rad + sigma_beta_2: 0.0035140981264683954 # rad + energy_4439: + a1: 0.60596747 + a2: 0.28269513 + sigma_beta_1: 7.811690920027629e-05 # rad + sigma_beta_2: 0.00273685099297266 # rad + energy_6130: + a1: 0.54893444 + a2: 0.23693385 + sigma_beta_1: 8.049907448898096e-05 # rad + sigma_beta_2: 0.002622310532968701 # rad + +materials: + Si: + density: 2.33 # g cm^-3 + moll_mass: 28.09 # g * cm^-3 + eff: 14 + # Comes from here https://www.physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Elem&Output2=Hand + NIST: [ [1.000E-03, 1.317E-02, 1.567E+03, 0.000E+00, 1.570E+03, 1.567E+03], [1.500E-03, 2.393E-02, 5.333E+02, 0.000E+00, 5.356E+02, 5.333E+02], [1.839E-03, 3.081E-02, 3.071E+02, 0.000E+00, 3.092E+02, 3.071E+02], [1.839E-03, 3.081E-02, 3.191E+03, 0.000E+00, 3.193E+03, 3.191E+03], [2.000E-03, 3.388E-02, 2.775E+03, 0.000E+00, 2.777E+03, 2.775E+03], [3.000E-03, 4.962E-02, 9.767E+02, 0.000E+00, 9.784E+02, 9.767E+02], [4.000E-03, 6.135E-02, 4.514E+02, 0.000E+00, 4.528E+02, 4.514E+02], [5.000E-03, 7.110E-02, 2.438E+02, 0.000E+00, 2.451E+02, 2.439E+02], [6.000E-03, 7.983E-02, 1.458E+02, 0.000E+00, 1.470E+02, 1.459E+02], [8.000E-03, 9.507E-02, 6.379E+01, 0.000E+00, 6.469E+01, 6.389E+01], [1.000E-02, 1.076E-01, 3.315E+01, 0.000E+00, 3.388E+01, 3.326E+01], [1.500E-02, 1.289E-01, 9.848E+00, 0.000E+00, 1.034E+01, 9.977E+00], [2.000E-02, 1.402E-01, 4.089E+00, 0.000E+00, 4.463E+00, 4.229E+00], [3.000E-02, 1.501E-01, 1.161E+00, 0.000E+00, 1.436E+00, 1.311E+00], [4.000E-02, 1.534E-01, 4.687E-01, 0.000E+00, 7.010E-01, 6.221E-01], [5.000E-02, 1.538E-01, 2.307E-01, 0.000E+00, 4.385E-01, 3.845E-01], [6.000E-02, 1.526E-01, 1.288E-01, 0.000E+00, 3.206E-01, 2.815E-01], [8.000E-02, 1.483E-01, 5.120E-02, 0.000E+00, 2.228E-01, 1.995E-01], [1.000E-01, 1.432E-01, 2.498E-02, 0.000E+00, 1.835E-01, 1.682E-01], [1.500E-01, 1.309E-01, 6.806E-03, 0.000E+00, 1.448E-01, 1.377E-01], [2.000E-01, 1.207E-01, 2.736E-03, 0.000E+00, 1.275E-01, 1.235E-01], [3.000E-01, 1.055E-01, 7.882E-04, 0.000E+00, 1.082E-01, 1.063E-01], [4.000E-01, 9.475E-02, 3.409E-04, 0.000E+00, 9.614E-02, 9.509E-02], [5.000E-01, 8.663E-02, 1.849E-04, 0.000E+00, 8.748E-02, 8.681E-02], [6.000E-01, 8.019E-02, 1.156E-04, 0.000E+00, 8.078E-02, 8.031E-02], [8.000E-01, 7.050E-02, 5.854E-05, 0.000E+00, 7.082E-02, 7.056E-02], [1.000E+00, 6.340E-02, 3.639E-05, 0.000E+00, 6.361E-02, 6.344E-02], [1.022E+00, 6.274E-02, 3.433E-05, 0.000E+00, 6.293E-02, 6.277E-02], [1.250E+00, 5.671E-02, 2.326E-05, 3.521E-05, 5.688E-02, 5.677E-02], [1.500E+00, 5.155E-02, 1.684E-05, 1.912E-04, 5.183E-02, 5.175E-02], [2.000E+00, 4.400E-02, 1.050E-05, 7.533E-04, 4.481E-02, 4.476E-02], [2.044E+00, 4.346E-02, 1.016E-05, 8.118E-04, 4.433E-02, 4.429E-02], [3.000E+00, 3.461E-02, 5.800E-06, 2.139E-03, 3.678E-02, 3.677E-02], [4.000E+00, 2.888E-02, 3.952E-06, 3.456E-03, 3.240E-02, 3.239E-02], [5.000E+00, 2.494E-02, 2.983E-06, 4.629E-03, 2.967E-02, 2.967E-02], [6.000E+00, 2.204E-02, 2.389E-06, 5.678E-03, 2.788E-02, 2.787E-02], [7.000E+00, 1.981E-02, 1.990E-06, 6.611E-03, 2.663E-02, 2.663E-02], [8.000E+00, 1.803E-02, 1.705E-06, 7.451E-03, 2.574E-02, 2.574E-02], [9.000E+00, 1.657E-02, 1.490E-06, 8.212E-03, 2.509E-02, 2.509E-02], [1.000E+01, 1.536E-02, 1.323E-06, 8.905E-03, 2.462E-02, 2.461E-02], [1.100E+01, 1.432E-02, 1.189E-06, 9.531E-03, 2.425E-02, 2.425E-02], [1.200E+01, 1.343E-02, 1.080E-06, 1.011E-02, 2.398E-02, 2.398E-02], [1.300E+01, 1.266E-02, 9.889E-07, 1.064E-02, 2.378E-02, 2.377E-02], [1.400E+01, 1.197E-02, 9.121E-07, 1.114E-02, 2.362E-02, 2.362E-02], [1.500E+01, 1.137E-02, 8.463E-07, 1.160E-02, 2.352E-02, 2.352E-02], [1.600E+01, 1.082E-02, 7.893E-07, 1.204E-02, 2.345E-02, 2.345E-02], [1.800E+01, 9.893E-03, 6.954E-07, 1.284E-02, 2.338E-02, 2.338E-02], [2.000E+01, 9.124E-03, 6.216E-07, 1.355E-02, 2.338E-02, 2.338E-02], [2.200E+01, 8.474E-03, 5.618E-07, 1.419E-02, 2.343E-02, 2.343E-02], [2.400E+01, 7.919E-03, 5.125E-07, 1.478E-02, 2.351E-02, 2.351E-02], [2.600E+01, 7.436E-03, 4.711E-07, 1.532E-02, 2.361E-02, 2.361E-02], [2.800E+01, 7.014E-03, 4.359E-07, 1.581E-02, 2.372E-02, 2.372E-02], [3.000E+01, 6.641E-03, 4.057E-07, 1.627E-02, 2.385E-02, 2.385E-02], [4.000E+01, 5.277E-03, 3.010E-07, 1.817E-02, 2.455E-02, 2.455E-02], [5.000E+01, 4.406E-03, 2.393E-07, 1.958E-02, 2.522E-02, 2.522E-02], [6.000E+01, 3.795E-03, 1.985E-07, 2.070E-02, 2.583E-02, 2.583E-02], [8.000E+01, 2.993E-03, 1.481E-07, 2.236E-02, 2.685E-02, 2.685E-02], [1.000E+02, 2.485E-03, 1.181E-07, 2.354E-02, 2.764E-02, 2.764E-02], [1.500E+02, 1.769E-03, 7.837E-08, 2.547E-02, 2.905E-02, 2.905E-02], [2.000E+02, 1.387E-03, 5.864E-08, 2.663E-02, 2.996E-02, 2.996E-02], [3.000E+02, 9.833E-04, 3.902E-08, 2.802E-02, 3.111E-02, 3.111E-02], [4.000E+02, 7.700E-04, 2.923E-08, 2.884E-02, 3.181E-02, 3.181E-02], [5.000E+02, 6.375E-04, 2.337E-08, 2.938E-02, 3.229E-02, 3.229E-02], [6.000E+02, 5.457E-04, 1.947E-08, 2.976E-02, 3.263E-02, 3.263E-02], [8.000E+02, 4.263E-04, 1.459E-08, 3.030E-02, 3.312E-02, 3.312E-02], [1.000E+03, 3.508E-04, 1.167E-08, 3.064E-02, 3.344E-02, 3.344E-02]] + BGO: + density: 7.13 # g cm^-3 + moll_mass: 1245.83 # g * cm^-3 + eff: 524 + # Comes from here https://www.physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Elem&Output2=Hand + NIST: [[1.000E-03, 4.741E-03, 4.680E+03, 0.000E+00, 4.690E+03, 4.680E+03], [1.103E-03, 5.573E-03, 3.813E+03, 0.000E+00, 3.823E+03, 3.813E+03], [1.217E-03, 6.540E-03, 3.109E+03, 0.000E+00, 3.118E+03, 3.109E+03], [1.232E-03, 6.674E-03, 3.640E+03, 0.000E+00, 3.650E+03, 3.640E+03], [1.248E-03, 6.811E-03, 3.623E+03, 0.000E+00, 3.632E+03, 3.623E+03], [1.328E-03, 7.526E-03, 3.483E+03, 0.000E+00, 3.493E+03, 3.483E+03], [1.414E-03, 8.301E-03, 3.099E+03, 0.000E+00, 3.108E+03, 3.099E+03], [1.500E-03, 9.082E-03, 2.843E+03, 0.000E+00, 2.852E+03, 2.843E+03], [2.000E-03, 1.370E-02, 1.477E+03, 0.000E+00, 1.485E+03, 1.477E+03], [2.580E-03, 1.901E-02, 8.092E+02, 0.000E+00, 8.174E+02, 8.092E+02], [2.633E-03, 1.949E-02, 1.496E+03, 0.000E+00, 1.504E+03, 1.496E+03], [2.688E-03, 1.998E-02, 1.506E+03, 0.000E+00, 1.515E+03, 1.506E+03], [3.000E-03, 2.271E-02, 1.571E+03, 0.000E+00, 1.579E+03, 1.571E+03], [3.177E-03, 2.420E-02, 1.355E+03, 0.000E+00, 1.363E+03, 1.355E+03], [3.427E-03, 2.629E-02, 1.275E+03, 0.000E+00, 1.283E+03, 1.275E+03], [3.696E-03, 2.842E-02, 1.057E+03, 0.000E+00, 1.064E+03, 1.057E+03], [3.845E-03, 2.953E-02, 1.012E+03, 0.000E+00, 1.019E+03, 1.012E+03], [3.999E-03, 3.075E-02, 9.201E+02, 0.000E+00, 9.271E+02, 9.201E+02], [4.000E-03, 3.076E-02, 9.558E+02, 0.000E+00, 9.628E+02, 9.558E+02], [5.000E-03, 3.781E-02, 5.529E+02, 0.000E+00, 5.592E+02, 5.529E+02], [6.000E-03, 4.402E-02, 3.509E+02, 0.000E+00, 3.565E+02, 3.509E+02], [8.000E-03, 5.448E-02, 1.689E+02, 0.000E+00, 1.734E+02, 1.689E+02], [1.000E-02, 6.292E-02, 9.493E+01, 0.000E+00, 9.872E+01, 9.499E+01], [1.110E-02, 6.687E-02, 7.222E+01, 0.000E+00, 7.566E+01, 7.228E+01], [1.221E-02, 7.042E-02, 7.982E+01, 0.000E+00, 8.296E+01, 7.989E+01], [1.342E-02, 7.394E-02, 6.245E+01, 0.000E+00, 6.530E+01, 6.252E+01], [1.500E-02, 7.805E-02, 9.156E+01, 0.000E+00, 9.409E+01, 9.163E+01], [1.571E-02, 7.972E-02, 8.089E+01, 0.000E+00, 8.330E+01, 8.097E+01], [1.605E-02, 8.047E-02, 1.017E+02, 0.000E+00, 1.040E+02, 1.018E+02], [1.639E-02, 8.122E-02, 9.662E+01, 0.000E+00, 9.892E+01, 9.670E+01], [2.000E-02, 8.804E-02, 6.576E+01, 0.000E+00, 6.758E+01, 6.584E+01], [3.000E-02, 1.004E-01, 2.251E+01, 0.000E+00, 2.363E+01, 2.261E+01], [4.000E-02, 1.070E-01, 1.037E+01, 0.000E+00, 1.116E+01, 1.048E+01], [5.000E-02, 1.103E-01, 5.646E+00, 0.000E+00, 6.238E+00, 5.756E+00], [6.000E-02, 1.117E-01, 3.422E+00, 0.000E+00, 3.894E+00, 3.534E+00], [8.000E-02, 1.118E-01, 1.546E+00, 0.000E+00, 1.884E+00, 1.658E+00], [9.053E-02, 1.111E-01, 1.098E+00, 0.000E+00, 1.393E+00, 1.209E+00], [9.053E-02, 1.111E-01, 4.804E+00, 0.000E+00, 5.099E+00, 4.915E+00], [1.000E-01, 1.102E-01, 3.705E+00, 0.000E+00, 3.971E+00, 3.815E+00], [1.500E-01, 1.039E-01, 1.281E+00, 0.000E+00, 1.462E+00, 1.385E+00], [2.000E-01, 9.742E-02, 5.983E-01, 0.000E+00, 7.415E-01, 6.957E-01], [3.000E-01, 8.659E-02, 2.073E-01, 0.000E+00, 3.157E-01, 2.939E-01], [4.000E-01, 7.842E-02, 1.004E-01, 0.000E+00, 1.915E-01, 1.788E-01], [5.000E-01, 7.206E-02, 5.852E-02, 0.000E+00, 1.389E-01, 1.306E-01], [6.000E-01, 6.693E-02, 3.833E-02, 0.000E+00, 1.111E-01, 1.053E-01], [8.000E-01, 5.906E-02, 2.037E-02, 0.000E+00, 8.280E-02, 7.943E-02], [1.000E+00, 5.322E-02, 1.284E-02, 0.000E+00, 6.824E-02, 6.606E-02], [1.022E+00, 5.267E-02, 1.229E-02, 0.000E+00, 6.704E-02, 6.496E-02], [1.250E+00, 4.767E-02, 8.293E-03, 2.784E-04, 5.765E-02, 5.624E-02], [1.500E+00, 4.337E-02, 5.907E-03, 1.339E-03, 5.160E-02, 5.062E-02], [2.000E+00, 3.706E-02, 3.572E-03, 4.095E-03, 4.528E-02, 4.472E-02], [2.044E+00, 3.660E-02, 3.443E-03, 4.340E-03, 4.491E-02, 4.438E-02], [3.000E+00, 2.918E-02, 1.868E-03, 9.118E-03, 4.042E-02, 4.017E-02], [4.000E+00, 2.435E-02, 1.223E-03, 1.321E-02, 3.897E-02, 3.883E-02], [5.000E+00, 2.103E-02, 8.960E-04, 1.667E-02, 3.877E-02, 3.868E-02], [6.000E+00, 1.859E-02, 7.020E-04, 1.966E-02, 3.914E-02, 3.908E-02], [7.000E+00, 1.671E-02, 5.748E-04, 2.231E-02, 3.981E-02, 3.977E-02], [8.000E+00, 1.521E-02, 4.855E-04, 2.468E-02, 4.062E-02, 4.059E-02], [9.000E+00, 1.398E-02, 4.196E-04, 2.684E-02, 4.153E-02, 4.150E-02], [1.000E+01, 1.295E-02, 3.692E-04, 2.882E-02, 4.246E-02, 4.244E-02], [1.100E+01, 1.208E-02, 3.291E-04, 3.066E-02, 4.342E-02, 4.340E-02], [1.200E+01, 1.133E-02, 2.969E-04, 3.235E-02, 4.436E-02, 4.435E-02], [1.300E+01, 1.068E-02, 2.703E-04, 3.395E-02, 4.531E-02, 4.530E-02], [1.400E+01, 1.010E-02, 2.480E-04, 3.544E-02, 4.623E-02, 4.622E-02], [1.500E+01, 9.590E-03, 2.291E-04, 3.682E-02, 4.710E-02, 4.709E-02], [1.600E+01, 9.133E-03, 2.128E-04, 3.813E-02, 4.796E-02, 4.796E-02]] + LaBr3: + density: 5.08 # g cm^-3 + moll_mass: 378.62 # g * cm^-3 + eff: 162 + # Comes from here https://www.physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Elem&Output2=Hand + NIST: [[1.00E-003, 5.80E-003, 4.99E+003, 0.00E+000, 5.00E+003, 4.988E+03], [1.06E-003, 6.34E-003, 4.37E+003, 0.00E+000, 4.37E+003, 4.366E+03], [1.12E-003, 6.92E-003, 3.82E+003, 0.00E+000, 3.83E+003, 3.822E+03], [1.12E-003, 6.92E-003, 4.20E+003, 0.00E+000, 4.20E+003, 4.197E+03], [1.16E-003, 7.29E-003, 3.89E+003, 0.00E+000, 3.89E+003, 3.885E+03], [1.20E-003, 7.67E-003, 3.60E+003, 0.00E+000, 3.60E+003, 3.596E+03], [1.20E-003, 7.67E-003, 3.75E+003, 0.00E+000, 3.75E+003, 3.748E+03], [1.28E-003, 8.38E-003, 3.28E+003, 0.00E+000, 3.28E+003, 3.276E+03], [1.36E-003, 9.14E-003, 2.87E+003, 0.00E+000, 2.87E+003, 2.866E+03], [1.36E-003, 9.14E-003, 2.96E+003, 0.00E+000, 2.97E+003, 2.958E+03], [1.50E-003, 1.05E-002, 2.38E+003, 0.00E+000, 2.39E+003, 2.379E+03], [1.55E-003, 1.09E-002, 2.20E+003, 0.00E+000, 2.21E+003, 2.204E+03], [1.55E-003, 1.09E-002, 4.35E+003, 0.00E+000, 4.35E+003, 4.348E+03], [1.57E-003, 1.12E-002, 4.06E+003, 0.00E+000, 4.06E+003, 4.056E+03], [1.60E-003, 1.14E-002, 3.79E+003, 0.00E+000, 3.79E+003, 3.787E+03], [1.60E-003, 1.14E-002, 4.74E+003, 0.00E+000, 4.75E+003, 4.742E+03], [1.69E-003, 1.23E-002, 4.18E+003, 0.00E+000, 4.19E+003, 4.181E+03], [1.78E-003, 1.32E-002, 3.69E+003, 0.00E+000, 3.69E+003, 3.685E+03], [1.78E-003, 1.32E-002, 4.02E+003, 0.00E+000, 4.02E+003, 4.018E+03], [2.00E-003, 1.52E-002, 3.06E+003, 0.00E+000, 3.06E+003, 3.055E+03], [3.00E-003, 2.43E-002, 1.11E+003, 0.00E+000, 1.12E+003, 1.114E+03], [4.00E-003, 3.24E-002, 5.30E+002, 0.00E+000, 5.34E+002, 5.297E+02], [5.00E-003, 3.95E-002, 2.94E+002, 0.00E+000, 2.98E+002, 2.940E+02], [5.48E-003, 4.26E-002, 2.30E+002, 0.00E+000, 2.34E+002, 2.298E+02], [5.48E-003, 4.26E-002, 3.66E+002, 0.00E+000, 3.70E+002, 3.664E+02], [5.68E-003, 4.39E-002, 3.34E+002, 0.00E+000, 3.38E+002, 3.344E+02], [5.89E-003, 4.51E-002, 3.05E+002, 0.00E+000, 3.09E+002, 3.051E+02], [5.89E-003, 4.51E-002, 3.68E+002, 0.00E+000, 3.71E+002, 3.676E+02], [6.00E-003, 4.57E-002, 3.53E+002, 0.00E+000, 3.57E+002, 3.529E+02], [6.27E-003, 4.73E-002, 3.15E+002, 0.00E+000, 3.18E+002, 3.148E+02], [6.27E-003, 4.73E-002, 3.46E+002, 0.00E+000, 3.50E+002, 3.464E+02], [8.00E-003, 5.63E-002, 1.84E+002, 0.00E+000, 1.87E+002, 1.838E+02], [1.00E-002, 6.52E-002, 1.01E+002, 0.00E+000, 1.03E+002, 1.010E+02], [1.35E-002, 7.78E-002, 4.49E+001, 0.00E+000, 4.66E+001, 4.493E+01], [1.35E-002, 7.78E-002, 1.24E+002, 0.00E+000, 1.26E+002, 1.243E+02], [1.50E-002, 8.23E-002, 9.40E+001, 0.00E+000, 9.55E+001, 9.408E+01], [2.00E-002, 9.38E-002, 4.37E+001, 0.00E+000, 4.48E+001, 4.378E+01], [3.00E-002, 1.07E-001, 1.43E+001, 0.00E+000, 1.50E+001, 1.439E+01], [3.89E-002, 1.13E-001, 6.83E+000, 0.00E+000, 7.33E+000, 6.943E+00], [3.89E-002, 1.13E-001, 1.51E+001, 0.00E+000, 1.56E+001, 1.518E+01], [4.00E-002, 1.13E-001, 1.40E+001, 0.00E+000, 1.45E+001, 1.409E+01], [5.00E-002, 1.17E-001, 7.63E+000, 0.00E+000, 8.01E+000, 7.750E+00], [6.00E-002, 1.18E-001, 4.61E+000, 0.00E+000, 4.92E+000, 4.730E+00], [8.00E-002, 1.17E-001, 2.06E+000, 0.00E+000, 2.29E+000, 2.173E+00], [1.00E-001, 1.15E-001, 1.09E+000, 0.00E+000, 1.28E+000, 1.204E+00], [1.50E-001, 1.08E-001, 3.39E-001, 0.00E+000, 4.86E-001, 4.471E-01], [2.00E-001, 1.01E-001, 1.49E-001, 0.00E+000, 2.72E-001, 2.491E-01], [3.00E-001, 8.90E-002, 4.74E-002, 0.00E+000, 1.47E-001, 1.364E-01], [4.00E-001, 8.04E-002, 2.18E-002, 0.00E+000, 1.08E-001, 1.022E-01], [5.00E-001, 7.37E-002, 1.23E-002, 0.00E+000, 9.00E-002, 8.602E-02], [6.00E-001, 6.84E-002, 7.86E-003, 0.00E+000, 7.90E-002, 7.626E-02], [8.00E-001, 6.03E-002, 4.07E-003, 0.00E+000, 6.59E-002, 6.433E-02], [1.00E+000, 5.43E-002, 2.54E-003, 0.00E+000, 5.78E-002, 5.682E-02], [1.02E+000, 5.37E-002, 2.43E-003, 0.00E+000, 5.71E-002, 5.612E-02], [1.25E+000, 4.86E-002, 1.64E-003, 1.41E-004, 5.10E-002, 5.035E-02], [1.50E+000, 4.42E-002, 1.17E-003, 6.68E-004, 4.65E-002, 4.602E-02], [2.00E+000, 3.77E-002, 7.14E-004, 2.32E-003, 4.10E-002, 4.076E-02], [2.04E+000, 3.73E-002, 6.89E-004, 2.48E-003, 4.07E-002, 4.043E-02], [3.00E+000, 2.97E-002, 3.80E-004, 6.00E-003, 3.62E-002, 3.609E-02], [4.00E+000, 2.48E-002, 2.52E-004, 9.31E-003, 3.44E-002, 3.437E-02], [5.00E+000, 2.14E-002, 1.86E-004, 1.22E-002, 3.39E-002, 3.385E-02], [6.00E+000, 1.89E-002, 1.47E-004, 1.47E-002, 3.39E-002, 3.386E-02], [7.00E+000, 1.70E-002, 1.21E-004, 1.69E-002, 3.42E-002, 3.420E-02], [8.00E+000, 1.55E-002, 1.03E-004, 1.89E-002, 3.47E-002, 3.469E-02], [9.00E+000, 1.42E-002, 8.92E-005, 2.07E-002, 3.53E-002, 3.528E-02], [1.00E+001, 1.32E-002, 7.87E-005, 2.24E-002, 3.59E-002, 3.593E-02]] diff --git a/opengate/contrib/compton_camera/coresi_default_config.yaml b/opengate/contrib/compton_camera/coresi_default_config.yaml new file mode 100644 index 000000000..5bce2c805 --- /dev/null +++ b/opengate/contrib/compton_camera/coresi_default_config.yaml @@ -0,0 +1,141 @@ +# SPDX-FileCopyrightText: 2024 Vincent Lequertier , Voichita Maxim +# +# SPDX-License-Identifier: MIT + +# data_file: CC_CLARYS_7P_IECPhantom_Cylinders_4Energies_t5000s_Ideal_140keVPeak.dat +data_file: simulation.raw +data_type: 'GATE' +log_dir: 'logs' +# Energy of the source. If set to -1 the gamma energy is the deposited energy i.e. E1 + E2 +# E0: [140, 245, 364, 511] +E0: [511.0] +remove_out_of_range_energies: False +# Only needed if E0 is unknown and remove_out_of_range_energies is True, ignored +# otherwise +energy_range: [100, 600] +# Ethres = 5 pour données idéales, Ethres = 50 pour resolution finie des d +# tecteurs. +# (valeurs Enrique, mais devrait être choisi en fonction de la résolution +# des détecteurs) +energy_threshold: 5 +n_events: 10 +starts_at: 0 + + +# Cameras have scatterer layers normal on z axis, centered on the z axis, below the source +# Up to 5 absorbers are implemented +cameras: + n_cameras: 6 + # Cameras scatterer have the same characteristics besides their positions + common_attributes: + n_sca_layers: 7 + sca_material: Si + abs_material: BGO + sca_layer_0: + center: [0, 0, -10] + size: [9.0, 9.0, 0.2] + sca_layer_1: + center: [0, 0, -11] + size: [9.0, 9.0, 0.2] + sca_layer_2: + center: [0, 0, -12] + size: [9.0, 9.0, 0.2] + sca_layer_3: + center: [0, 0, -13] + size: [9.0, 9.0, 0.2] + sca_layer_4: + center: [0, 0, -14] + size: [9.0, 9.0, 0.2] + sca_layer_5: + center: [0, 0, -15] + size: [9.0, 9.0, 0.2] + sca_layer_6: + center: [0, 0, -16] + size: [9.0, 9.0, 0.2] + + n_absorbers: 1 + + # if an absorber layer is missing, define it as None + # For example, abs_layer_1: None + abs_layer_0: # bottom + center: [0, 0, -31.0] + size: [28.0, 21.0, 3] + abs_layer_1: None # front + abs_layer_2: None # rear + abs_layer_3: None # right + abs_layer_4: None # left + + position_0: + frame_origin: [0, 0, 0] + Ox: [1, 0, 0] # parallel to scatterer edge + Oy: [0, 1, 0] # parallel to scatterer edge + Oz: [0, 0, 1] # orthogonal to the camera, tw the source + + position_1: # Second camera + frame_origin: [0, 0, 0] + Ox: [0, 0, 1] # parallel to scatterer edge + Oy: [0, 1, 0] # parallel to scatterer edge + Oz: [-1, 0, 0] # orthogonal to the camera, tw the source + + position_2: # Third camera + frame_origin: [0, 0, 0] + Ox: [0, 0, -1] # parallel to scatterer edge + Oy: [0, 1, 0] # parallel to scatterer edge + Oz: [1, 0, 0] # orthogonal to the camera, tw the source + + position_3: # Fourth camera + frame_origin: [0, 0, 0] + Ox: [-1, 0, 0] # parallel to scatterer edge + Oy: [0, 1, 0] # parallel to scatterer edge + Oz: [0, 0, -1] # orthogonal to the camera, tw the source + + position_4: # Fifth camera + frame_origin: [0, 0, 0] + Ox: [0.707106781187, 0, -0.707106781187] # parallel to scatterer edge + Oy: [0, 1, 0] # parallel to scatterer edge + Oz: [0.707106781187, 0, 0.707106781187] # orthogonal to the camera, tw the source + + position_5: # Sixth camera + frame_origin: [0, 0, 0] + Ox: [0.707106781187, 0, 0.707106781187] # parallel to scatterer edge + Oy: [0, 1, 0] # parallel to scatterer edge + Oz: [-0.707106781187, 0, 0.707106781187] # orthogonal to the camera, tw the source + + +volume: + volume_dimensions: [20.25, 20.25, 10.25] # In centimeters + n_voxels: [81, 81, 41] + # volume_dimensions: [20.25, 20.25, 0.25] # In centimeters + # n_voxels: [81,81,1] + # volume_dimensions: [20,20,0.4] # In centimeters + # n_voxels: [50,50,1] + volume_centre: [0, 0, 0] + +simulation: + n_events: 20000 + # n_events: 20 + n_V2: 100 # samples of V2 positions + phantom: misc/phantom.pth + output_file: simulation.raw + visualize_generated_source: True + angle_threshold: 5 # In degrees. It's hard to compute the exact V2 (second hit in the absorber) that corresponds to the randomly chosen angle. So draw candidates and pick one that loosely match the angle according to this threshold + +lm_mlem: + cone_thickness: 'parallel' # angular or parallel + last_iter: 3 + model: 'cos0rho0' # Either cos0rho0, cos0rho2 or cos1rho2. Use cos0rho0 if sensitivity is disabled. + first_iter: 0 + save_every: 1 + checkpoint_dir: checkpoints_test + # force sperctral algorithm even if there is only one energy + force_spectral: false + n_sigma: 2 # Skip the Gaussian above n_sigma * Gaussian std + width_factor: 1 # Used to artificially modify the spread of the cone + sensitivity: false + sensitivity_model: like_system_matrix #from less precise to more precise: solid_angle, solid_angle_with_attn, like_system_matrix + mc_samples: 50000 + tv: false + alpha_tv: 0.000 + # Supply a sensitivity file. Needs to have the same size as the number of + # voxels in the volume + sensitivity_file: none diff --git a/opengate/contrib/compton_camera/coresi_helpers.py b/opengate/contrib/compton_camera/coresi_helpers.py new file mode 100644 index 000000000..4f12e54b6 --- /dev/null +++ b/opengate/contrib/compton_camera/coresi_helpers.py @@ -0,0 +1,192 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import yaml +import uproot +from pathlib import Path +import opengate_core as g4 +from opengate.exception import fatal +from opengate.utility import g4_units +from opengate.geometry.utility import vec_g4_as_np, rot_g4_as_np + + +# --- Custom List for Inline YAML Formatting --- +class FlowList(list): + """A custom list that will be dumped as [x, y, z] in YAML.""" + + pass + + +def flow_list_representer(dumper, data): + return dumper.represent_sequence("tag:yaml.org,2002:seq", data, flow_style=True) + + +yaml.add_representer(FlowList, flow_list_representer) + + +def convert_to_flowlist(data): + """ + Recursively traverse the dictionary. + Convert any list containing only simple scalars (int, float, str) to FlowList + so they appear as [a, b, c] in the YAML output. + """ + if isinstance(data, dict): + return {k: convert_to_flowlist(v) for k, v in data.items()} + elif isinstance(data, list): + # Check if list is "simple" (contains only primitives, no dicts or nested lists) + is_simple = all( + isinstance(i, (int, float, str, bool)) or i is None for i in data + ) + + if is_simple: + return FlowList(data) + else: + # If list contains complex objects, process them recursively but keep the list as is + return [convert_to_flowlist(i) for i in data] + else: + return data + + +def coresi_new_config(): + # get current path of the script + current_dir = Path(__file__).parent.resolve() + print(current_dir) + # read the yaml file + with open(current_dir / "coresi_default_config.yaml", "r") as f: + config = yaml.safe_load(f) + return config + + +def set_hook_coresi_config(sim, cameras, filename): + """ + Prepare everything to create the coresi config file at the init of the simulation. + The param structure allows retrieving the coresi config at the end of the simulation. + """ + # create the param structure + param = { + "cameras": cameras, + "filename": filename, + "coresi_config": coresi_new_config(), + } + sim.user_hook_after_init = create_coresi_config + sim.user_hook_after_init_arg = param + return param + + +def create_coresi_config(simulation_engine, param): + # (note: simulation_engine is not used here but must be the first param) + coresi_config = param["coresi_config"] + cameras = param["cameras"] + + for camera in cameras.values(): + c = coresi_config["cameras"] + c["n_cameras"] += 1 + scatter_layer_names = camera["scatter_layer_names"] + absorber_layer_names = camera["absorber_layer_names"] + + for layer_name in scatter_layer_names: + coresi_add_scatterer(coresi_config, layer_name) + for layer_name in absorber_layer_names: + coresi_add_absorber(coresi_config, layer_name) + + +def coresi_add_scatterer(coresi_config, layer_name): + # find all volumes ('touchable' in Geant4 terminology) + touchables = g4.FindAllTouchables(layer_name) + if len(touchables) != 1: + fatal(f"Cannot find unique volume for layer {layer_name}: {touchables}") + touchable = touchables[0] + + # current nb of scatterers + id = coresi_config["cameras"]["common_attributes"]["n_sca_layers"] + coresi_config["cameras"]["common_attributes"]["n_sca_layers"] += 1 + layer = { + "center": [0, 0, 0], + "size": [0, 0, 0], + } + coresi_config["cameras"]["common_attributes"][f"sca_layer_{id}"] = layer + + # Get the information: WARNING in cm! + cm = g4_units.cm + translation = vec_g4_as_np(touchable.GetTranslation(0)) / cm + solid = touchable.GetSolid(0) + pMin_local = g4.G4ThreeVector() + pMax_local = g4.G4ThreeVector() + solid.BoundingLimits(pMin_local, pMax_local) + size = [ + (pMax_local.x - pMin_local.x) / cm, + (pMax_local.y - pMin_local.y) / cm, + (pMax_local.z - pMin_local.z) / cm, + ] + layer["center"] = translation.tolist() + layer["size"] = size + + +def coresi_add_absorber(coresi_config, layer_name): + # find all volumes ('touchable' in Geant4 terminology) + touchables = g4.FindAllTouchables(layer_name) + if len(touchables) != 1: + fatal(f"Cannot find unique volume for layer {layer_name}: {touchables}") + touchable = touchables[0] + + # current nb of scatterers + id = coresi_config["cameras"]["common_attributes"]["n_absorbers"] + coresi_config["cameras"]["common_attributes"]["n_absorbers"] += 1 + layer = { + "center": [0, 0, 0], + "size": [0, 0, 0], + } + coresi_config["cameras"]["common_attributes"][f"abs_layer_{id}"] = layer + + # Get the information: WARNING in cm! + cm = g4_units.cm + translation = vec_g4_as_np(touchable.GetTranslation(0)) / cm + solid = touchable.GetSolid(0) + pMin_local = g4.G4ThreeVector() + pMax_local = g4.G4ThreeVector() + solid.BoundingLimits(pMin_local, pMax_local) + size = [ + (pMax_local.x - pMin_local.x) / cm, + (pMax_local.y - pMin_local.y) / cm, + (pMax_local.z - pMin_local.z) / cm, + ] + layer["center"] = translation.tolist() + layer["size"] = size + + +def coresi_write_config(coresi_config, filename): + # Convert vectors to FlowList just before writing + formatted_config = convert_to_flowlist(coresi_config) + + with open(filename, "w") as f: + yaml.dump( + formatted_config, f, default_flow_style=False, sort_keys=False, indent=2 + ) + + +def coresi_convert_root_data(root_filename, branch_name, output_filename): + root_file = uproot.open(root_filename) + tree = root_file[branch_name] + # Load all needed branches at once + arrays = tree.arrays( + ["X1", "Y1", "Z1", "Energy1", "X2", "Y2", "Z2", "EnergyRest"], library="np" + ) + + with open(output_filename, "w") as fout: + n = len(arrays["X1"]) + for i in range(n): + # Energy in keV and position in mm, as expected by coresi + line = ( + f"2\t1\t" + f"{arrays['X1'][i]:.2f}\t" + f"{arrays['Y1'][i]:.2f}\t" + f"{arrays['Z1'][i]:.2f}\t" + f"{arrays['Energy1'][i]*1000:.2f}\t" + f"2\t" + f"{arrays['X2'][i]:.2f}\t" + f"{arrays['Y2'][i]:.2f}\t" + f"{arrays['Z2'][i]:.2f}\t" + f"{arrays['EnergyRest'][i]*1000:.2f}\t" + f"3\t0\t0\t0\t0\n" + ) + fout.write(line) diff --git a/opengate/contrib/compton_camera/macaco.py b/opengate/contrib/compton_camera/macaco.py new file mode 100644 index 000000000..44b29d1ae --- /dev/null +++ b/opengate/contrib/compton_camera/macaco.py @@ -0,0 +1,431 @@ +from pathlib import Path +import uproot +from opengate.utility import g4_units +from opengate.actors.coincidences import * +import numpy as np +import matplotlib.pyplot as plt + + +def add_macaco1_materials(sim): + """ + Adds the Macaco materials database to the simulation if not already present. + """ + db_folder = Path(__file__).parent.resolve() + db_filename = db_folder / "macaco_materials.db" + if not db_filename in sim.volume_manager.material_database.filenames: + sim.volume_manager.add_material_database(db_filename) + + +def add_macaco1_camera(sim, name="macaco1"): + """ + Adds a MACACO1 camera to the simulation. + - Bounding box (BB_box) + - Plane 1 (scatterer side): holder, crystal carcass, PCB, SiPM, LaBr3_VLC crystal + - Plane 2 (absorber side): holder, crystal carcass, PCB, SiPM, LaBr3_VLC crystal + """ + + # units + mm = g4_units.mm + cm = g4_units.cm + + # material + add_macaco1_materials(sim) + + # BB box (acts as camera envelope) + camera = sim.add_volume("Box", f"{name}_BB_box") + camera.mother = sim.world + camera.material = "G4_AIR" + camera.size = [16 * cm, 40 * cm, 7.6 * cm] + camera.color = [0.1, 0.1, 0.1, 0.1] + + # Scatterer + holder1 = sim.add_volume("Box", f"{name}_Holder1") + holder1.mother = camera.name + holder1.material = "Plastic" + holder1.size = [6.2 * cm, 6.2 * cm, 0.56 * cm] + holder1.translation = [0, 0, -2.74 * cm] + holder1.color = [0.1, 0.1, 0.1, 0.9] + + crys_carcase_scatt = sim.add_volume("Box", f"{name}_crysCarcaseScatt") + crys_carcase_scatt.mother = holder1.name + crys_carcase_scatt.material = "G4_Al" + crys_carcase_scatt.size = [2.72 * cm, 2.68 * cm, 0.01 * cm] + crys_carcase_scatt.translation = [0, 0, -0.265 * cm] + + pcb_scatt = sim.add_volume("Box", f"{name}_PCBScatt") + pcb_scatt.mother = camera.name + pcb_scatt.material = "PCB" + pcb_scatt.size = [10.89 * cm, 20.7 * cm, 0.4 * cm] + # pcb_scatt.translation = [0, 6.25 * cm, -2.46 * cm] + pcb_scatt.translation = [0, 6.25 * cm, -2.26 * cm] + pcb_scatt.color = [0.0, 0.5, 0.0, 0.9] + + sipm_scatt = sim.add_volume("Box", f"{name}_SiPMScatt") + sipm_scatt.mother = holder1.name + sipm_scatt.material = "G4_Si" + sipm_scatt.size = [2.72 * cm, 2.68 * cm, 0.04 * cm] + sipm_scatt.translation = [0, 0, 0.26 * cm] + sipm_scatt.color = [1.0, 0.5, 0.0, 0.9] + + scatterer = sim.add_volume("Box", f"{name}_scatterer") + scatterer.mother = holder1.name + scatterer.material = "LaBr3_VLC" + scatterer.size = [2.72 * cm, 2.68 * cm, 0.50 * cm] + scatterer.translation = [0, 0, -0.01 * cm] + scatterer.color = [0.4, 0.7, 1.0, 1.0] + + # Absorber + holder2 = sim.add_volume("Box", f"{name}_Holder2") + holder2.mother = camera.name + holder2.material = "Plastic" + holder2.size = [8.0 * cm, 8.0 * cm, 1.06 * cm] + holder2.translation = [0, 0, 2.51 * cm] + holder2.color = [0.1, 0.1, 0.1, 0.9] + + crys_carcase_abs = sim.add_volume("Box", f"{name}_crysCarcaseAbs") + crys_carcase_abs.mother = holder2.name + crys_carcase_abs.material = "G4_Al" + crys_carcase_abs.size = [3.24 * cm, 3.60 * cm, 0.01 * cm] + crys_carcase_abs.translation = [0, 0, -0.515 * cm] + + absorber = sim.add_volume("Box", f"{name}_absorber") + absorber.mother = holder2.name + absorber.material = "LaBr3_VLC" + absorber.size = [3.24 * cm, 3.60 * cm, 1.00 * cm] + absorber.translation = [0, 0, -0.01 * cm] + absorber.color = [0.4, 0.7, 1.0, 1.0] + + sipm_abs = sim.add_volume("Box", f"{name}_SiPMAbs") + sipm_abs.mother = holder2.name + sipm_abs.material = "G4_Si" + sipm_abs.size = [3.24 * cm, 3.60 * cm, 0.04 * cm] + sipm_abs.translation = [0, 0, 0.51 * cm] + sipm_abs.color = [1.0, 0.5, 0.0, 0.9] + + pcb_abs = sim.add_volume("Box", f"{name}_PCBAbs") + pcb_abs.mother = camera.name + pcb_abs.material = "PCB" + pcb_abs.size = [9.54 * cm, 16.0 * cm, 0.4 * cm] + pcb_abs.translation = [0, 4.50 * cm, 3.24 * cm] + pcb_abs.color = [0.0, 0.5, 0.0, 0.9] + + return { + "camera": camera, + "scatterer": scatterer, + "absorber": absorber, + "holder1": holder1, + "holder2": holder2, + } + + +def add_macaco1_camera_digitizer(sim, scatterer, absorber): + # Units + keV = g4_units.keV + MeV = g4_units.MeV + mm = g4_units.mm + ns = g4_units.ns + + # 1) Collect step-level hits for both scatterer and absorber + hits_scatt = sim.add_actor("DigitizerHitsCollectionActor", "HitsScatt") + hits_scatt.attached_to = scatterer.name + hits_scatt.attributes = [ + "EventID", + "TrackID", + "TotalEnergyDeposit", + "GlobalTime", + "PrePosition", + "PostPosition", + "PreStepUniqueVolumeID", + ] + + hits_abs = sim.add_actor("DigitizerHitsCollectionActor", "HitsAbs") + hits_abs.attached_to = absorber.name + hits_abs.attributes = hits_scatt.attributes + scatt_collection = hits_scatt.name + abs_collection = hits_abs.name + + # 2) Process Hits into Singles + sing_scatt = sim.add_actor("DigitizerAdderActor", "SinglesScatt") + sing_scatt.input_digi_collection = scatt_collection + sing_scatt.policy = "EnergyWeightedCentroidPosition" + sing_scatt.attributes = [ + "EventID", + "GlobalTime", + "TotalEnergyDeposit", + "PostPosition", + ] + scatt_collection = sing_scatt.name + + sing_abs = sim.add_actor("DigitizerAdderActor", "SinglesAbs") + sing_abs.input_digi_collection = abs_collection + sing_abs.policy = "EnergyWeightedCentroidPosition" + sing_abs.attributes = sing_scatt.attributes + abs_collection = sing_abs.name + + # Spatial blurring + spat_scatt = sim.add_actor("DigitizerSpatialBlurringActor", "SpatScatt") + spat_scatt.attached_to = scatterer.name + spat_scatt.input_digi_collection = scatt_collection + spat_scatt.blur_attribute = "PostPosition" + spat_scatt.blur_fwhm = [2 * mm, 2 * mm, 2 * mm] + scatt_collection = spat_scatt.name + + spat_abs = sim.add_actor("DigitizerSpatialBlurringActor", "SpatAbs") + spat_abs.attached_to = absorber.name + spat_abs.input_digi_collection = abs_collection + spat_abs.blur_attribute = "PostPosition" + spat_abs.blur_fwhm = [2 * mm, 2 * mm, 2 * mm] + abs_collection = spat_abs.name + + # Energy blurring + reference_energy = 511 * keV + scatt_resolution = 0.085 # FWHM/E at 511 keV + abs_resolution = 0.125 # FWHM/E at 511 keV + + blur_scatt = sim.add_actor("DigitizerBlurringActor", "BlurScatt") + blur_scatt.attached_to = scatterer.name + blur_scatt.input_digi_collection = scatt_collection + blur_scatt.blur_attribute = "TotalEnergyDeposit" + blur_scatt.blur_method = "InverseSquare" + blur_scatt.blur_reference_value = reference_energy + blur_scatt.blur_resolution = scatt_resolution + scatt_collection = blur_scatt.name + + blur_abs = sim.add_actor("DigitizerBlurringActor", "BlurAbs") + blur_abs.attached_to = absorber.name + blur_abs.input_digi_collection = abs_collection + blur_abs.blur_attribute = "TotalEnergyDeposit" + blur_abs.blur_method = "InverseSquare" + blur_abs.blur_reference_value = reference_energy + blur_abs.blur_resolution = abs_resolution + abs_collection = blur_abs.name + + # Time blurring + time_fwhm = 10 * ns + time_sigma = time_fwhm / 2.355 + time_scatt = sim.add_actor("DigitizerBlurringActor", "TimeBlurScatt") + time_scatt.attached_to = scatterer.name + time_scatt.input_digi_collection = scatt_collection + time_scatt.blur_attribute = "GlobalTime" + time_scatt.blur_method = "Gaussian" + time_scatt.blur_sigma = time_sigma + scatt_collection = time_scatt.name + + time_abs = sim.add_actor("DigitizerBlurringActor", "TimeBlurAbs") + time_abs.attached_to = absorber.name + time_abs.input_digi_collection = abs_collection + time_abs.blur_attribute = "GlobalTime" + time_abs.blur_method = "Gaussian" + time_abs.blur_sigma = time_sigma + abs_collection = time_abs.name + + # Energy windows (thresholds) + threshold_min_scatt = 70 * keV + threshold_min_abs = 70 * keV + threshold_max = 2.0 * MeV + thr_scatt = sim.add_actor("DigitizerEnergyWindowsActor", "ThrScatt") + thr_scatt.attached_to = scatterer.name + thr_scatt.input_digi_collection = scatt_collection + thr_scatt.channels = [ + {"name": thr_scatt.name, "min": threshold_min_scatt, "max": threshold_max} + ] + scatt_collection = thr_scatt.name + + thr_abs = sim.add_actor("DigitizerEnergyWindowsActor", "ThrAbs") + thr_abs.attached_to = absorber.name + thr_abs.input_digi_collection = abs_collection + thr_abs.channels = [ + {"name": thr_abs.name, "min": threshold_min_abs, "max": threshold_max} + ] + + # Saving root files + thr_scatt.output_filename = f"{thr_scatt.name}.root" + thr_abs.output_filename = f"{thr_abs.name}.root" + scatt_file = thr_scatt.get_output_path() + abs_file = thr_abs.get_output_path() + + return scatt_file, abs_file + + +def macaco1_merge_and_compute_coincidences( + scatt_root_filename, + abs_root_filename, + output_root_filename, + time_windows, + scatt_tree_name, + abs_tree_name, + merged_tree_name, +): + """ + Merge scatterer/absorber singles, run the Gate coincidence sorter, + then filter good coincidences. + """ + + ccmod_merge_several_singles_root_into_one( + scatt_root_filename, + abs_root_filename, + output_root_filename, + scatt_tree_name=scatt_tree_name, + abs_tree_name=abs_tree_name, + output_tree_name=merged_tree_name, + overwrite=True, + ) + + with uproot.open(output_root_filename) as f: + singles_tree = f[merged_tree_name] + coincidences = cc_coincidences_sorter(singles_tree, time_windows) + + if coincidences is None or coincidences.empty: + return coincidences + + coincidences = kill_multiple_coinc(coincidences) + coincidences = kill_same_volume_pairs(coincidences) + return coincidences + + +def load_exp_histograms( + path: Path, scatter_name, abs_name +) -> dict[str, tuple[np.ndarray, np.ndarray]]: + if not path.exists(): + raise FileNotFoundError(f"Experimental data file not found: {path}") + histograms: dict[str, tuple[np.ndarray, np.ndarray]] = {} + with uproot.open(path) as f: + for name in (scatter_name, abs_name): + obj = f[name] + counts, edges = obj.to_numpy() + histograms[name] = ( + np.asarray(counts, dtype=float), + np.asarray(edges, dtype=float), + ) + return histograms + + +def gaussian_fit_hist(centers: np.ndarray, counts: np.ndarray) -> tuple[float, float]: + mask = counts > 0 + if mask.sum() < 3: # need at least 3 points to fit a quadratic + return float("nan"), float("nan") + x = centers[mask].astype(float) # x values for the fit (bin centers) + y = np.log(counts[mask].astype(float)) # log of counts to linearize Gaussian + a, b, _ = np.polyfit(x, y, 2) # quadratic fit: y = a x^2 + b x + c + mu = -b / (2.0 * a) # mean of Gaussian from quadratic coefficients + sigma = float(np.sqrt(-1.0 / (2.0 * a))) # sigma from curvature + return float(mu), sigma # return fitted mean and sigma + + +def compare_peak_gaussian( + sim_energies: np.ndarray, + exp_counts: np.ndarray, + exp_edges: np.ndarray, + expected_energy: float, + label: str, + tol_frac: float = 0.20, + output_plot_path: Path = None, # New parameter +): + centers = 0.5 * (exp_edges[:-1] + exp_edges[1:]) + win_lo, win_hi = peak_window_from_exp( + centers, exp_counts, expected_energy=expected_energy + ) + exp_mask = (centers >= win_lo) & (centers <= win_hi) + exp_counts_win = exp_counts[exp_mask] + exp_centers_win = centers[exp_mask] + + if exp_counts_win.sum() <= 0.0: + print(f"✗ {label}: no experimental counts in window", flush=True) + return False + + sim_counts, _ = np.histogram(sim_energies, bins=exp_edges) + sim_counts_win = sim_counts[exp_mask] + + mu_exp, sigma_exp = gaussian_fit_hist(exp_centers_win, exp_counts_win) + mu_sim, sigma_sim = gaussian_fit_hist(exp_centers_win, sim_counts_win) + + # --- Plotting Logic --- + if output_plot_path: + plt.figure(figsize=(8, 6)) + + # Plot Data + plt.step( + exp_centers_win, + exp_counts_win, + where="mid", + label="Exp Data", + color="black", + alpha=0.6, + ) + plt.step( + exp_centers_win, + sim_counts_win, + where="mid", + label="Sim Data", + color="blue", + alpha=0.6, + ) + + # Plot Fits + x_fit = np.linspace(win_lo, win_hi, 100) + if not np.isnan(mu_exp): + y_exp = np.max(exp_counts_win) * np.exp( + -((x_fit - mu_exp) ** 2) / (2 * sigma_exp**2) + ) + plt.plot(x_fit, y_exp, "k--", label=f"Exp Fit ($\mu$={mu_exp:.1f})") + + if not np.isnan(mu_sim): + y_sim = np.max(sim_counts_win) * np.exp( + -((x_fit - mu_sim) ** 2) / (2 * sigma_sim**2) + ) + plt.plot(x_fit, y_sim, "b--", label=f"Sim Fit ($\mu$={mu_sim:.1f})") + + plt.title(f"Gaussian Fit Comparison: {label}") + plt.xlabel("Energy (keV)") + plt.ylabel("Counts") + plt.legend() + plt.grid(True, linestyle=":", alpha=0.7) + plt.savefig(output_plot_path) + plt.close() + print("Figure saved in: ", output_plot_path) + # ---------------------- + + if ( + np.isnan(mu_exp) + or np.isnan(sigma_exp) + or np.isnan(mu_sim) + or np.isnan(sigma_sim) + ): + print(f"✗ {label}: Gaussian fit failed", flush=True) + return False + + # (Keep the rest of your existing logic for validation...) + mean_pct = 100.0 * abs(mu_sim - mu_exp) / abs(mu_exp) + sigma_pct = 100.0 * abs(sigma_sim - sigma_exp) / abs(sigma_exp) + fwhm_pct = ( + 100.0 * abs((2.355 * sigma_sim) - (2.355 * sigma_exp)) / abs(2.355 * sigma_exp) + ) + tol_pct = tol_frac * 100.0 + + max_diff = max(mean_pct, sigma_pct, fwhm_pct) + if mean_pct < tol_pct and sigma_pct < tol_pct and fwhm_pct < tol_pct: + print(f"✓ {label} pass (max diff={max_diff:.1f}%)", flush=True) + return True + + print(f"✗ {label} fail (max diff={max_diff:.1f}%)", flush=True) + return False + + +# builds the energy window around the expected peak in the experimental histogram. +def peak_window_from_exp( + centers: np.ndarray, + counts: np.ndarray, + expected_energy: float, + rel_span: float = 0.1, +) -> tuple[float, float]: + if expected_energy < 700.0: + rel_span = 0.12 + else: + rel_span = 0.075 + low = expected_energy * (1.0 - rel_span) + high = expected_energy * (1.0 + rel_span) + mask = (centers > low) & (centers < high) + sel_counts = counts[mask] + if sel_counts.sum() <= 0.0: + raise ValueError(f"No experimental counts near {expected_energy} keV.") + return low, high diff --git a/opengate/contrib/compton_camera/macaco_materials.db b/opengate/contrib/compton_camera/macaco_materials.db new file mode 100644 index 000000000..9bd8fd933 --- /dev/null +++ b/opengate/contrib/compton_camera/macaco_materials.db @@ -0,0 +1,103 @@ +[Elements] +Hydrogen: S= H ; Z= 1. ; A= 1.01 g/mole +Carbon: S= C ; Z= 6. ; A= 12.01 g/mole +Nitrogen: S= N ; Z= 7. ; A= 14.01 g/mole +Oxygen: S= O ; Z= 8. ; A= 16.00 g/mole +Aluminium: S= Al ; Z= 13. ; A= 26.98 g/mole +Argon: S= Ar ; Z= 18. ; A= 39.95 g/mole +Titanium: S= Ti ; Z= 22. ; A= 47.867 g/mole +Manganese: S= Mn ; Z= 25. ; A= 54.938 g/mole +Iron: S= Fe ; Z= 26. ; A= 55.845 g/mole +Lead: S= Pb ; Z= 82. ; A= 207.20 g/mole +Silicon: S= Si ; Z= 14. ; A= 28.09 g/mole +Cadmium: S= Cd ; Z= 48. ; A= 112.41 g/mole +Tellurium: S= Te ; Z= 52. ; A= 127.6 g/mole +Zinc: S= Zn ; Z= 30. ; A= 65.39 g/mole +Tungsten: S= W ; Z= 74. ; A= 183.84 g/mole +Gallium: S= Ga ; Z= 31. ; A= 69.723 g/mole +Germanium: S= Ge ; Z= 32. ; A= 72.61 g/mole +Bromine: S= Br ; Z= 35. ; A= 79.90 g/mole +Yttrium: S= Y ; Z= 39. ; A= 88.91 g/mole +Cesium: S= Cs ; Z= 55. ; A= 132.905 g/mole +Lanthanum: S= La ; Z= 57. ; A= 138.91 g/mole +Cerium: S= Ce ; Z= 58. ; A= 140.116 g/mole +Gadolinium: S= Gd ; Z= 64. ; A= 157.25 g/mole +Lutetium: S= Lu ; Z= 71. ; A= 174.97 g/mole +Thallium: S= Tl ; Z= 81. ; A= 204.37 g/mole +Bismuth: S= Bi ; Z= 83. ; A= 208.98 g/mole +Uranium: S= U ; Z= 92. ; A= 238.03 g/mole + +[Materials] +Vacuum: d=0.000001 mg/cm3 ; n=1 + +el: name=Hydrogen ; n=1 + +Tungsten: d=19.3 g/cm3 ; n=1 ; state=solid + +el: name=auto ; n=1 + +Lead: d=11.4 g/cm3 ; n=1 ; state=solid + +el: name=auto ; n=1 + + +Aluminium: d=2.7 g/cm3 ; n=1 ; state=solid + +el: name=auto ; n=1 + +Plastic: d=1.18 g/cm3 ; n=3; state=solid + +el: name=Carbon ; n=5 + +el: name=Hydrogen ; n=8 + +el: name=Oxygen ; n=2 + +CZT: d=5.68 g/cm3 ; n=3; state=solid + +el: name=Cadmium ; n=9 + +el: name=Zinc ; n=1 + +el: name=Tellurium ; n=10 + + +LYSO-Ce-Hilger: d=7.10 g/cm3; n=5 ; state=Solid + +el: name=Lutetium ; f=0.713838658203075 + +el: name=Yttrium ; f=0.040302477781781 + +el: name=Silicon; f=0.063721807284236 + +el: name=Oxygen; f=0.181501252152072 + +el: name=Cerium; f=0.000635804578835201 + +Polyethylene: d=0.96 g/cm3 ; n=2 + +el: name=Hydrogen ; n=2 + +el: name=Carbon ; n=1 + +PVC: d=1.65 g/cm3 ; n=3 ; state=solid + +el: name=Hydrogen ; n=3 + +el: name=Carbon ; n=2 + +el: name=Chlorine ; n=1 + +SS304: d=7.92 g/cm3 ; n=4 ; state=solid + +el: name=Iron ; f=0.695 + +el: name=Chromium ; f=0.190 + +el: name=Nickel ; f=0.095 + +el: name=Manganese ; f=0.020 + +PTFE: d= 2.18 g/cm3 ; n=2 ; state=solid + +el: name=Carbon ; n=1 + +el: name=Fluorine ; n=2 + +Plexiglass: d=1.19 g/cm3; n=3; state=solid + +el: name=Hydrogen; f=0.080538 + +el: name=Carbon; f=0.599848 + +el: name=Oxygen; f=0.319614 + +LaBr3: d=5.29 g/cm3 ; n=2 + +el: name=Bromine; n=3 + +el: name=Lanthanum; n=1 + + +CeBr3: d=5.1 g/cm3 ; n=2 + +el: name=Bromine; n=3 + +el: name=Cerium; n=1 + + +LaBr3_VLC: d=5.06 g/cm3 ; n=2 + +el: name=Bromine; n=3 + +el: name=Lanthanum; n=1 + +PCB: d=1.20 g/cm3 ; n=3 + +el: name=Hydrogen; n=8 + +el: name=Carbon; n=5 + +el: name=Oxygen; n=2 diff --git a/opengate/tests/data b/opengate/tests/data index e0dc12194..93fd96a2f 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit e0dc12194098002f8ec8390a86961072cf63557d +Subproject commit 93fd96a2f7016b40c55f0f422781db0364f3b84d diff --git a/opengate/tests/src/chemistry/test101_chemical_counting_actor.py b/opengate/tests/src/chemistry/test101_chemical_counting_actor.py old mode 100644 new mode 100755 index 852b81cb4..75f328696 --- a/opengate/tests/src/chemistry/test101_chemical_counting_actor.py +++ b/opengate/tests/src/chemistry/test101_chemical_counting_actor.py @@ -283,6 +283,11 @@ def compare_results(result1, result2): if __name__ == "__main__": + + paths = utility.get_default_test_paths( + __file__, gate_folder="", output_folder="test101_chemical_counting_actor" + ) + sim_global, stats_global, chem_actor_global = create_simulation( use_actor_requested_dna_em=False ) @@ -290,11 +295,8 @@ def compare_results(result1, result2): results_global = chem_actor_global.results.get_data() print_results("Explicit region DNA EM results:", results_global) - plot_species_counts( - results_global, - "Explicit region DNA EM", - "test101_species_explicit_region_dna_em.png", - ) + fn = paths.output / "test101_species_explicit_region_dna_em.png" + plot_species_counts(results_global, "Explicit region DNA EM", fn) # sim_actor, stats_actor, chem_actor_actor = create_simulation( # use_actor_requested_dna_em=True diff --git a/opengate/tests/src/chemistry/test103_chemistry_world.py b/opengate/tests/src/chemistry/test103_chemistry_world.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/chemistry/test104_scavenger_minimal.py b/opengate/tests/src/chemistry/test104_scavenger_minimal.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/chemistry/test105_chemistry_confinement.py b/opengate/tests/src/chemistry/test105_chemistry_confinement.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/chemistry/test106_reaction_counter_timing.py b/opengate/tests/src/chemistry/test106_reaction_counter_timing.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/external/coresi/constants.yaml b/opengate/tests/src/external/coresi/constants.yaml new file mode 100644 index 000000000..ab7c57642 --- /dev/null +++ b/opengate/tests/src/external/coresi/constants.yaml @@ -0,0 +1,72 @@ +# SPDX-FileCopyrightText: 2024 Vincent Lequertier , Voichita Maxim +# +# SPDX-License-Identifier: MIT + +avogadro: 6.02214179e+23 +m_e: 510.99895000 # electron rest mass in keV +r_e: 2.8179403227e-13 # electron classical radius in cm + +# TODO: This depends on the material (Silicium currently). These values were +# computed by former students, Enrique or Yuemeng. Add new keys if we add other +# materials +doppler_broadening: + energy_140: + a1: 0.3773 + a2: 0.1443 + sigma_beta_1: 0.002090928993768552 # rad + sigma_beta_2: 0.01835770181680156 # rad + energy_245: + a1: 0.6291 + a2: 0.2335 + sigma_beta_1: 0.0014484966569412875 # rad + sigma_beta_2: 0.011749867543634142 # rad + energy_364: + a1: 0.8393 + a2: 0.2735 + sigma_beta_1: 0.0013006032426011388 # rad + sigma_beta_2: 0.009809405280466551 # rad + energy_511: + a1: 1.0000 + a2: 0.2980 + sigma_beta_1: 0.0012218929989592862 # rad + sigma_beta_2: 0.008513331113645892 # rad + energy_1275: + a1: 1.00000000 + a2: 0.21889611 + sigma_beta_1: 0.0010469077761151514 # rad + sigma_beta_2: 0.006182086979184491 # rad + energy_2170: + a1: 0.42537069 + a2: 0.33833700 + sigma_beta_1: 0.00040536484092147874 # rad + sigma_beta_2: 0.0035140981264683954 # rad + energy_4439: + a1: 0.60596747 + a2: 0.28269513 + sigma_beta_1: 7.811690920027629e-05 # rad + sigma_beta_2: 0.00273685099297266 # rad + energy_6130: + a1: 0.54893444 + a2: 0.23693385 + sigma_beta_1: 8.049907448898096e-05 # rad + sigma_beta_2: 0.002622310532968701 # rad + +materials: + Si: + density: 2.33 # g cm^-3 + moll_mass: 28.09 # g * cm^-3 + eff: 14 + # Comes from here https://www.physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Elem&Output2=Hand + NIST: [ [1.000E-03, 1.317E-02, 1.567E+03, 0.000E+00, 1.570E+03, 1.567E+03], [1.500E-03, 2.393E-02, 5.333E+02, 0.000E+00, 5.356E+02, 5.333E+02], [1.839E-03, 3.081E-02, 3.071E+02, 0.000E+00, 3.092E+02, 3.071E+02], [1.839E-03, 3.081E-02, 3.191E+03, 0.000E+00, 3.193E+03, 3.191E+03], [2.000E-03, 3.388E-02, 2.775E+03, 0.000E+00, 2.777E+03, 2.775E+03], [3.000E-03, 4.962E-02, 9.767E+02, 0.000E+00, 9.784E+02, 9.767E+02], [4.000E-03, 6.135E-02, 4.514E+02, 0.000E+00, 4.528E+02, 4.514E+02], [5.000E-03, 7.110E-02, 2.438E+02, 0.000E+00, 2.451E+02, 2.439E+02], [6.000E-03, 7.983E-02, 1.458E+02, 0.000E+00, 1.470E+02, 1.459E+02], [8.000E-03, 9.507E-02, 6.379E+01, 0.000E+00, 6.469E+01, 6.389E+01], [1.000E-02, 1.076E-01, 3.315E+01, 0.000E+00, 3.388E+01, 3.326E+01], [1.500E-02, 1.289E-01, 9.848E+00, 0.000E+00, 1.034E+01, 9.977E+00], [2.000E-02, 1.402E-01, 4.089E+00, 0.000E+00, 4.463E+00, 4.229E+00], [3.000E-02, 1.501E-01, 1.161E+00, 0.000E+00, 1.436E+00, 1.311E+00], [4.000E-02, 1.534E-01, 4.687E-01, 0.000E+00, 7.010E-01, 6.221E-01], [5.000E-02, 1.538E-01, 2.307E-01, 0.000E+00, 4.385E-01, 3.845E-01], [6.000E-02, 1.526E-01, 1.288E-01, 0.000E+00, 3.206E-01, 2.815E-01], [8.000E-02, 1.483E-01, 5.120E-02, 0.000E+00, 2.228E-01, 1.995E-01], [1.000E-01, 1.432E-01, 2.498E-02, 0.000E+00, 1.835E-01, 1.682E-01], [1.500E-01, 1.309E-01, 6.806E-03, 0.000E+00, 1.448E-01, 1.377E-01], [2.000E-01, 1.207E-01, 2.736E-03, 0.000E+00, 1.275E-01, 1.235E-01], [3.000E-01, 1.055E-01, 7.882E-04, 0.000E+00, 1.082E-01, 1.063E-01], [4.000E-01, 9.475E-02, 3.409E-04, 0.000E+00, 9.614E-02, 9.509E-02], [5.000E-01, 8.663E-02, 1.849E-04, 0.000E+00, 8.748E-02, 8.681E-02], [6.000E-01, 8.019E-02, 1.156E-04, 0.000E+00, 8.078E-02, 8.031E-02], [8.000E-01, 7.050E-02, 5.854E-05, 0.000E+00, 7.082E-02, 7.056E-02], [1.000E+00, 6.340E-02, 3.639E-05, 0.000E+00, 6.361E-02, 6.344E-02], [1.022E+00, 6.274E-02, 3.433E-05, 0.000E+00, 6.293E-02, 6.277E-02], [1.250E+00, 5.671E-02, 2.326E-05, 3.521E-05, 5.688E-02, 5.677E-02], [1.500E+00, 5.155E-02, 1.684E-05, 1.912E-04, 5.183E-02, 5.175E-02], [2.000E+00, 4.400E-02, 1.050E-05, 7.533E-04, 4.481E-02, 4.476E-02], [2.044E+00, 4.346E-02, 1.016E-05, 8.118E-04, 4.433E-02, 4.429E-02], [3.000E+00, 3.461E-02, 5.800E-06, 2.139E-03, 3.678E-02, 3.677E-02], [4.000E+00, 2.888E-02, 3.952E-06, 3.456E-03, 3.240E-02, 3.239E-02], [5.000E+00, 2.494E-02, 2.983E-06, 4.629E-03, 2.967E-02, 2.967E-02], [6.000E+00, 2.204E-02, 2.389E-06, 5.678E-03, 2.788E-02, 2.787E-02], [7.000E+00, 1.981E-02, 1.990E-06, 6.611E-03, 2.663E-02, 2.663E-02], [8.000E+00, 1.803E-02, 1.705E-06, 7.451E-03, 2.574E-02, 2.574E-02], [9.000E+00, 1.657E-02, 1.490E-06, 8.212E-03, 2.509E-02, 2.509E-02], [1.000E+01, 1.536E-02, 1.323E-06, 8.905E-03, 2.462E-02, 2.461E-02], [1.100E+01, 1.432E-02, 1.189E-06, 9.531E-03, 2.425E-02, 2.425E-02], [1.200E+01, 1.343E-02, 1.080E-06, 1.011E-02, 2.398E-02, 2.398E-02], [1.300E+01, 1.266E-02, 9.889E-07, 1.064E-02, 2.378E-02, 2.377E-02], [1.400E+01, 1.197E-02, 9.121E-07, 1.114E-02, 2.362E-02, 2.362E-02], [1.500E+01, 1.137E-02, 8.463E-07, 1.160E-02, 2.352E-02, 2.352E-02], [1.600E+01, 1.082E-02, 7.893E-07, 1.204E-02, 2.345E-02, 2.345E-02], [1.800E+01, 9.893E-03, 6.954E-07, 1.284E-02, 2.338E-02, 2.338E-02], [2.000E+01, 9.124E-03, 6.216E-07, 1.355E-02, 2.338E-02, 2.338E-02], [2.200E+01, 8.474E-03, 5.618E-07, 1.419E-02, 2.343E-02, 2.343E-02], [2.400E+01, 7.919E-03, 5.125E-07, 1.478E-02, 2.351E-02, 2.351E-02], [2.600E+01, 7.436E-03, 4.711E-07, 1.532E-02, 2.361E-02, 2.361E-02], [2.800E+01, 7.014E-03, 4.359E-07, 1.581E-02, 2.372E-02, 2.372E-02], [3.000E+01, 6.641E-03, 4.057E-07, 1.627E-02, 2.385E-02, 2.385E-02], [4.000E+01, 5.277E-03, 3.010E-07, 1.817E-02, 2.455E-02, 2.455E-02], [5.000E+01, 4.406E-03, 2.393E-07, 1.958E-02, 2.522E-02, 2.522E-02], [6.000E+01, 3.795E-03, 1.985E-07, 2.070E-02, 2.583E-02, 2.583E-02], [8.000E+01, 2.993E-03, 1.481E-07, 2.236E-02, 2.685E-02, 2.685E-02], [1.000E+02, 2.485E-03, 1.181E-07, 2.354E-02, 2.764E-02, 2.764E-02], [1.500E+02, 1.769E-03, 7.837E-08, 2.547E-02, 2.905E-02, 2.905E-02], [2.000E+02, 1.387E-03, 5.864E-08, 2.663E-02, 2.996E-02, 2.996E-02], [3.000E+02, 9.833E-04, 3.902E-08, 2.802E-02, 3.111E-02, 3.111E-02], [4.000E+02, 7.700E-04, 2.923E-08, 2.884E-02, 3.181E-02, 3.181E-02], [5.000E+02, 6.375E-04, 2.337E-08, 2.938E-02, 3.229E-02, 3.229E-02], [6.000E+02, 5.457E-04, 1.947E-08, 2.976E-02, 3.263E-02, 3.263E-02], [8.000E+02, 4.263E-04, 1.459E-08, 3.030E-02, 3.312E-02, 3.312E-02], [1.000E+03, 3.508E-04, 1.167E-08, 3.064E-02, 3.344E-02, 3.344E-02]] + BGO: + density: 7.13 # g cm^-3 + moll_mass: 1245.83 # g * cm^-3 + eff: 524 + # Comes from here https://www.physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Elem&Output2=Hand + NIST: [[1.000E-03, 4.741E-03, 4.680E+03, 0.000E+00, 4.690E+03, 4.680E+03], [1.103E-03, 5.573E-03, 3.813E+03, 0.000E+00, 3.823E+03, 3.813E+03], [1.217E-03, 6.540E-03, 3.109E+03, 0.000E+00, 3.118E+03, 3.109E+03], [1.232E-03, 6.674E-03, 3.640E+03, 0.000E+00, 3.650E+03, 3.640E+03], [1.248E-03, 6.811E-03, 3.623E+03, 0.000E+00, 3.632E+03, 3.623E+03], [1.328E-03, 7.526E-03, 3.483E+03, 0.000E+00, 3.493E+03, 3.483E+03], [1.414E-03, 8.301E-03, 3.099E+03, 0.000E+00, 3.108E+03, 3.099E+03], [1.500E-03, 9.082E-03, 2.843E+03, 0.000E+00, 2.852E+03, 2.843E+03], [2.000E-03, 1.370E-02, 1.477E+03, 0.000E+00, 1.485E+03, 1.477E+03], [2.580E-03, 1.901E-02, 8.092E+02, 0.000E+00, 8.174E+02, 8.092E+02], [2.633E-03, 1.949E-02, 1.496E+03, 0.000E+00, 1.504E+03, 1.496E+03], [2.688E-03, 1.998E-02, 1.506E+03, 0.000E+00, 1.515E+03, 1.506E+03], [3.000E-03, 2.271E-02, 1.571E+03, 0.000E+00, 1.579E+03, 1.571E+03], [3.177E-03, 2.420E-02, 1.355E+03, 0.000E+00, 1.363E+03, 1.355E+03], [3.427E-03, 2.629E-02, 1.275E+03, 0.000E+00, 1.283E+03, 1.275E+03], [3.696E-03, 2.842E-02, 1.057E+03, 0.000E+00, 1.064E+03, 1.057E+03], [3.845E-03, 2.953E-02, 1.012E+03, 0.000E+00, 1.019E+03, 1.012E+03], [3.999E-03, 3.075E-02, 9.201E+02, 0.000E+00, 9.271E+02, 9.201E+02], [4.000E-03, 3.076E-02, 9.558E+02, 0.000E+00, 9.628E+02, 9.558E+02], [5.000E-03, 3.781E-02, 5.529E+02, 0.000E+00, 5.592E+02, 5.529E+02], [6.000E-03, 4.402E-02, 3.509E+02, 0.000E+00, 3.565E+02, 3.509E+02], [8.000E-03, 5.448E-02, 1.689E+02, 0.000E+00, 1.734E+02, 1.689E+02], [1.000E-02, 6.292E-02, 9.493E+01, 0.000E+00, 9.872E+01, 9.499E+01], [1.110E-02, 6.687E-02, 7.222E+01, 0.000E+00, 7.566E+01, 7.228E+01], [1.221E-02, 7.042E-02, 7.982E+01, 0.000E+00, 8.296E+01, 7.989E+01], [1.342E-02, 7.394E-02, 6.245E+01, 0.000E+00, 6.530E+01, 6.252E+01], [1.500E-02, 7.805E-02, 9.156E+01, 0.000E+00, 9.409E+01, 9.163E+01], [1.571E-02, 7.972E-02, 8.089E+01, 0.000E+00, 8.330E+01, 8.097E+01], [1.605E-02, 8.047E-02, 1.017E+02, 0.000E+00, 1.040E+02, 1.018E+02], [1.639E-02, 8.122E-02, 9.662E+01, 0.000E+00, 9.892E+01, 9.670E+01], [2.000E-02, 8.804E-02, 6.576E+01, 0.000E+00, 6.758E+01, 6.584E+01], [3.000E-02, 1.004E-01, 2.251E+01, 0.000E+00, 2.363E+01, 2.261E+01], [4.000E-02, 1.070E-01, 1.037E+01, 0.000E+00, 1.116E+01, 1.048E+01], [5.000E-02, 1.103E-01, 5.646E+00, 0.000E+00, 6.238E+00, 5.756E+00], [6.000E-02, 1.117E-01, 3.422E+00, 0.000E+00, 3.894E+00, 3.534E+00], [8.000E-02, 1.118E-01, 1.546E+00, 0.000E+00, 1.884E+00, 1.658E+00], [9.053E-02, 1.111E-01, 1.098E+00, 0.000E+00, 1.393E+00, 1.209E+00], [9.053E-02, 1.111E-01, 4.804E+00, 0.000E+00, 5.099E+00, 4.915E+00], [1.000E-01, 1.102E-01, 3.705E+00, 0.000E+00, 3.971E+00, 3.815E+00], [1.500E-01, 1.039E-01, 1.281E+00, 0.000E+00, 1.462E+00, 1.385E+00], [2.000E-01, 9.742E-02, 5.983E-01, 0.000E+00, 7.415E-01, 6.957E-01], [3.000E-01, 8.659E-02, 2.073E-01, 0.000E+00, 3.157E-01, 2.939E-01], [4.000E-01, 7.842E-02, 1.004E-01, 0.000E+00, 1.915E-01, 1.788E-01], [5.000E-01, 7.206E-02, 5.852E-02, 0.000E+00, 1.389E-01, 1.306E-01], [6.000E-01, 6.693E-02, 3.833E-02, 0.000E+00, 1.111E-01, 1.053E-01], [8.000E-01, 5.906E-02, 2.037E-02, 0.000E+00, 8.280E-02, 7.943E-02], [1.000E+00, 5.322E-02, 1.284E-02, 0.000E+00, 6.824E-02, 6.606E-02], [1.022E+00, 5.267E-02, 1.229E-02, 0.000E+00, 6.704E-02, 6.496E-02], [1.250E+00, 4.767E-02, 8.293E-03, 2.784E-04, 5.765E-02, 5.624E-02], [1.500E+00, 4.337E-02, 5.907E-03, 1.339E-03, 5.160E-02, 5.062E-02], [2.000E+00, 3.706E-02, 3.572E-03, 4.095E-03, 4.528E-02, 4.472E-02], [2.044E+00, 3.660E-02, 3.443E-03, 4.340E-03, 4.491E-02, 4.438E-02], [3.000E+00, 2.918E-02, 1.868E-03, 9.118E-03, 4.042E-02, 4.017E-02], [4.000E+00, 2.435E-02, 1.223E-03, 1.321E-02, 3.897E-02, 3.883E-02], [5.000E+00, 2.103E-02, 8.960E-04, 1.667E-02, 3.877E-02, 3.868E-02], [6.000E+00, 1.859E-02, 7.020E-04, 1.966E-02, 3.914E-02, 3.908E-02], [7.000E+00, 1.671E-02, 5.748E-04, 2.231E-02, 3.981E-02, 3.977E-02], [8.000E+00, 1.521E-02, 4.855E-04, 2.468E-02, 4.062E-02, 4.059E-02], [9.000E+00, 1.398E-02, 4.196E-04, 2.684E-02, 4.153E-02, 4.150E-02], [1.000E+01, 1.295E-02, 3.692E-04, 2.882E-02, 4.246E-02, 4.244E-02], [1.100E+01, 1.208E-02, 3.291E-04, 3.066E-02, 4.342E-02, 4.340E-02], [1.200E+01, 1.133E-02, 2.969E-04, 3.235E-02, 4.436E-02, 4.435E-02], [1.300E+01, 1.068E-02, 2.703E-04, 3.395E-02, 4.531E-02, 4.530E-02], [1.400E+01, 1.010E-02, 2.480E-04, 3.544E-02, 4.623E-02, 4.622E-02], [1.500E+01, 9.590E-03, 2.291E-04, 3.682E-02, 4.710E-02, 4.709E-02], [1.600E+01, 9.133E-03, 2.128E-04, 3.813E-02, 4.796E-02, 4.796E-02]] + LaBr3: + density: 5.08 # g cm^-3 + moll_mass: 378.62 # g * cm^-3 + eff: 162 + # Comes from here https://www.physics.nist.gov/cgi-bin/Xcom/xcom2?Method=Elem&Output2=Hand + NIST: [[1.00E-003, 5.80E-003, 4.99E+003, 0.00E+000, 5.00E+003, 4.988E+03], [1.06E-003, 6.34E-003, 4.37E+003, 0.00E+000, 4.37E+003, 4.366E+03], [1.12E-003, 6.92E-003, 3.82E+003, 0.00E+000, 3.83E+003, 3.822E+03], [1.12E-003, 6.92E-003, 4.20E+003, 0.00E+000, 4.20E+003, 4.197E+03], [1.16E-003, 7.29E-003, 3.89E+003, 0.00E+000, 3.89E+003, 3.885E+03], [1.20E-003, 7.67E-003, 3.60E+003, 0.00E+000, 3.60E+003, 3.596E+03], [1.20E-003, 7.67E-003, 3.75E+003, 0.00E+000, 3.75E+003, 3.748E+03], [1.28E-003, 8.38E-003, 3.28E+003, 0.00E+000, 3.28E+003, 3.276E+03], [1.36E-003, 9.14E-003, 2.87E+003, 0.00E+000, 2.87E+003, 2.866E+03], [1.36E-003, 9.14E-003, 2.96E+003, 0.00E+000, 2.97E+003, 2.958E+03], [1.50E-003, 1.05E-002, 2.38E+003, 0.00E+000, 2.39E+003, 2.379E+03], [1.55E-003, 1.09E-002, 2.20E+003, 0.00E+000, 2.21E+003, 2.204E+03], [1.55E-003, 1.09E-002, 4.35E+003, 0.00E+000, 4.35E+003, 4.348E+03], [1.57E-003, 1.12E-002, 4.06E+003, 0.00E+000, 4.06E+003, 4.056E+03], [1.60E-003, 1.14E-002, 3.79E+003, 0.00E+000, 3.79E+003, 3.787E+03], [1.60E-003, 1.14E-002, 4.74E+003, 0.00E+000, 4.75E+003, 4.742E+03], [1.69E-003, 1.23E-002, 4.18E+003, 0.00E+000, 4.19E+003, 4.181E+03], [1.78E-003, 1.32E-002, 3.69E+003, 0.00E+000, 3.69E+003, 3.685E+03], [1.78E-003, 1.32E-002, 4.02E+003, 0.00E+000, 4.02E+003, 4.018E+03], [2.00E-003, 1.52E-002, 3.06E+003, 0.00E+000, 3.06E+003, 3.055E+03], [3.00E-003, 2.43E-002, 1.11E+003, 0.00E+000, 1.12E+003, 1.114E+03], [4.00E-003, 3.24E-002, 5.30E+002, 0.00E+000, 5.34E+002, 5.297E+02], [5.00E-003, 3.95E-002, 2.94E+002, 0.00E+000, 2.98E+002, 2.940E+02], [5.48E-003, 4.26E-002, 2.30E+002, 0.00E+000, 2.34E+002, 2.298E+02], [5.48E-003, 4.26E-002, 3.66E+002, 0.00E+000, 3.70E+002, 3.664E+02], [5.68E-003, 4.39E-002, 3.34E+002, 0.00E+000, 3.38E+002, 3.344E+02], [5.89E-003, 4.51E-002, 3.05E+002, 0.00E+000, 3.09E+002, 3.051E+02], [5.89E-003, 4.51E-002, 3.68E+002, 0.00E+000, 3.71E+002, 3.676E+02], [6.00E-003, 4.57E-002, 3.53E+002, 0.00E+000, 3.57E+002, 3.529E+02], [6.27E-003, 4.73E-002, 3.15E+002, 0.00E+000, 3.18E+002, 3.148E+02], [6.27E-003, 4.73E-002, 3.46E+002, 0.00E+000, 3.50E+002, 3.464E+02], [8.00E-003, 5.63E-002, 1.84E+002, 0.00E+000, 1.87E+002, 1.838E+02], [1.00E-002, 6.52E-002, 1.01E+002, 0.00E+000, 1.03E+002, 1.010E+02], [1.35E-002, 7.78E-002, 4.49E+001, 0.00E+000, 4.66E+001, 4.493E+01], [1.35E-002, 7.78E-002, 1.24E+002, 0.00E+000, 1.26E+002, 1.243E+02], [1.50E-002, 8.23E-002, 9.40E+001, 0.00E+000, 9.55E+001, 9.408E+01], [2.00E-002, 9.38E-002, 4.37E+001, 0.00E+000, 4.48E+001, 4.378E+01], [3.00E-002, 1.07E-001, 1.43E+001, 0.00E+000, 1.50E+001, 1.439E+01], [3.89E-002, 1.13E-001, 6.83E+000, 0.00E+000, 7.33E+000, 6.943E+00], [3.89E-002, 1.13E-001, 1.51E+001, 0.00E+000, 1.56E+001, 1.518E+01], [4.00E-002, 1.13E-001, 1.40E+001, 0.00E+000, 1.45E+001, 1.409E+01], [5.00E-002, 1.17E-001, 7.63E+000, 0.00E+000, 8.01E+000, 7.750E+00], [6.00E-002, 1.18E-001, 4.61E+000, 0.00E+000, 4.92E+000, 4.730E+00], [8.00E-002, 1.17E-001, 2.06E+000, 0.00E+000, 2.29E+000, 2.173E+00], [1.00E-001, 1.15E-001, 1.09E+000, 0.00E+000, 1.28E+000, 1.204E+00], [1.50E-001, 1.08E-001, 3.39E-001, 0.00E+000, 4.86E-001, 4.471E-01], [2.00E-001, 1.01E-001, 1.49E-001, 0.00E+000, 2.72E-001, 2.491E-01], [3.00E-001, 8.90E-002, 4.74E-002, 0.00E+000, 1.47E-001, 1.364E-01], [4.00E-001, 8.04E-002, 2.18E-002, 0.00E+000, 1.08E-001, 1.022E-01], [5.00E-001, 7.37E-002, 1.23E-002, 0.00E+000, 9.00E-002, 8.602E-02], [6.00E-001, 6.84E-002, 7.86E-003, 0.00E+000, 7.90E-002, 7.626E-02], [8.00E-001, 6.03E-002, 4.07E-003, 0.00E+000, 6.59E-002, 6.433E-02], [1.00E+000, 5.43E-002, 2.54E-003, 0.00E+000, 5.78E-002, 5.682E-02], [1.02E+000, 5.37E-002, 2.43E-003, 0.00E+000, 5.71E-002, 5.612E-02], [1.25E+000, 4.86E-002, 1.64E-003, 1.41E-004, 5.10E-002, 5.035E-02], [1.50E+000, 4.42E-002, 1.17E-003, 6.68E-004, 4.65E-002, 4.602E-02], [2.00E+000, 3.77E-002, 7.14E-004, 2.32E-003, 4.10E-002, 4.076E-02], [2.04E+000, 3.73E-002, 6.89E-004, 2.48E-003, 4.07E-002, 4.043E-02], [3.00E+000, 2.97E-002, 3.80E-004, 6.00E-003, 3.62E-002, 3.609E-02], [4.00E+000, 2.48E-002, 2.52E-004, 9.31E-003, 3.44E-002, 3.437E-02], [5.00E+000, 2.14E-002, 1.86E-004, 1.22E-002, 3.39E-002, 3.385E-02], [6.00E+000, 1.89E-002, 1.47E-004, 1.47E-002, 3.39E-002, 3.386E-02], [7.00E+000, 1.70E-002, 1.21E-004, 1.69E-002, 3.42E-002, 3.420E-02], [8.00E+000, 1.55E-002, 1.03E-004, 1.89E-002, 3.47E-002, 3.469E-02], [9.00E+000, 1.42E-002, 8.92E-005, 2.07E-002, 3.53E-002, 3.528E-02], [1.00E+001, 1.32E-002, 7.87E-005, 2.24E-002, 3.59E-002, 3.593E-02]] diff --git a/opengate/tests/src/external/coresi/test097a_ccmod_simulation.py b/opengate/tests/src/external/coresi/test097a_ccmod_simulation.py new file mode 100755 index 000000000..e968f8327 --- /dev/null +++ b/opengate/tests/src/external/coresi/test097a_ccmod_simulation.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.utility import g4_units +import opengate.tests.utility as utility +from scipy.spatial.transform import Rotation +import opengate.contrib.compton_camera.macaco as macaco + +if __name__ == "__main__": + + # get tests paths + paths = utility.get_default_test_paths( + __file__, gate_folder="", output_folder="test097_coresi_ccmod" + ) + output_folder = paths.output + data_folder = paths.data + + # units + m = g4_units.m + mm = g4_units.mm + cm = g4_units.cm + keV = g4_units.keV + Bq = g4_units.Bq + MBq = 1e6 * g4_units.Bq + sec = g4_units.s + ns = g4_units.ns + + # sim + sim = gate.Simulation() + sim.visu = False + sim.visu_type = "vrml" + sim.random_seed = "auto" + sim.output_dir = output_folder + sim.number_of_threads = 4 + sim.progress_bar = True + + # world + world = sim.world + world.size = [0.5 * m, 0.5 * m, 0.7 * m] + sim.world.material = "G4_AIR" + + # add one camera + name1 = "macaco1" + macaco1 = macaco.add_macaco1_camera(sim, name1) + camera1 = macaco1["camera"] + scatterer = macaco1["scatterer"] + absorber = macaco1["absorber"] + # Coresi requires camera in the -Z direction + camera1.rotation = Rotation.from_euler("y", 180, degrees=True).as_matrix() + camera1.translation = [0, 0, -10 * cm] + + # add the digitizer (output singles) + scatt_file, abs_file = macaco.add_macaco1_camera_digitizer(sim, scatterer, absorber) + print(f"Scatt file: {scatt_file}") + print(f"Abs file: {abs_file}") + + # stats + stats = sim.add_actor("SimulationStatisticsActor", "Stats") + stats.output_filename = "stats.txt" + + # physics + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.set_production_cut("world", "all", 1 * mm) + sim.physics_manager.set_production_cut(camera1.name, "all", 0.1 * mm) + # sim.physics_manager.set_production_cut(camera2.name, "all", 0.1 * mm) + + # source + source = sim.add_source("GenericSource", "src") + source.particle = "gamma" + source.energy.mono = 662 * keV + source.position.type = "sphere" + source.position.radius = 0.25 * mm + source.position.translation = [0, 0, 0] + source.direction.type = "iso" + source.activity = 0.847 * MBq / sim.number_of_threads + + # acquisition time + if sim.visu: + source.activity = 10 * Bq + sim.run_timing_intervals = [[0 * sec, 30 * sec]] + + # go + sim.run() + + # print stats + print(stats) diff --git a/opengate/tests/src/external/coresi/test097b_ccmod_coincidences.py b/opengate/tests/src/external/coresi/test097b_ccmod_coincidences.py new file mode 100755 index 000000000..0b85bf719 --- /dev/null +++ b/opengate/tests/src/external/coresi/test097b_ccmod_coincidences.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from opengate.utility import g4_units +import opengate.tests.utility as utility +from opengate.actors.coincidences import * +import opengate.contrib.compton_camera.macaco as macaco +import opengate.contrib.compton_camera.coresi_helpers as ch + +if __name__ == "__main__": + + # get tests paths + paths = utility.get_default_test_paths( + __file__, gate_folder="", output_folder="test097_coresi_ccmod" + ) + output_folder = paths.output + data_folder = paths.data + + # units + mm = g4_units.mm + cm = g4_units.cm + ns = g4_units.ns + + # compute the coincidences + scatt_file = output_folder / "ThrScatt.root" + abs_file = output_folder / "ThrAbs.root" + print(f"Computing coincidences from {scatt_file.name} and {abs_file.name}") + merge_file = output_folder / "singles.root" + coinc_file = output_folder / "coincidences.root" + coincidences = macaco.macaco1_merge_and_compute_coincidences( + scatt_file, + abs_file, + time_windows=12 * ns, + output_root_filename=merge_file, + scatt_tree_name="ThrScatt", + abs_tree_name="ThrAbs", + merged_tree_name="Singles", + ) + root_write_trees(coinc_file, ["Coincidences"], [coincidences]) + print(f"Merged singles file: {merge_file}") + print(f"Coincidences file: {coinc_file}") + print(f"Number of singles in coincidences {len(coincidences)}") + + # Consider the coincidences from the macaco simulation + # (no need to re-read the file, but this is to show an example) + coinc_file = output_folder / "coincidences.root" + if not coinc_file.exists(): + raise Exception(f"File {coinc_file} not found") + coincidences = uproot.open(coinc_file)["Coincidences"].arrays(library="pd") + print(f"Reading the coincidences file: {coinc_file.name}") + print(f"Number of singles in coincidences : {len(coincidences)}") + + # cones + print() + print("Computing cones from the coincidences ...") + data_cones = ccmod_make_cones(coincidences, energy_key_name="TotalEnergyDeposit") + cones_filename = output_folder / "cones.root" + root_write_trees(cones_filename, ["cones"], [data_cones]) + print(f"Cones file: {cones_filename.name}") + print(f"Number of cones (coincidences) : {len(data_cones)}") + + # CORESI stage1: convert the root file + print() + print(f"Converting root cones to coresi data file format ...") + data_filename = output_folder / "coincidences.dat" + ch.coresi_convert_root_data(cones_filename, "cones", data_filename) + print(f"Data file: {data_filename.name}") diff --git a/opengate/tests/src/external/coresi/test097c_coresi_reconstruction.py b/opengate/tests/src/external/coresi/test097c_coresi_reconstruction.py new file mode 100755 index 000000000..fee34cc6e --- /dev/null +++ b/opengate/tests/src/external/coresi/test097c_coresi_reconstruction.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.tests.utility as utility +from opengate.actors.coincidences import * +import opengate.contrib.compton_camera.coresi_helpers as ch + +if __name__ == "__main__": + + # get tests paths + paths = utility.get_default_test_paths( + __file__, gate_folder="", output_folder="test097_coresi_ccmod" + ) + output_folder = paths.output + data_folder = paths.data + + # Create default coresi config file + coresi_config = ch.coresi_new_config() + + # set the main coresi parameters + data_filename = output_folder / "coincidences.dat" + coresi_config["data_file"] = str(data_filename) + coresi_config["n_events"] = 1000 # max number of events to process + coresi_config["E0"] = [662] # in keV + coresi_config["energy_range"] = [100, 700] # in keV + + # set the camera parameters from macaco1 + coresi_config["cameras"]["n_cameras"] = 1 + + # scatterer layer (cm) volume name = "*_scatterer" (in holder1) + ca = coresi_config["cameras"]["common_attributes"] + ca["n_sca_layers"] = 1 + ca["sca_material"] = "LaBr3" + ca["sca_layer_0"]["center"] = [0, 0, -2.74 - 0.01] + ca["sca_layer_0"]["size"] = [2.72, 2.68, 0.50] + + # absorber layer (cm) volume name = "*_absorber" (in holder2) + ca["n_absorbers"] = 1 + ca["abs_material"] = "LaBr3" + ca["abs_layer_0"]["center"] = [0, 0, 2.51 - 0.01] + ca["abs_layer_0"]["size"] = [3.24, 3.60, 1] + + # (single) camera position + p = coresi_config["cameras"]["position_0"] + p["frame_origin"] = [0, 0, -10] # camera at -10 cm from the source + p["Ox"] = [1, 0, 0] # parallel to scatterer edge + p["Oy"] = [0, 1, 0] # parallel to scatterer edge + p["Oz"] = [0, 0, 1] # orthogonal to the camera, tw the source + + # reconstructed volume dimension + v = coresi_config["volume"] + v["volume_dimensions"] = [5, 5, 5] # total size in cm + v["n_voxels"] = [21, 21, 21] # nb of voxels in each dimension + v["volume_centre"] = [0, 0, 0] + + config_filename = output_folder / "coresi_config.yaml" + ch.coresi_write_config(coresi_config, config_filename) + print(f"Coresi config file: {config_filename}") + + # Mock the command line arguments + sys.argv = ["coresi", "-c", str(config_filename)] + + # prepare coresi run + try: + import coresi.main + import sys + except ModuleNotFoundError: + print( + "Coresi module not available. Please install coresi first (https://github.com/CoReSi-SPECT/coresi)" + ) + exit(1) + + # Run the reconstruction + print(f'Running coresi with command: "coresi -c {config_filename}"') + coresi.main.run() diff --git a/opengate/tests/src/geometry/test090_helpers.py b/opengate/tests/src/geometry/test090_helpers.py old mode 100755 new mode 100644 diff --git a/opengate/tests/src/geometry/test099_fields_analytical_B.py b/opengate/tests/src/geometry/test099_fields_analytical_B.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_analytical_E.py b/opengate/tests/src/geometry/test099_fields_analytical_E.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_api.py b/opengate/tests/src/geometry/test099_fields_api.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_custom_vs_native_B.py b/opengate/tests/src/geometry/test099_fields_custom_vs_native_B.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_custom_vs_native_E.py b/opengate/tests/src/geometry/test099_fields_custom_vs_native_E.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_mapped_multi_volume_refresh.py b/opengate/tests/src/geometry/test099_fields_mapped_multi_volume_refresh.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_mapped_repeated_placements.py b/opengate/tests/src/geometry/test099_fields_mapped_repeated_placements.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_mapped_rotated_volume.py b/opengate/tests/src/geometry/test099_fields_mapped_rotated_volume.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_multi_volume_refresh.py b/opengate/tests/src/geometry/test099_fields_multi_volume_refresh.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_repeated_placements.py b/opengate/tests/src/geometry/test099_fields_repeated_placements.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_rotated_volume.py b/opengate/tests/src/geometry/test099_fields_rotated_volume.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_serialization.py b/opengate/tests/src/geometry/test099_fields_serialization.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test099_fields_stepper.py b/opengate/tests/src/geometry/test099_fields_stepper.py old mode 100644 new mode 100755 index 95a57749d..5fbf6b8c4 --- a/opengate/tests/src/geometry/test099_fields_stepper.py +++ b/opengate/tests/src/geometry/test099_fields_stepper.py @@ -22,7 +22,7 @@ from opengate.geometry import fields from opengate.tests import utility -from test099_fields_helpers import ( +from opengate.tests.src.geometry.test099_fields_helpers import ( g4_m, g4_cm, g4_mm, diff --git a/opengate/tests/src/geometry/test099_fields_stepper_em.py b/opengate/tests/src/geometry/test099_fields_stepper_em.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/geometry/test107_macaco1_mt.py b/opengate/tests/src/geometry/test107_macaco1_mt.py new file mode 100755 index 000000000..5307fece9 --- /dev/null +++ b/opengate/tests/src/geometry/test107_macaco1_mt.py @@ -0,0 +1,796 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np +import pandas as pd +import uproot +import matplotlib +import matplotlib.pyplot as plt +from pathlib import Path +import time +import sys + +import opengate as gate +from opengate.tests import utility +from opengate.utility import g4_units +from opengate.contrib.compton_camera.macaco import ( + add_macaco1_camera, + add_macaco1_camera_digitizer, + macaco1_merge_and_compute_coincidences, + load_exp_histograms, + compare_peak_gaussian, +) +from opengate.actors.coincidences import cc_coincidences_sorter +from opengate.contrib.root_helpers import ( + root_tree_get_branch_data, + root_tree_get_branch_types, + root_write_tree, +) +from opengate.actors.coincidences import kill_multiple_coinc, kill_same_volume_pairs + +try: + from scipy.optimize import curve_fit +except Exception: + curve_fit = None + +# ----------------------------- +# Configuration +# ----------------------------- +ME_C2_KEV = 511.0 +ARM_UNITS = "rad" +ARM_RANGE = (-2.0, 2.0) +SOURCE_POS_MM = np.array([0.0, 0.0, 0.0], dtype=float) +ESUM_BIN_WIDTH_KEV = 5.0 +ESUM_MAX_KEV = 2000.0 +E_BIN_WIDTH_KEV = 5.0 +E_MAX_KEV = 1500.0 +ARM_REBIN_FACTOR = 2 +FWHM_PASS_THRESHOLD_PCT = 40.0 +COINC_TIME_WINDOW_NS = 50.0 + + +# ----------------------------- +# Helper functions (from working local test) +# ----------------------------- +def _open_root_tree(path: Path, *, retries: int = 5, delay_s: float = 0.5): + for attempt in range(1, retries + 1): + if not path.exists(): + time.sleep(delay_s) + continue + with uproot.open(str(path)) as f: + keys = list(f.keys()) + if keys: + tname = keys[0].split(";")[0] + return f[tname] + time.sleep(delay_s) + raise RuntimeError( + f"ROOT file has no trees after {retries} attempts: {path} " + f"(size={path.stat().st_size if path.exists() else 'missing'})" + ) + + +def load_root_as_dataframe(root_path: Path, *, out_dir: Path) -> pd.DataFrame: + root_path = root_path if root_path.exists() else out_dir / root_path.name + tree = _open_root_tree(root_path) + return tree.arrays(library="pd") + + +def load_digitizer_singles( + scatt_file: Path, abs_file: Path, *, out_dir: Path +) -> tuple[pd.DataFrame, pd.DataFrame]: + return ( + load_root_as_dataframe(scatt_file, out_dir=out_dir), + load_root_as_dataframe(abs_file, out_dir=out_dir), + ) + + +def prepare_layered_singles( + scatter_df: pd.DataFrame, absorber_df: pd.DataFrame +) -> pd.DataFrame: + def _tag_layer(df: pd.DataFrame, label: str) -> pd.DataFrame: + tagged = df.copy() + tagged["PreStepUniqueVolumeID"] = label + return tagged + + combined = pd.concat( + [_tag_layer(scatter_df, "scatterer"), _tag_layer(absorber_df, "absorber")], + ignore_index=True, + ) + return combined.sort_values(["GlobalTime", "EventID"]).reset_index(drop=True) + + +def build_singles_dataframe(layered: pd.DataFrame) -> pd.DataFrame: + required = [ + "EventID", + "GlobalTime", + "PreStepUniqueVolumeID", + "TotalEnergyDeposit", + "PostPosition_X", + "PostPosition_Y", + "PostPosition_Z", + ] + missing = [col for col in required if col not in layered.columns] + if missing: + raise ValueError(f"Missing required singles columns: {missing}") + return layered[required] + + +def write_singles_root(path: Path, df: pd.DataFrame) -> Path: + df_to_write = df.copy() + for col in df_to_write.columns: + if df_to_write[col].dtype == object: + df_to_write[col] = pd.Categorical(df_to_write[col]) + data = df_to_write.to_dict(orient="list") + data = root_tree_get_branch_data(data) + types = root_tree_get_branch_types(data) + with uproot.recreate(path) as f: + root_write_tree(f, "Singles", types, data) + return path + + +def extract_sorted_coincidences( + coincidences: pd.DataFrame, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + labels1 = coincidences["PreStepUniqueVolumeID1"].astype(str).str.lower() + is_scatt1 = labels1.str.contains("scatt") + scatt_pos = np.column_stack( + [ + np.where( + is_scatt1, + coincidences[f"PostPosition_{ax}1"], + coincidences[f"PostPosition_{ax}2"], + ) + for ax in "XYZ" + ] + ).astype(float) + abs_pos = np.column_stack( + [ + np.where( + is_scatt1, + coincidences[f"PostPosition_{ax}2"], + coincidences[f"PostPosition_{ax}1"], + ) + for ax in "XYZ" + ] + ).astype(float) + e_scatt = np.where( + is_scatt1, + coincidences["TotalEnergyDeposit1"], + coincidences["TotalEnergyDeposit2"], + ).astype(float) + e_abs = np.where( + is_scatt1, + coincidences["TotalEnergyDeposit2"], + coincidences["TotalEnergyDeposit1"], + ).astype(float) + return scatt_pos, abs_pos, e_scatt, e_abs, is_scatt1.to_numpy(dtype=bool) + + +def pairs_from_coincident_singles(df: pd.DataFrame) -> pd.DataFrame: + if "CoincID" not in df.columns: + raise ValueError("Missing CoincID in coincident singles.") + rows = [] + for _, group in df.groupby("CoincID"): + if len(group) < 2: + continue + s1, s2 = group.iloc[0], group.iloc[1] + rows.append( + { + "PreStepUniqueVolumeID1": s1["PreStepUniqueVolumeID"], + "PreStepUniqueVolumeID2": s2["PreStepUniqueVolumeID"], + **{f"PostPosition_{ax}1": s1[f"PostPosition_{ax}"] for ax in "XYZ"}, + **{f"PostPosition_{ax}2": s2[f"PostPosition_{ax}"] for ax in "XYZ"}, + "TotalEnergyDeposit1": s1["TotalEnergyDeposit"], + "TotalEnergyDeposit2": s2["TotalEnergyDeposit"], + } + ) + return pd.DataFrame(rows) + + +def compute_theta_c(e_scatt: np.ndarray, e0: np.ndarray) -> np.ndarray: + e_scattered = e0 - e_scatt + valid = (e0 > 0.0) & (e_scattered > 0.0) + cos_theta = np.full_like(e0, np.nan, dtype=float) + cos_theta[valid] = 1.0 - ( + ME_C2_KEV * e_scatt[valid] / (e_scattered[valid] * e0[valid]) + ) + in_range = (cos_theta >= -1.0) & (cos_theta <= 1.0) + cos_theta = np.where(in_range, cos_theta, np.nan) + return np.arccos(cos_theta) + + +def compute_theta_g( + scatter_pos: np.ndarray, absorber_pos: np.ndarray, source_pos: np.ndarray +) -> np.ndarray: + incoming = scatter_pos - source_pos + outgoing = absorber_pos - scatter_pos + incoming_norm = np.linalg.norm(incoming, axis=1) + outgoing_norm = np.linalg.norm(outgoing, axis=1) + valid = (incoming_norm > 0.0) & (outgoing_norm > 0.0) + cos_theta = np.full(incoming_norm.shape, np.nan, dtype=float) + dot = np.einsum("ij,ij->i", incoming, outgoing) + cos_theta[valid] = dot[valid] / (incoming_norm[valid] * outgoing_norm[valid]) + in_range = (cos_theta >= -1.0) & (cos_theta <= 1.0) + cos_theta = np.where(in_range, cos_theta, np.nan) + return np.arccos(cos_theta) + + +def compute_arm( + scatter_pos: np.ndarray, + absorber_pos: np.ndarray, + e_scatt: np.ndarray, + e_abs: np.ndarray, + source_pos: np.ndarray, +) -> np.ndarray: + e0 = np.where(e_scatt + e_abs < 600.0, 511.0, 1275.0) + theta_c = compute_theta_c(e_scatt, e0) + theta_g = compute_theta_g(scatter_pos, absorber_pos, source_pos) + valid = ~np.isnan(theta_c) & ~np.isnan(theta_g) + return (theta_g - theta_c)[valid] + + +def normalize_hist_counts(counts: np.ndarray, edges: np.ndarray) -> np.ndarray: + total = float(np.sum(counts)) + if total <= 0: + return np.zeros_like(counts, dtype=float) + bin_width = float(np.diff(edges)[0]) + return counts.astype(float) / (total * bin_width) + + +def rebin_hist( + counts: np.ndarray, edges: np.ndarray, factor: int +) -> tuple[np.ndarray, np.ndarray]: + if factor <= 1: + return counts, edges + n = (len(counts) // factor) * factor + if n == 0: + return counts, edges + return counts[:n].reshape(-1, factor).sum(axis=1), edges[: n + 1 : factor] + + +def lorentzian( + x: np.ndarray, amp: float, x0: float, gamma: float, offset: float +) -> np.ndarray: + half_gamma = 0.5 * gamma + return amp * half_gamma**2 / ((x - x0) ** 2 + half_gamma**2) + offset + + +def fit_lorentzian( + x: np.ndarray, y: np.ndarray +) -> tuple[np.ndarray, np.ndarray] | tuple[None, None]: + if curve_fit is None or np.count_nonzero(y) < 2: + return None, None + amp0 = float(np.max(y)) + x0 = float(x[np.argmax(y)]) + above_half = x[y >= amp0 / 2.0] + gamma0 = ( + float(above_half[-1] - above_half[0]) + if len(above_half) >= 2 + else float((x[-1] - x[0]) / 10.0) + ) + gamma0 = max(gamma0, float(np.diff(x).mean())) + offset0 = float(np.min(y)) + p0 = [amp0, x0, gamma0, offset0] + bin_w = float(np.diff(x).mean()) if len(x) > 1 else 1e-3 + amp_max = float(np.max(y)) + try: + popt, pcov = curve_fit( + lorentzian, + x, + y, + p0=p0, + bounds=( + [0, x.min(), bin_w, -np.inf], + [amp_max, x.max(), np.ptp(x), np.inf], + ), + maxfev=10000, + ) + except RuntimeError: + return None, None + return popt, pcov + + +def load_experimental_arm_histograms( + path: Path, +) -> dict[str, tuple[np.ndarray, np.ndarray]]: + path = Path(path) + if not path.exists(): + raise FileNotFoundError(f"Experimental ROOT file not found: {path}") + mapping = { + "all": "AngularUncertaintyReal", + "511": "AngularUncertaintyReal511", + "1275": "AngularUncertaintyReal1275", + } + histograms: dict[str, tuple[np.ndarray, np.ndarray]] = {} + with uproot.open(path) as f: + for key, name in mapping.items(): + if name not in f: + continue + counts, edges = f[name].to_numpy() + centers = 0.5 * (edges[:-1] + edges[1:]) + histograms[key] = ( + np.asarray(centers, dtype=float), + np.asarray(counts, dtype=float), + ) + if not histograms: + raise ValueError(f"No AngularUncertainty* histograms found in {path}") + return histograms + + +def plot_arm_overlay( + label: str, + exp_centers: np.ndarray, + exp_counts: np.ndarray, + sim_arm: np.ndarray, + out_dir: Path, + *, + sim_bin_factor: int = 1, +) -> tuple | None: + if exp_centers.size < 2: + print(f"Skipping {label} ARM (empty exp histogram).") + return None + bin_width = float(np.median(np.diff(exp_centers))) + edges = np.concatenate( + [exp_centers - bin_width / 2.0, [exp_centers[-1] + bin_width / 2.0]] + ) + exp_counts_rb, edges_rb = rebin_hist(exp_counts, edges, ARM_REBIN_FACTOR) + exp_hist = normalize_hist_counts(exp_counts_rb, edges_rb) + sim_edges = edges_rb[::sim_bin_factor] + if sim_edges[-1] != edges_rb[-1]: + sim_edges = np.append(sim_edges, edges_rb[-1]) + sim_counts_rb, _ = np.histogram(sim_arm, bins=sim_edges) + sim_hist = normalize_hist_counts(sim_counts_rb, sim_edges) + exp_centers_rb = 0.5 * (edges_rb[:-1] + edges_rb[1:]) + sim_centers = 0.5 * (sim_edges[:-1] + sim_edges[1:]) + + fit_mask_exp = np.abs(exp_centers_rb) <= 0.2 + fit_mask_sim = np.abs(sim_centers) <= 0.2 + exp_fit, _ = fit_lorentzian(exp_centers_rb[fit_mask_exp], exp_hist[fit_mask_exp]) + sim_fit, _ = fit_lorentzian(sim_centers[fit_mask_sim], sim_hist[fit_mask_sim]) + + exp_mu = float(exp_fit[1]) if exp_fit is not None else None + exp_fwhm = float(abs(exp_fit[2])) if exp_fit is not None else None + sim_mu = float(sim_fit[1]) if sim_fit is not None else None + sim_fwhm = float(abs(sim_fit[2])) if sim_fit is not None else None + + exp_label = ( + f"Experimental (μ={exp_mu:.3f}, FWHM={exp_fwhm:.3f} {ARM_UNITS})" + if exp_fit is not None + else "Experimental" + ) + sim_label = ( + f"Simulated (μ={sim_mu:.3f}, FWHM={sim_fwhm:.3f} {ARM_UNITS})" + if sim_fit is not None + else "Simulated" + ) + + matplotlib.use("Agg") + fig, ax = plt.subplots(figsize=(7.0, 5.5)) + ax.plot(exp_centers_rb, exp_hist, color="tab:red", linewidth=1.2, label=exp_label) + ax.fill_between( + exp_centers_rb, exp_hist, 0.0, color="tab:red", alpha=0.25, linewidth=0 + ) + ax.plot( + sim_centers, + sim_hist, + color="black", + linestyle="--", + linewidth=1.1, + label=sim_label, + ) + ax.fill_between(sim_centers, sim_hist, 0.0, color="0.6", alpha=0.2, linewidth=0) + if exp_fit is not None: + xfit = np.linspace(exp_centers_rb.min(), exp_centers_rb.max(), 600) + ax.plot( + xfit, lorentzian(xfit, *exp_fit), color="tab:red", linewidth=1.4, alpha=0.9 + ) + if sim_fit is not None: + xfit = np.linspace(sim_centers.min(), sim_centers.max(), 600) + ax.plot( + xfit, + lorentzian(xfit, *sim_fit), + color="black", + linewidth=1.2, + alpha=0.9, + linestyle=":", + ) + ax.set_xlabel(r"$\theta_G - \theta_C$ ({})".format(ARM_UNITS)) + ax.set_ylabel("Normalised counts") + ax.set_title(f"ARM overlay ({label})") + ax.set_xlim(ARM_RANGE) + ax.set_ylim(bottom=0) + ax.grid(True, linestyle="--", alpha=0.2) + ax.legend(frameon=False) + fig.tight_layout() + out_path = out_dir / f"arm_overlay_{label.replace(' ', '_')}.png" + fig.savefig(out_path, dpi=200) + plt.close(fig) + print(f"Saved ARM overlay: {out_path}") + + return exp_centers_rb, exp_hist, sim_mu, sim_fwhm, exp_fwhm + + +# ----------------------------- +# Main +# ----------------------------- +def main(): + paths = utility.get_default_test_paths( + __file__, gate_folder="", output_folder="test099_macaco1" + ) + output_folder = paths.output + output_ref = paths.output_ref + + exp_singles_file = output_ref / "singles_experimental.root" + print(exp_singles_file) + + # ====================================================== + # 1) Create simulation + # ====================================================== + sim = gate.Simulation() + sim.visu = False + sim.number_of_threads = 1 + sim.check_volumes_overlap = False + sim.progress_bar = True + sim.output_dir = output_folder + sim.random_seed = "auto" + + m = g4_units.m + mm = g4_units.mm + keV = g4_units.keV + Bq = g4_units.Bq + sec = g4_units.s + ns = g4_units.ns + + # ====================================================== + # 2) Geometry + # ====================================================== + sim.world.size = [1 * m, 1 * m, 1 * m] + cam = add_macaco1_camera(sim) + scatterer = cam["scatterer"] + absorber = cam["absorber"] + camera = cam["camera"] + camera.translation = [0, 0, 83 * mm] + + # ====================================================== + # 3) Source + # ====================================================== + src_holder = sim.add_volume("Sphere", "test_source_holder") + src_holder.mother = sim.world + src_holder.material = "Plastic" + src_holder.rmax = 0.25 * mm + + src = sim.add_source("GenericSource", "Na22_decay") + src.particle = "ion 11 22" + src.attached_to = src_holder + src.activity = 847e3 * Bq + src.position.type = "point" + src.direction.type = "iso" + src.user_particle_life_time = 0 + + sim.physics_manager.enable_decay = True + + # ====================================================== + # 4) Timing + # ====================================================== + if sim.visu: + sim.run_timing_intervals = [[0, 0.000003 * sec]] + else: + sim.run_timing_intervals = [[0, 7 * sec]] + + # ====================================================== + # 5) Digitizer + # ====================================================== + scatt_file, abs_file = add_macaco1_camera_digitizer(sim, scatterer, absorber) + print(f"Scatt file: {scatt_file}") + print(f"Abs file: {abs_file}") + + # ====================================================== + # 6) Run simulation + # ====================================================== + sim.run() + + # ====================================================== + # 7) VALIDATION : SINGLES + # ====================================================== + with uproot.open(scatt_file) as f: + scatt = f[f.keys()[0]] + E_scatt = np.asarray(scatt["TotalEnergyDeposit"].array()) / keV + + with uproot.open(abs_file) as f: + absr = f[f.keys()[0]] + E_abs = np.asarray(absr["TotalEnergyDeposit"].array()) / keV + + scatt_name = "E_ly1" + abs_name = "E_ly2" + exp_hists = load_exp_histograms(exp_singles_file, scatt_name, abs_name) + exp_abs_counts, exp_abs_edges = exp_hists[abs_name] + exp_scatt_counts, exp_scatt_edges = exp_hists[scatt_name] + + is_ok = True + print() + b = compare_peak_gaussian( + E_abs, + exp_abs_counts, + exp_abs_edges, + expected_energy=1274.5, + label="Absorber 1274.5 keV", + output_plot_path=output_folder / "test099_abs_peak_1274.5keV.png", + ) + is_ok = is_ok and b + print() + b = compare_peak_gaussian( + E_scatt, + exp_scatt_counts, + exp_scatt_edges, + expected_energy=1274.5, + label="Scatterer 1274.5 keV", + output_plot_path=output_folder / "test099_scatt_peak_1274.5keV.png", + ) + is_ok = is_ok and b + print() + compare_peak_gaussian( + E_abs, + exp_abs_counts, + exp_abs_edges, + expected_energy=511.0, + label="Absorber 511 keV", + output_plot_path=output_folder / "test099_abs_peak_511keV.png", + ) + is_ok = is_ok and b + print() + compare_peak_gaussian( + E_scatt, + exp_scatt_counts, + exp_scatt_edges, + expected_energy=511.0, + label="Scatterer 511 keV", + output_plot_path=output_folder / "test099_scatt_peak_511keV.png", + ) + is_ok = is_ok and b + print("✓ MACACO1 singles energy test completed") + + print(f"FIXME the absolute counts is not correct ? normalization issue ? ") + + +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +from opengate.contrib.compton_camera.macaco import * + + +def main(): + # get tests paths + paths = utility.get_default_test_paths( + __file__, gate_folder="", output_folder="test107_macaco1" + ) + output_folder = paths.output + output_ref = paths.output_ref + + exp_singles_file = output_ref / "singles_experimental.root" + print(exp_singles_file) + + # ====================================================== + # 1) Create simulation + # ====================================================== + sim = gate.Simulation() + sim.visu = False + sim.number_of_threads = 4 + sim.check_volumes_overlap = False + sim.progress_bar = True + sim.output_dir = output_folder + sim.random_seed = 123456789 + + m = g4_units.m + mm = g4_units.mm + keV = g4_units.keV + Bq = g4_units.Bq + sec = g4_units.s + ns = g4_units.ns + + # ====================================================== + # 2) Geometry + # ====================================================== + sim.world.size = [1 * m, 1 * m, 1 * m] + cam = add_macaco1_camera(sim) + scatterer = cam["scatterer"] + absorber = cam["absorber"] + camera = cam["camera"] + # 53 mm from the first plane + camera.translation = [0, 0, 83 * mm] + + # ====================================================== + # 3) Source + # ====================================================== + src_holder = sim.add_volume("Sphere", "test_source_holder") + src_holder.mother = sim.world + src_holder.material = "Plastic" + src_holder.rmax = 0.25 * mm + + src = sim.add_source("GenericSource", "test_gammas") + src.particle = "gamma" + src.attached_to = src_holder + src.activity = 847e3 * Bq + src.position.type = "sphere" + src.position.radius = 0.25 * mm + src.direction.type = "iso" + + # Na-22-like spectrum + src.energy.type = "spectrum_discrete" + src.energy.spectrum_energies = [1274.5 * keV, 511 * keV] + src.energy.spectrum_weights = [0.9994, 1.807] + + # ====================================================== + # 4) Timing + # ====================================================== + if sim.visu: + sim.run_timing_intervals = [[0, 0.000003 * sec]] + else: + sim.run_timing_intervals = [[0, 5 * sec]] + + # ====================================================== + # 5) Digitizer + # ====================================================== + scatt_file, abs_file = add_macaco1_camera_digitizer(sim, scatterer, absorber) + print(f"Scatt file: {scatt_file}") + print(f"Abs file: {abs_file}") + + # ====================================================== + # 6) Run simulation + # ====================================================== + sim.run() + + # ====================================================== + # 7) VALIDATION : SINGLES + # ====================================================== + # Load data + with uproot.open(scatt_file) as f: + scatt = f[f.keys()[0]] + E_scatt = np.asarray(scatt["TotalEnergyDeposit"].array()) / keV + + with uproot.open(abs_file) as f: + absr = f[f.keys()[0]] + E_abs = np.asarray(absr["TotalEnergyDeposit"].array()) / keV + + # Energy-only tests (Gaussian peaks) + scatt_name = "E_ly1" # Scatterer histogram name + abs_name = "E_ly2" # Absorber histogram name + exp_hists = load_exp_histograms(exp_singles_file, scatt_name, abs_name) + exp_abs_counts, exp_abs_edges = exp_hists[abs_name] + exp_scatt_counts, exp_scatt_edges = exp_hists[scatt_name] + + is_ok = True + print() + b = compare_peak_gaussian( + E_abs, + exp_abs_counts, + exp_abs_edges, + expected_energy=1274.5, + label="Absorber 1274.5 keV", + output_plot_path=output_folder / "test099_abs_peak_1274.5keV.png", + ) + is_ok = is_ok and b + print() + b = compare_peak_gaussian( + E_scatt, + exp_scatt_counts, + exp_scatt_edges, + expected_energy=1274.5, + label="Scatterer 1274.5 keV", + output_plot_path=output_folder / "test099_scatt_peak_1274.5keV.png", + ) + is_ok = is_ok and b + print() + compare_peak_gaussian( + E_abs, + exp_abs_counts, + exp_abs_edges, + expected_energy=511.0, + label="Absorber 511 keV", + output_plot_path=output_folder / "test099_abs_peak_511keV.png", + ) + is_ok = is_ok and b + print() + compare_peak_gaussian( + E_scatt, + exp_scatt_counts, + exp_scatt_edges, + expected_energy=511.0, + label="Scatterer 511 keV", + output_plot_path=output_folder / "test099_scatt_peak_511keV.png", + ) + is_ok = is_ok and b + print("✓ MACACO1 singles energy test completed") + + # ====================================================== + # 8) VALIDATION : COINCIDENCES + # ====================================================== + coinc_file = output_folder / "coincidences.root" + coincidences = macaco1_merge_and_compute_coincidences( + scatt_file, + abs_file, + time_windows=12 * ns, + output_root_filename=coinc_file, + scatt_tree_name="ThrScatt", + abs_tree_name="ThrAbs", + merged_tree_name="Singles", + ) + + # ====================================================== + # 8b) VALIDATION : COINCIDENCES ARM + # ====================================================== + + exp_coinc_root = output_ref / "new_coinc_exp_data.root" + + if not exp_coinc_root.exists(): + print( + f"Experimental coincidence ROOT not found: {exp_coinc_root}, skipping ARM validation." + ) + else: + scatt_df, abs_df = load_digitizer_singles( + scatt_file, abs_file, out_dir=output_folder + ) + layered = prepare_layered_singles(scatt_df, abs_df) + singles_path = write_singles_root( + output_folder / "singles_sim.root", build_singles_dataframe(layered) + ) + + time_window = COINC_TIME_WINDOW_NS * ns + with uproot.open(singles_path) as f: + coinc_singles = cc_coincidences_sorter(f["Singles"], time_window) + + if coinc_singles is None or coinc_singles.empty: + raise AssertionError("No coincidences after sorter") + + coinc_singles = kill_multiple_coinc(coinc_singles, group_col="CoincID") + coinc_singles = kill_same_volume_pairs( + coinc_singles, group_col="CoincID", volume_col="PreStepUniqueVolumeID" + ) + coinc_pairs = pairs_from_coincident_singles(coinc_singles) + + scatt_pos, abs_pos, e_scatt, e_abs, _ = extract_sorted_coincidences(coinc_pairs) + e_scatt /= keV + e_abs /= keV + mask_1275 = (e_scatt + e_abs) >= 600.0 + + exp_arm_hists = load_experimental_arm_histograms(exp_coinc_root) + sim_arm = compute_arm( + scatt_pos[mask_1275], + abs_pos[mask_1275], + e_scatt[mask_1275], + e_abs[mask_1275], + SOURCE_POS_MM, + ) + centers, counts = exp_arm_hists["1275"] + + fit_stats = plot_arm_overlay( + "1275", centers, counts, sim_arm, output_folder, sim_bin_factor=2 + ) + + if fit_stats is None: + print("✗ ARM plot failed") + is_ok = False + else: + _, _, sim_mu, sim_fwhm, exp_fwhm = fit_stats + if sim_fwhm is not None and exp_fwhm is not None: + diff_pct = abs(sim_fwhm - exp_fwhm) / exp_fwhm * 100.0 + print( + f"1275 FWHM: sim={sim_fwhm:.3f} {ARM_UNITS}, " + f"exp={exp_fwhm:.3f} {ARM_UNITS}, diff={diff_pct:.1f}% " + f"(threshold={FWHM_PASS_THRESHOLD_PCT:.1f}%)" + ) + b = diff_pct < FWHM_PASS_THRESHOLD_PCT + print( + f"{'✓' if b else '✗'} 1275 FWHM {'within' if b else 'exceeds'} tolerance" + ) + is_ok = is_ok and b + else: + print("ARM FWHM fit unavailable — FAIL") + is_ok = False + + utility.test_ok(is_ok) + + +if __name__ == "__main__": + main() diff --git a/opengate/tests/src/source/test035b_dose_rate_vrt.py b/opengate/tests/src/source/test035b_dose_rate_vrt.py old mode 100644 new mode 100755