Open justingm opened 7 months ago
Hi,
I've also tried your script to add polyethylene and polystyrene but I get errors when I run the compute_coefficients.py file.
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
Dear,
Sorry for the late reply.
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.
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
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!!
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.
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.
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:
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
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!!