bwheelz36 / ParticlePhaseSpace

Import, analysis, and export of particle phase space data
https://bwheelz36.github.io/ParticlePhaseSpace/
GNU General Public License v3.0
11 stars 8 forks source link

Problem with conversion of particle name to PDG codes in TOPAS phase space export #169

Closed MarcBHahn closed 3 months ago

MarcBHahn commented 3 months ago

Hi, I have to report a small bug in the Topas phase space export code:

When I read a phase space based on the example provided in the documentation and export it as a csv with the build in function everything works fine: data = NewDataLoader(fileName, particle_type='electrons') exportName = fileName.strip(".dat") DataExporters.CSV_Exporter(PS, output_location=Path('.'), output_name=exportName)

When I replace the CSV export with the Topas Phase space expert DataExporters.Topas_Exporter(PS, output_location=Path('.'), output_name=exportName)

I get the following error message:

File /usr/lib/python3/dist-packages/numpy/lib/npyio.py:1610 in savetxt raise TypeError("Mismatch between array dtype ('%s') and " TypeError: Mismatch between array dtype ('object') and format specifier ('%11.5f %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f %2d %2d %2d')

This seems to be related to the fact that the line in the topas phase export code (in the Topas_Exporter class of DataExporters.py)

self._PS.ps_data['particle type [pdg_code]'].to_numpy(), third_direction_flag, first_particle_flag] Data = np.transpose(Data) FormatSpec = ['%11.5f', '%11.5f', '%11.5f', '%11.5f', '%11.5f', '%11.5f', '%11.5f', '%2d', '%2d', '%2d'] np.savetxt(WritefilePath, Data, fmt=FormatSpec, delimiter=' ')

returns a string ("electrons"), instead of the PDG code (11).

Because, when i exchange the related part "'%2d'" the FormatSpect with '%2s',

FormatSpec = ['%11.5f', '%11.5f', '%11.5f', '%11.5f', '%11.5f', '%11.5f', '%11.5f', '%2d', '%2d', '%2d']

the error message dissapears and the phase space is written with "electrons" instead of the PDG code of 11. But Topas expects PDG codes instead of a descriptive string.

So I wonder how can I get in this line

self._PS.ps_data['particle type [pdg_code]'].to_numpy()

the PDG code as integer instead of the string?

Thanks a lot.

Bests Marc

bwheelz36 commented 3 months ago

thanks for using the software and reporting the bug, I will look into it :-)

bwheelz36 commented 3 months ago

OK @MarcBHahn - I think you found two things here:

  1. There is a small error in the example of creating a new data loader - particles should be read in using pdg codes, not names. This doesn't effect the core code, it's the example which has the error in it.
  2. This should have been caught by an automatic check carried out at the data read in stage, but that check was failing to detect the issue.

I've fixed both things, and released a new version of the code (version 0.6.0). if you look at the docs again, you'll find the following line is changed in the example you were looking at:

old:

self.data[self._columns['particle type']]= particle_cfg.particle_properties[self._particle_type]['name']   

new:

self.data[self._columns['particle type']]= particle_cfg.particle_properties[self._particle_type]['pdg_code'] 

thanks again for reporting! if this fixes your problem can you please close the issue?

MarcBHahn commented 3 months ago

Hi @bwheelz36 ,

thanks for so quickly taking care of that.

I just tried it out, now properly using the pdg as int with data = NewDataLoader(data_loc, particle_type=11) which leads me to the following error message:

InputBeamMarc_for_simulations/./ParticlePhaseSpace/ParticlePhaseSpace/utilities.py:20 in _check_particle_types raise AttributeError(f'unknown particle type {particle_cfg.particle_properties[code]["name"]}') AttributeError: unknown particle type electrons

which is related to that code block

`

for code in pdg_codes:

    try:
        assert isinstance(code, (int, np.integer))
        test = particle_cfg.particle_properties[code]
    except:
        raise AttributeError(f'unknown particle type {particle_cfg.particle_properties[code]["name"]}')`

where assert isinstance fails and the test variable is never used.

Thanks for looking into it again. Bests Marc

bwheelz36 commented 3 months ago

sorry for the confusion. You should still instantiate with 'particle_type="electrons"'.

If you follow the (updated) example it shows how this works: https://bwheelz36.github.io/ParticlePhaseSpace/new_data_loader.html

Of course, this is just an example of a DataLoader - you could write it so it accepts pdg code instead.

MarcBHahn commented 3 months ago

Sorry for the delay in responding. I will try to test it as soon as possible - hopefully during the next week. Bests Marc

MarcBHahn commented 3 months ago

I can confirm now that the combination of self.data[self._columns['particle type']]= particle_cfg.particle_properties[self._particle_type]['pdg_code'] and 'particle_type="electrons" works as expected. Thanks a lot!