Skip to content

Commit

Permalink
Fixed sequence_to_stills script to account for multi-pass data
Browse files Browse the repository at this point in the history
  • Loading branch information
@PrinceWalnut
PrinceWalnut committed Nov 13, 2023
1 parent ee6d3d8 commit 844bbc5
Showing 1 changed file with 70 additions and 78 deletions.
148 changes: 70 additions & 78 deletions src/laue_dials/command_line/sequence_to_stills.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,100 +60,92 @@ def sequence_to_stills(experiments, reflections, params):

new_experiments = []
new_reflections = []

first_loop = True
crystals = []
imagesets = []
id_counter = 0

print("Building experiment lists.")
for expt_id, experiment in enumerate(experiments):
# Generate crystals and imagesets for each scan point on first loop
if first_loop:
for i_scan_point in range(*experiment.scan.get_array_range()):
# Get the goniometer setting matrix
goniometer_setting_matrix = matrix.sqr(
experiment.goniometer.get_setting_rotation()
for i_scan_point in range(*experiment.scan.get_array_range()):
# Get the goniometer setting matrix
goniometer_setting_matrix = matrix.sqr(
experiment.goniometer.get_setting_rotation()
)
goniometer_axis = matrix.col(experiment.goniometer.get_rotation_axis())
step = experiment.scan.get_oscillation()[1]

# The A matrix is the goniometer setting matrix for this scan point
# times the scan varying A matrix at this scan point. Note, the
# goniometer setting matrix for scan point zero will be the identity
# matrix and represents the beginning of the oscillation.
# For stills, the A matrix needs to be positioned in the midpoint of an
# oscillation step. Hence, here the goniometer setting matrixis rotated
# by a further half oscillation step.
A = (
goniometer_axis.axis_and_angle_as_r3_rotation_matrix(
angle=experiment.scan.get_angle_from_array_index(i_scan_point)
+ (step / 2),
deg=True,
)
goniometer_axis = matrix.col(experiment.goniometer.get_rotation_axis())
step = experiment.scan.get_oscillation()[1]

# The A matrix is the goniometer setting matrix for this scan point
# times the scan varying A matrix at this scan point. Note, the
# goniometer setting matrix for scan point zero will be the identity
# matrix and represents the beginning of the oscillation.
# For stills, the A matrix needs to be positioned in the midpoint of an
# oscillation step. Hence, here the goniometer setting matrixis rotated
# by a further half oscillation step.
A = (
goniometer_axis.axis_and_angle_as_r3_rotation_matrix(
angle=experiment.scan.get_angle_from_array_index(i_scan_point)
+ (step / 2),
deg=True,
)
* goniometer_setting_matrix
* matrix.sqr(experiment.crystal.get_A())
* goniometer_setting_matrix
* matrix.sqr(experiment.crystal.get_A())
)
crystal = MosaicCrystalSauter2014(experiment.crystal)
crystal.set_A(A)

# Copy in mosaic parameters if available
if params.output.domain_size_ang is None and hasattr(
experiment.crystal, "get_domain_size_ang"
):
crystal.set_domain_size_ang(
experiment.crystal.get_domain_size_ang()
)
crystal = MosaicCrystalSauter2014(experiment.crystal)
crystal.set_A(A)

# Copy in mosaic parameters if available
if params.output.domain_size_ang is None and hasattr(
experiment.crystal, "get_domain_size_ang"
):
crystal.set_domain_size_ang(
experiment.crystal.get_domain_size_ang()
)
elif params.output.domain_size_ang is not None:
crystal.set_domain_size_ang(params.output.domain_size_ang)

if params.output.half_mosaicity_deg is None and hasattr(
experiment.crystal, "get_half_mosaicity_deg"
):
crystal.set_half_mosaicity_deg(
experiment.crystal.get_half_mosaicity_deg()
)
elif params.output.half_mosaicity_deg is not None:
crystal.set_half_mosaicity_deg(params.output.half_mosaicity_deg)

# Add to list of crystals
crystals.append(crystal)

# Add to list of imagesets
imagesets.append(
experiment.imageset.as_imageset()[i_scan_point : i_scan_point + 1]
elif params.output.domain_size_ang is not None:
crystal.set_domain_size_ang(params.output.domain_size_ang)

if params.output.half_mosaicity_deg is None and hasattr(
experiment.crystal, "get_half_mosaicity_deg"
):
crystal.set_half_mosaicity_deg(
experiment.crystal.get_half_mosaicity_deg()
)
elif params.output.half_mosaicity_deg is not None:
crystal.set_half_mosaicity_deg(params.output.half_mosaicity_deg)

# Don't regenerate this list
first_loop = False

# Split experiment by scan points
for i_scan_point in range(*experiment.scan.get_array_range()):
# Append experiment
new_experiment = Experiment(
identifier=str(i_scan_point),
identifier=str(id_counter),
detector=experiment.detector,
beam=experiment.beam,
crystal=crystals[i_scan_point],
imageset=imagesets[i_scan_point],
crystal=crystal,
imageset=experiment.imageset.as_imageset()[i_scan_point : i_scan_point + 1],
)
new_experiments.append(new_experiment)

# Increment ID for experiment
id_counter = id_counter + 1

# ----------------EXPERIMENTS CREATED---------------------------------
print("Building reflection table.")
for i_scan_point in trange(len(crystals)):
# Get subset of reflections on this image
_, _, _, _, z1, z2 = reflections["bbox"].parts()
subrefls = reflections.select((i_scan_point >= z1) & (i_scan_point < z2))
new_refls = subrefls.copy()
centroids = np.asarray(subrefls["xyzobs.px.value"] - [0.0, 0.0, 0.5])
centroids[:, 2] = i_scan_point
new_refls["xyzobs.px.value"] = flex.vec3_double(centroids)
new_refls["imageset_id"] = flex.int(len(new_refls), i_scan_point)
x, y, _ = subrefls["xyzobs.mm.value"].parts()
new_refls["xyzobs.mm.value"] = flex.vec3_double(
x, y, flex.double(len(new_refls), 0)
)
new_refls["id"] = flex.int(len(new_refls), i_scan_point)
new_reflections.append(new_refls)
id_counter = 0
for expt_id, experiment in enumerate(experiments):
# Get subset of reflections for this pass
pass_refls = reflections.select(reflections['imageset_id'] == expt_id)
for i_scan_point in range(*experiment.scan.get_array_range()):
# Get subset of reflections on this image
_, _, _, _, z1, z2 = pass_refls["bbox"].parts()
subrefls = pass_refls.select((i_scan_point >= z1) & (i_scan_point < z2))
new_refls = subrefls.copy()
centroids = np.asarray(subrefls["xyzobs.px.value"] - [0.0, 0.0, 0.5])
centroids[:, 2] = i_scan_point
new_refls["xyzobs.px.value"] = flex.vec3_double(centroids)
new_refls["imageset_id"] = flex.int(len(new_refls), i_scan_point)
x, y, _ = subrefls["xyzobs.mm.value"].parts()
new_refls["xyzobs.mm.value"] = flex.vec3_double(
x, y, flex.double(len(new_refls), 0)
)
new_refls["id"] = flex.int(len(new_refls), id_counter)
id_counter = id_counter + 1
new_reflections.append(new_refls)
return (new_experiments, new_reflections)


Expand Down

0 comments on commit 844bbc5

Please sign in to comment.