diff options
| author | Thomas White <taw@physics.org> | 2024-02-21 11:02:13 +0100 |
|---|---|---|
| committer | Thomas White <taw@physics.org> | 2024-02-21 11:02:13 +0100 |
| commit | e1d77b345190f43eec2cc53684e0ded6b1939cd7 (patch) | |
| tree | 32f97e29bc15809b855f8006bf04242f1b5d2d8c | |
| parent | 0da7d9d1208f273aeb557adaed3f3e0f6dcdfe5a (diff) | |
julia/process_hkl.jl: Break into routines, add correction function
| -rw-r--r-- | julia/process_hkl.jl | 60 |
1 files changed, 34 insertions, 26 deletions
diff --git a/julia/process_hkl.jl b/julia/process_hkl.jl index ac6ec662..6d0dbd40 100644 --- a/julia/process_hkl.jl +++ b/julia/process_hkl.jl @@ -1,38 +1,46 @@ using CrystFEL -st = Stream("input.stream", "r") -sym = SymOpList("2/m_uab") -merged = RefList{MergedReflection}(sym) -for (cr,reflections) in allcrystals(st) +function mergereflections(correction, crystalrefls, sym) - for refl in reflections + merged = RefList{MergedReflection}(sym) - indices = asymmetricindices(sym, refl.indices) - model_version = get!(merged, indices) + for (cr,reflections) in crystalrefls - w = 1.0 - mean = model_version.intensity - sumweight = model_version.temp1 - M2 = model_version.temp2 - temp = w + sumweight - delta = refl.intensity - mean - R = delta * w / temp - model_version.intensity = mean + R - model_version.temp1 = temp - model_version.temp2 = M2 + sumweight * delta * R - model_version.nmeasurements += 1 + for refl in reflections - end + indices = asymmetricindices(sym, refl.indices) + model_version = get!(merged, indices) -end + w = 1.0 + mean = model_version.intensity + sumweight = model_version.temp1 + M2 = model_version.temp2 + temp = w + sumweight + delta = correction(refl.intensity, cr) - mean + R = delta * w / temp + model_version.intensity = mean + R + model_version.temp1 = temp + model_version.temp2 = M2 + sumweight * delta * R + model_version.nmeasurements += 1 + + end + end -for refl in merged - if refl.nmeasurements > 1 - refl.sigintensity = sqrt(refl.temp2/refl.temp1)/sqrt(refl.nmeasurements) - else - refl.nmeasurements = 0 + for refl in merged + if refl.nmeasurements > 1 + refl.sigintensity = sqrt(refl.temp2/refl.temp1)/sqrt(refl.nmeasurements) + else + refl.nmeasurements = 0 + end end + + return merged + end -savereflist!(merged, "merged.hkl") + +let st = Stream("input.stream", "r") + merged = mergereflections((i,cr)->i, allcrystals(st), SymOpList("2/m_uab")) + savereflist!(merged, "merged.hkl") +end |