diff options
Diffstat (limited to 'julia/CrystFEL/src/mergeutils.jl')
| -rw-r--r-- | julia/CrystFEL/src/mergeutils.jl | 57 |
1 files changed, 57 insertions, 0 deletions
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 |