aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/merge.c45
-rw-r--r--src/merge.h6
-rw-r--r--src/partialator.c5
-rw-r--r--src/post-refinement.c5
-rw-r--r--src/rejection.c11
5 files changed, 39 insertions, 33 deletions
diff --git a/src/merge.c b/src/merge.c
index 787afca3..8017182b 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -187,7 +187,7 @@ static void run_merge_job(void *vwargs, int cookie)
Reflection *f;
signed int h, k, l;
double mean, sumweight, M2, temp, delta, R;
- double corr, res, w;
+ double res, w;
struct reflection_contributions *c;
if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
@@ -217,25 +217,16 @@ static void run_merge_job(void *vwargs, int cookie)
continue;
}
- /* Total (multiplicative) correction factor */
- corr = 1.0/G * exp(B*res*res) * get_lorentz(refl)
- / get_partiality(refl);
- if ( isnan(corr) ) {
- ERROR("NaN corr:\n");
- ERROR("G = %f, B = %e\n", G, B);
- ERROR("res = %e\n", res);
- ERROR("p = %f\n", get_partiality(refl));
- unlock_reflection(f);
- continue;
- }
-
/* Reflections count less the more they have to be scaled up */
w = get_partiality(refl);
/* Running mean and variance calculation */
temp = w + sumweight;
- if (ln_merge) delta = log(get_intensity(refl)*corr) - mean;
- else delta = get_intensity(refl)*corr - mean;
+ if ( ln_merge ) {
+ delta = log(correct_reflection(refl, G, B, res)) - mean;
+ } else {
+ delta = correct_reflection(refl, G, B, res) - mean;
+ }
R = delta * w / temp;
set_intensity(f, mean + R);
set_temp2(f, M2 + sumweight * delta * R);
@@ -352,6 +343,25 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
}
+/* Correct intensity in pattern for scaling and Lorentz factors,
+ * but not partiality nor polarisation */
+double correct_reflection_nopart(Reflection *refl, double osf, double Bfac,
+ double res)
+{
+ double corr = osf * exp(-Bfac*res*res);
+ return (get_intensity(refl) / corr) / get_lorentz(refl);
+}
+
+
+/* Correct intensity in pattern for scaling, partiality and Lorentz factors,
+ * but not polarisation */
+double correct_reflection(Reflection *refl, double osf, double Bfac, double res)
+{
+ double Ipart = correct_reflection_nopart(refl, osf, Bfac, res);
+ return Ipart / get_partiality(refl);
+}
+
+
double residual(Crystal *cr, const RefList *full, int free,
int *pn_used, const char *filename)
{
@@ -368,7 +378,7 @@ double residual(Crystal *cr, const RefList *full, int free,
refl != NULL;
refl = next_refl(refl, iter) )
{
- double w, corr, res;
+ double w, res;
signed int h, k, l;
Reflection *match;
double I_full;
@@ -384,8 +394,7 @@ double residual(Crystal *cr, const RefList *full, int free,
if ( get_redundancy(match) < 2 ) continue;
- corr = get_lorentz(refl) / (G * exp(-B*res*res));
- int1 = get_intensity(refl) * corr;
+ int1 = correct_reflection_nopart(refl, G, B, res);
int2 = get_partiality(refl)*I_full;
w = 1.0;
diff --git a/src/merge.h b/src/merge.h
index 476f8755..2c2cfeda 100644
--- a/src/merge.h
+++ b/src/merge.h
@@ -47,6 +47,12 @@ extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
int min_meas, double push_res, int use_weak,
int ln_merge);
+extern double correct_reflection_nopart(Reflection *refl, double osf,
+ double Bfac, double res);
+
+extern double correct_reflection(Reflection *refl, double osf, double Bfac,
+ double res);
+
extern double residual(Crystal *cr, const RefList *full, int free,
int *pn_used, const char *filename);
diff --git a/src/partialator.c b/src/partialator.c
index 10befba8..8cf1dad9 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -638,7 +638,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
{
signed int h, k, l;
double pobs, pcalc;
- double res, corr, Ipart;
+ double res, Ipart;
Reflection *match;
if ( !get_flag(refl) ) continue; /* Not free-flagged */
@@ -660,8 +660,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
pcalc = get_partiality(refl);
/* Observed partiality */
- corr = G * exp(B*res*res) * get_lorentz(refl);
- Ipart = get_intensity(refl) * corr;
+ Ipart = correct_reflection_nopart(refl, G, B, res);
pobs = Ipart / get_intensity(match);
fprintf(fh, "%5i %4i %4i %4i %e %e %8.3f %8.3f %s\n",
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 0a799b9e..b1af92db 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -559,7 +559,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
refl != NULL;
refl = next_refl(refl, iter) )
{
- double corr, Ipart, Ifull, pobs, pcalc;
+ double Ipart, Ifull, pobs, pcalc;
double res;
signed int h, k, l;
Reflection *match;
@@ -573,8 +573,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
/* Don't calculate pobs if reference reflection is weak */
if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
- corr = G * exp(B*res*res) * get_lorentz(refl);
- Ipart = get_intensity(refl) * corr;
+ Ipart = correct_reflection_nopart(refl, G, B, res);
Ifull = get_intensity(match);
pobs = Ipart / Ifull;
pcalc = get_partiality(refl);
diff --git a/src/rejection.c b/src/rejection.c
index 6ec59151..5d1041fa 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -143,11 +143,7 @@ static int calculate_refl_mean_var(RefList *full)
G = crystal_get_osf(c->contrib_crystals[j]);
B = crystal_get_Bfac(c->contrib_crystals[j]);
- Ii = get_intensity(c->contribs[j]);
-
- /* Total (multiplicative) correction factor */
- Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j])
- / get_partiality(c->contribs[j]);
+ Ii = correct_reflection(c->contribs[j], G, B, res);
Ex += Ii - K;
Ex2 += (Ii - K) * (Ii - K);
@@ -232,16 +228,13 @@ static double calculate_cchalf(RefList *template, RefList *full,
while ( exrefl != NULL ) {
double G, B;
- double Ii = get_intensity(exrefl);
G = crystal_get_osf(exclude);
B = crystal_get_Bfac(exclude);
if ( get_partiality(exrefl) > MIN_PART_MERGE ) {
- /* Total (multiplicative) correction factor */
- Ii *= 1.0/G * exp(B*res*res) * get_lorentz(exrefl)
- / get_partiality(exrefl);
+ double Ii = correct_reflection(exrefl, G, B, res);
/* Remove contribution of this reflection */
Ex -= Ii - K;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-17 11:31:56 +0000