Update DataExporters.py to deal with an issue when all particles have the same kinetic energy.

[MarcBHahn]

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.

[bwheelz36]

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)

[MarcBHahn]

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. But in my understanding all phase spaces where the 'Ek [Mev]' is exactly the same for all particles should reproduce it, I believe.

Thanks for looking into it.

Bests Marc

[MarcBHahn]

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. Maybe that's somehow related to the python version and a change in behavior of min()/max() for a list of strictly equal int's or even floats similar within a certain 'delta_epsilon'.

[bwheelz36]

Hey @MarcBHahn Yes, I imagine it might definitely be python version specific. Let me know if you can reproduce - I'm happy to accept the fix but i'd like to be able to reproduce first to make sure we understand the underlying cause.

[brendan-whelan-seetreat]

@MarcBHahn will close this for now, please reopen if you are able to reproduce the bug