Skip to content

Commit

Permalink
Kassiopeia: KassiopeiaReader example for updated user guide
Browse files Browse the repository at this point in the history
  • Loading branch information
zykure committed Apr 16, 2021
1 parent 5642006 commit f012ccb
Show file tree
Hide file tree
Showing 2 changed files with 35 additions and 38 deletions.
69 changes: 33 additions & 36 deletions Kassiopeia/Python/Examples/QuadrupoleTrapAnalysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,61 +13,58 @@
# create reader instance
reader = KassiopeiaReader.Iterator(sys.argv[1])

# retrieve list of valid step ranges
step_presence = reader.getTree('component_step_cell_PRESENCE')
step_valid = zip(*[step_presence['INDEX'], step_presence['LENGTH']])

# retrieve list of step ranges per track
track_step_index = list(zip(*[reader.getTracks('FIRST_STEP_INDEX'), reader.getTracks('LAST_STEP_INDEX')]))
#print("#tracks:", len(track_step_index))

# load step-data tree
reader.loadTree('component_step_cell')
#print("#tracks:", len(reader.getTracks('TRACK_INDEX')))
#print('#steps: ', len(reader))
#print('fields: ', dir(reader))

# select a single output field
reader.select('orbital_magnetic_moment')

# retrieve list of valid step ranges
step_presence = reader.getTree('component_step_cell_PRESENCE')
step_valid = zip(*[step_presence['INDEX'], step_presence['LENGTH']])

# create step iterator
step = iter(reader)
#print('fields: ', dir(step))

# iterate over track indices
for track_index in reader.getTracks('TRACK_INDEX'):
for first_step_index, last_step_index in track_step_index:

# retrieve index of first/last step for current track
first_step_index = reader.getTracks('FIRST_STEP_INDEX')[track_index]
last_step_index = reader.getTracks('LAST_STEP_INDEX')[track_index]
max_moment = -np.inf
min_moment = np.inf

# adjust first/last step range to include only valid entries
for first_valid,valid_length in step_valid:
last_valid = first_valid + valid_length - 1
if first_valid >= first_step_index:
first_step_index = first_valid
if last_valid < last_step_index:
last_step_index = last_valid
break
# iterate over steps in given range
for step in iter(reader):

#print("track #{} valid interval: {} .. {}".format(track_index, first_step_index, last_step_index))
step_index = reader.iev - 1

# advance iterator to first step
while reader.iev < first_step_index:
item = next(step)
# advance iterator to first step in this track
if step_index < first_step_index:
continue

max_moment = -float('Inf')
min_moment = float('Inf')
# only process steps in valid range within this track
for first_valid,valid_length in step_valid:

# iterate over steps in given range
while reader.iev <= last_step_index:
last_valid = first_valid + valid_length - 1
if step_index >= first_valid and step_index <= last_valid:

moment = float(step.orbital_magnetic_moment)
if moment > max_moment:
max_moment = moment
if moment < min_moment:
min_moment = moment

# advance step iterator
item = next(step)
if first_valid > first_step_index:
break

moment = float(item.orbital_magnetic_moment)
if moment > max_moment:
max_moment = moment
if moment < min_moment:
min_moment = moment
# stop iteration at the end of this track
if step_index >= last_step_index:
break

# calculate result for current track
deviation = 2.0 * (max_moment - min_moment) / (max_moment + min_moment)
print("extrema for track <{:g}>".format(deviation))

break
4 changes: 2 additions & 2 deletions Kassiopeia/XML/Examples/VTK/QuadrupoleTrapSimulation.xml
Original file line number Diff line number Diff line change
Expand Up @@ -367,9 +367,9 @@
<command parent="write_root" field="add_track_output" child="component_track_world"/>
<command parent="write_root" field="add_step_output" child="component_step_world"/>
<command parent="write_vtk" field="set_track_point" child="component_track_position"/>
<command parent="write_vtk" field="set_track_data" child="component_track_length"/>
<command parent="write_vtk" field="set_track_data" child="component_track_world"/>
<command parent="write_vtk" field="set_step_point" child="component_step_position"/>
<command parent="write_vtk" field="set_step_data" child="component_step_length"/>
<command parent="write_vtk" field="set_step_data" child="component_step_world"/>
<geo_surface name="surface_start" surfaces="world/quadrupole_trap/start">
<command parent="root_terminator" field="add_terminator" child="term_start"/>
</geo_surface>
Expand Down

0 comments on commit f012ccb

Please sign in to comment.