extract energy spread

[justingm]

Hi,

I've just started using georges and I am still a bit confused on how things work.
Is it possible to extract the energy spread of the beam after a degrader? If so, can you guide me on how?
Thanks!!

[justingm]

Hi,

I've also tried your script to add polyethylene and polystyrene but I get errors when I run the compute_coefficients.py file.

[justingm]

I am also wondering about the losses. I tried using 2.5 cm of water as degrader for 60 MeV protons and I lose 45% of the beam. I would expect that all protons would go through

[rtesse]

Dear,

Sorry for the late reply.

• Yes, it is possible to extract the energy spread. Using the BeamObserver, you have access to the momentum deviation and then, you can use the following relation to have the energy spread. Here is an example on how to get the momentum deviation.

import georges
from georges import ureg as _ureg
from georges.manzoni import Input, observers
from georges.manzoni.beam import MadXBeam

c1 = georges.Element.Degrader(NAME='DEG',
MATERIAL=georges.fermi.materials.Water,
L = 2.5*_ureg.cm,
WITH_LOSSES=True
)

sequence = georges.PlacementSequence(name="Sequence")
sequence.place(c1, at_entry=0 * _ureg.m);

kin = georges.Kinematics(60 * _ureg.MeV, particle=georges.particles.Proton, kinetic=True)
sequence.metadata.kinematics = kin

beam = MadXBeam(
kinematics=kin,
distribution=georges.Distribution.from_5d_multigaussian_distribution(
n=10000).distribution.values,
)

mi = Input.from_sequence(sequence=sequence)
mi.adjust_energy(input_energy=kin.ekin)
mi.freeze()
beam_observer_beam = mi.track(beam=beam, observers=observers.BeamObserver(with_input_beams=True))
dpp = beam_observer_beam.to_df().loc['DEG']['BEAM_OUT'][:,4]

• For the second point, do you have an issue with Georges? Could you post the error / the code?

• The losses in Georges are computed in such a way the Gaussian queue is removed. It means that after the degrader, the beam is Gaussian in each coordinate (x, x’, y, and y’). You can find the description of the algorithm here:

https://www.sciencedirect.com/science/article/pii/S0168583X22002336?via%3Dihub

For example, I computed the horizontal beam distribution with BDSIM and compared it with Georges. As you can see, the beam is cut in Georges and this explains the difference in the losses.

I hope it answers your question. If not, do not hesitate to comment this issue.

Regards.

[justingm]

Hi,

Thanks for the response.
I did figure out my first question eventually but I was using a different material (polyethylene) which was not on the database so I'm not getting any change in the energy spread. This led to my second question because I tried to follow your BDSIM simulation (bdsim_input.gmad) and ran compute_coefficients.py but I get the following error

File ".\compute_coefficients.py", line 105, in
data_transmission = results[["Epost", "Transmission_cut"]]
File "C:\Users\P290425\AppData\Roaming\Python\Python38\site-packages\pandas\core\frame.py", line 3810, in getitem
indexer = self.columns._get_indexer_strict(key, "columns")[1]
File "C:\Users\P290425\AppData\Roaming\Python\Python38\site-packages\pandas\core\indexes\base.py", line 6111, in _get_indexer_strict
self._raise_if_missing(keyarr, indexer, axis_name)
File "C:\Users\P290425\AppData\Roaming\Python\Python38\site-packages\pandas\core\indexes\base.py", line 6171, in _raise_if_missing
raise KeyError(f"None of [{key}] are in the [{axis_name}]")
KeyError: "None of [Index(['Epost', 'Transmission_cut'], dtype='object')] are in the [columns]"

I don't really understand what went wrong. Most likely it's from my side. I can also send you the bdsim output if that would help.

As for the losses, I get it now why the transmission was so low.
Thanks!!

[rtesse]

Hello,

Could you send me your BDSIM input file ? You can also send me the output file, I will try to take a look next week.

Regards.

[justingm]

Hi,

Here's the input filebdsim_for_georges.zip. The output is a bit too big to upload here. I can send it through email if that's more convenient.

[rtesse]

Hello,

I looked into your input. It has no issue. The problem is coming from the Python script.
Here is a corrected version: 





import os
import sys

import georges_core.kinematics as gkin
import numpy as np
import pandas as pd
from pybdsim.DataUproot import BDSIMOutput
from compute_quantiles import compute_quantile_2dim, compute_quantile_4dim
from lmfit import Model, Parameters

