aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2024-02-27 11:35:46 +0100
committerThomas White <taw@physics.org>2024-02-27 11:35:46 +0100
commit0ebaf0eb0465be23d8fee4088271cb92154a606d (patch)
tree74a085215b78ad1b2e7df9b1e504510a795fcca0
parent92e4b1fcfedb9315bdd002fb5258ce67860519ae (diff)
Julia: Move merging utils to separate module
-rw-r--r--julia/CrystFEL/src/CrystFEL.jl4
-rw-r--r--julia/CrystFEL/src/mergeutils.jl57
-rw-r--r--julia/process_hkl.jl52
3 files changed, 61 insertions, 52 deletions
diff --git a/julia/CrystFEL/src/CrystFEL.jl b/julia/CrystFEL/src/CrystFEL.jl
index 5ecfab91..897d2614 100644
--- a/julia/CrystFEL/src/CrystFEL.jl
+++ b/julia/CrystFEL/src/CrystFEL.jl
@@ -87,4 +87,8 @@ include("peaksearch.jl")
using .PeakSearch
export zaefpeaks, peakfinder8, peakfinder9
+include("mergeutils.jl")
+using .MergeUtils
+export @addmeasurement, cstddev, mergereflections
+
end # of module
diff --git a/julia/CrystFEL/src/mergeutils.jl b/julia/CrystFEL/src/mergeutils.jl
new file mode 100644
index 00000000..b75bd2e8
--- /dev/null
+++ b/julia/CrystFEL/src/mergeutils.jl
@@ -0,0 +1,57 @@
+module MergeUtils
+
+using ..CrystFEL.RefLists
+using ..CrystFEL.Symmetry
+export @addmeasurement, cstddev, mergereflections
+
+macro addmeasurement(measurement, weight,
+ mean, sumweight, wksp)
+ return quote
+ delta = $(esc(measurement)) - $(esc(mean))
+ newsumweight = $(esc(sumweight)) + $(esc(weight))
+ R = delta * $(esc(weight)) / newsumweight
+ $(esc(mean)) += R
+ $(esc(wksp)) += $(esc(sumweight)) * delta * R
+ $(esc(sumweight)) = newsumweight
+ end
+end
+
+cstddev(nmeas, work1, work2) = sqrt(work2/work1)/sqrt(nmeas)
+
+
+function mergereflections(correction, crystalrefls, sym)
+
+ merged = RefList{MergedReflection}(sym)
+
+ for (cr,reflections) in crystalrefls
+
+ for refl in reflections
+
+ indices = asymmetricindices(sym, refl.indices)
+ model_version = get!(merged, indices)
+ @addmeasurement(correction(refl, cr), 1.0,
+ model_version.intensity,
+ model_version.temp1,
+ model_version.temp2)
+ model_version.nmeasurements += 1
+
+ end
+ end
+
+ for refl in merged
+ if refl.nmeasurements > 1
+ refl.sigintensity = cstddev(refl.nmeasurements, refl.temp1, refl.temp2)
+ else
+ # Cannot delete a reflection from a list (especially not
+ # while iterating), but setting nmeasurements to zero
+ # prevents it from being written to the output file.
+ refl.nmeasurements = 0
+ end
+ end
+
+ return merged
+
+end
+
+
+end # of module
diff --git a/julia/process_hkl.jl b/julia/process_hkl.jl
index d2cebc86..80b78e52 100644
--- a/julia/process_hkl.jl
+++ b/julia/process_hkl.jl
@@ -1,57 +1,6 @@
using CrystFEL
using LinearAlgebra
-
-macro addmeasurement(measurement, weight,
- mean, sumweight, wksp)
- return quote
- delta = $(esc(measurement)) - $(esc(mean))
- newsumweight = $(esc(sumweight)) + $(esc(weight))
- R = delta * $(esc(weight)) / newsumweight
- $(esc(mean)) += R
- $(esc(wksp)) += $(esc(sumweight)) * delta * R
- $(esc(sumweight)) = newsumweight
- end
-end
-
-cstddev(nmeas, work1, work2) = sqrt(work2/work1)/sqrt(nmeas)
-
-
-function mergereflections(correction, crystalrefls, sym)
-
- merged = RefList{MergedReflection}(sym)
-
- for (cr,reflections) in crystalrefls
-
- for refl in reflections
-
- indices = asymmetricindices(sym, refl.indices)
- model_version = get!(merged, indices)
- @addmeasurement(correction(refl, cr), 1.0,
- model_version.intensity,
- model_version.temp1,
- model_version.temp2)
- model_version.nmeasurements += 1
-
- end
- end
-
- for refl in merged
- if refl.nmeasurements > 1
- refl.sigintensity = cstddev(refl.nmeasurements, refl.temp1, refl.temp2)
- else
- # Cannot delete a reflection from a list (especially not
- # while iterating), but setting nmeasurements to zero
- # prevents it from being written to the output file.
- refl.nmeasurements = 0
- end
- end
-
- return merged
-
-end
-
-
function anglebetween(v1, v2)
let v1n = norm(v1), v2n = norm(v2)
return 2*atan(norm(v1*v2n - v2*v1n),
@@ -59,7 +8,6 @@ function anglebetween(v1, v2)
end
end
-
let st = Stream("/home/twhite/experiments/cxidb-193/short.stream", "r"),
merged = mergereflections(allcrystals(st), SymOpList("mmm")) do refl,crystal