-
Notifications
You must be signed in to change notification settings - Fork 8
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Update DataExporters.py to deal with an issue when all particles have the same kinetic energy. #168
base: main
Are you sure you want to change the base?
Conversation
Fix the problem that a min() and max during header creation on 'Ek [MeV]' creates an empty list when all particles have the exact same kinetic energy by setting as default value that energy which euqals the average in this case.
Use mean instead of average for Pandas Series
Hey @MarcBHahn, thanks for reporting this. Can you provide an example that reproduces the bug? I've tried and I can't. my attempt is below: from pathlib import Path
from ParticlePhaseSpace import DataLoaders, DataExporters
from ParticlePhaseSpace import PhaseSpace
from ParticlePhaseSpace import ParticlePhaseSpaceUnits
from get_all_files import get_all_files
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np
# generate some data:
Nparticles = 10000
x = np.random.randn(Nparticles) # normal distributed random data
y = np.random.randn(Nparticles) # normal distributed random data
z = np.ones(Nparticles) * 100 # let's say z = 100 mm for arguments sake
px = np.zeros(x.shape) # zero divergence
py = np.zeros(x.shape)
pz = np.ones(Nparticles) * 11 # mono energetic
weight = np.ones(Nparticles) # all particles weighted equally
demo_data = pd.DataFrame(
{'x [mm]': x,
'y [mm]': y,
'z [mm]': z,
'px [MeV/c]': px,
'py [MeV/c]': py,
'pz [MeV/c]': pz,
'particle type [pdg_code]': (np.ones(x.shape, dtype=int) * 11),
'weight': np.ones(x.shape),
'particle id': np.arange(len(x)),
'time [ps]': np.zeros(x.shape)})
data = DataLoaders.Load_PandasData(demo_data)
PS = PhaseSpace(data)
export_path = Path('~/Documents').expanduser()
DataExporters.Topas_Exporter(PS, output_location=export_path,
output_name='test.phsp')
# now test reading it back in:
data_loc = export_path / 'test.phsp'
data = DataLoaders.Load_TopasData(data_loc)
PS2 = PhaseSpace(data)
# now you can do whatever tests you want, PS and PS2 seem to be the same to me (within rounding errors) |
It happened to me for a purely (synthetic) monoenergetic beam when all electrons hat exactly the same kinetic energy. So for all electrons the value electron_PS.ps_data['Ek [MeV]'] was exactly the same. In this case the min() and max() function return an empty list. Leading to a subsequent error. This behaviour is fixed by the PR mentioned above. I will try to provide a minimal working example in the next days, as soon as I find the time. Thanks for looking into it. Bests Marc |
It seems I can't reproduce it on that device here, which has a more recent (maybe slightly different) python environment, compared to the device i was getting the error on - which i currently don't have access to. I will try to reproduce it with the old setup. |
Hey @MarcBHahn |
@MarcBHahn will close this for now, please reopen if you are able to reproduce the bug |
The problem that a min() and max during header creation on 'Ek [MeV]' creates an empty list when all particles have the exact same kinetic energy and aborts.
This is fixed by setting a default value for the output of the min and max function.
This default values is set to the average Ekin which is for an monoenergetic beam equal to max and min.