from georges import ureg as _ureg


def cut_data(filepath: str = None, epos: float = 100, matname: str = None, nprimary: int = 10) -> pd.DataFrame:
    """

    Args:
        filepath (str): Path to the file
        epos (float): Theoritical energy at the exit of the degrader
        matname (str): Name of the material
        nprimary (int): Number of particles launched in BDSIM.

    Returns:

    """

    element = "DEG"
    coordinates = ["x", "y", "xp", "yp", "energy", "p", "partID", "parentID", "S"]

    data = BDSIMOutput(filepath).event.samplers[element].df
    data = data[coordinates].copy(deep=True)

    data.reset_index(drop=True, inplace=True)

    data.query("parentID == 0 and partID == 2212", inplace=True)  # only get the primaries
    data["energy"] = gkin.etot_to_ekin(data["energy"].values * _ureg.GeV).m_as("MeV")
    data["momentum"] = data["p"] * 1000  # In Mev/c
    data.loc[:, "xp"] *= 1000
    data.loc[:, "yp"] *= 1000
    data.loc[:, "x"] *= 1000
    data.loc[:, "y"] *= 1000

    # Cut the data
    if epos < 220:
        quantile, data_cutted = compute_quantile_4dim(data)

    else:
        quantile, data_cutted = compute_quantile_2dim(data)

    # Properties of the cutted data
    momentum_cut = data_cutted["momentum"].mean()
    deviation_cut = (data_cutted["momentum"] - momentum_cut) / momentum_cut
    dpp_mean_cut = np.mean(deviation_cut)
    dpp_rms_cut = np.std(deviation_cut)
    transmission_cut = len(data_cutted) / nprimary

    # Return the results
    return pd.DataFrame(
        data={
            "MatName": [matname],
            "Epost": [epos],
            "DPPmean_cut": [dpp_mean_cut],
            "DPPrms_cut": [dpp_rms_cut],
            "Transmission_cut": [transmission_cut],
        },
    )


def fit_values(x, **params):
    val = 0.0
    parnames = sorted(params.keys())
    for i, pname in enumerate(parnames):
        val += params[pname] * x**i
    return val


if __name__ == "__main__":
    path = sys.argv[1]
    nparticles = int(sys.argv[2])

    results = []
    for file in os.listdir(path):
        if "root" in file:
            print(f"Process {file}")
            print(file.split("_"))
            results.append(
                    cut_data(
                        filepath=os.path.join(path, file),
                        epos=float(file.split("_")[-2].replace("E", "").replace("root", "")),
                        matname=file.split("_")[2],
                        nprimary=nparticles,
            ))
    results = pd.concat(results, axis=0)

    if len(results) < 2:
        print(results)
        sys.exit()

    params_transmission = Parameters()
    params_transmission.add("C0", value=1)
    params_transmission.add("C1", value=1)
    params_transmission.add("C2", value=1)
    params_transmission.add("C3", value=0)

    fit_data = Model(fit_values)
    data_transmission = results[["Epost", "Transmission_cut"]]
    fit_transmission = fit_data.fit(
        data_transmission["Transmission_cut"],
        params_transmission,
        x=data_transmission["Epost"],
    )
    coeff_transmission = fit_transmission.best_values

    params_dpp = Parameters()
    params_dpp.add("C0", value=1)
    params_dpp.add("C1", value=1)
    params_dpp.add("C2", value=1)
    params_dpp.add("C3", value=0)
    params_dpp.add("C4", value=0)
    params_dpp.add("C5", value=0)
    data_dpp = results[["Epost", "DPPrms_cut"]]
    fit_dpp = fit_data.fit(data_dpp["DPPrms_cut"], params_dpp, x=data_dpp["Epost"])
    coeff_dpp = fit_dpp.best_values

    print(
        f"""
    materiaName,
    {coeff_transmission['C0']},{coeff_transmission['C1']},{coeff_transmission['C2']},{coeff_transmission['C3']}
    {coeff_dpp['C0']},{coeff_dpp['C1']},{coeff_dpp['C2']},{coeff_dpp['C3']},{coeff_dpp['C4']},{coeff_dpp['C5']}
    """,
    )

To test, I put all the results in a folder res with the following structure:

• res
  • G4Al
    • output_G4_Al_E60_bic.root

Feel free to modify the scripts for your convenance.

If you have more than one result file, the script will make an interpolation that you have to add to the file georges/fermi/bdsim/data.csv

Regards