aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2012-07-02 17:00:20 +0200
committerThomas White <taw@physics.org>2012-07-02 17:00:20 +0200
commitbbf7cea5c7aade69c48ca888170e12537accdd29 (patch)
tree96a8dc4d248c5211af06655594f2c570207116ea /src
parentd7cacb7871a7d49057453b768e420dfb0ab52077 (diff)
compare_hkl: Add --sigma-cutoff option
Diffstat (limited to 'src')
-rw-r--r--src/compare_hkl.c44
1 files changed, 36 insertions, 8 deletions
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index b926c3ae..7e3a1634 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -85,6 +85,7 @@ static void show_help(const char *s)
" Choose from: 'Rsplit', 'R1f' and 'R1i'.\n"
" --rmin=<res> Fix lower resolution limit for --shells (m^-1).\n"
" --rmax=<res> Fix upper resolution limit for --shells (m^-1).\n"
+" --sigma-cutoff=<n> Discard reflections with I/sigma(I) < n.\n"
"\n");
}
@@ -298,7 +299,7 @@ int main(int argc, char *argv[])
SymOpList *sym;
double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
double scale_rintint, scale_r1i, scale_r1, scale_r1fi;
- int ncom;
+ int ncom, nrej;
RefList *list1;
RefList *list2;
RefList *list1_raw;
@@ -312,16 +313,18 @@ int main(int argc, char *argv[])
RefListIterator *iter;
int config_unity = 0;
double scale_for_shells;
+ float sigma_cutoff = -INFINITY;
/* Long options */
const struct option longopts[] = {
{"help", 0, NULL, 'h'},
{"ratio" , 1, NULL, 'o'},
{"symmetry", 1, NULL, 'y'},
- {"shells", 1, NULL, 4},
{"pdb", 1, NULL, 'p'},
- {"rmin", 1, NULL, 2},
- {"rmax", 1, NULL, 3},
+ {"rmin", 1, NULL, 2},
+ {"rmax", 1, NULL, 3},
+ {"shells", 1, NULL, 4},
+ {"sigma-cutoff", 1, NULL, 5},
{0, 0, NULL, 0}
};
@@ -373,6 +376,14 @@ int main(int argc, char *argv[])
config_shells = get_r_shell(optarg);
break;
+ case 5 :
+ if ( sscanf(optarg, "%f", &sigma_cutoff) != 1 ) {
+ ERROR("Invalid value for --sigma-cutoff\n");
+ return 1;
+ }
+ break;
+
+
default :
return 1;
@@ -431,12 +442,14 @@ int main(int argc, char *argv[])
/* Find common reflections and calculate ratio */
ratio = reflist_new();
ncom = 0;
+ nrej = 0;
for ( refl1 = first_refl(list1, &iter);
refl1 != NULL;
refl1 = next_refl(refl1, iter) )
{
signed int h, k, l;
double val1, val2;
+ double esd1, esd2;
Reflection *refl2;
Reflection *tr;
@@ -445,11 +458,21 @@ int main(int argc, char *argv[])
refl2 = find_refl(list2, h, k, l);
if ( refl2 == NULL ) continue;
- ncom++;
-
val1 = get_intensity(refl1);
val2 = get_intensity(refl2);
+ esd1 = get_esd_intensity(refl1);
+ esd2 = get_esd_intensity(refl2);
+
+ if ( (val1 < sigma_cutoff * esd1)
+ || (val2 < sigma_cutoff * esd2) )
+ {
+ nrej++;
+ continue;
+ }
+
+ ncom++;
+
/* Add divided version to 'output' list */
tr = add_refl(ratio, h, k, l);
set_intensity(tr, val1/val2);
@@ -463,8 +486,13 @@ int main(int argc, char *argv[])
gsl_set_error_handler_off();
- STATUS("%i,%i reflections: %i in common\n",
- num_reflections(list1), num_reflections(list2), ncom);
+ STATUS("%i,%i reflections: %i in common where both versions have "
+ "I/sigma(I) >= %f.\n",
+ num_reflections(list1), num_reflections(list2), ncom,
+ sigma_cutoff);
+
+ STATUS("Discarded %i reflections because either or both versions "
+ " had I/sigma(I) < %f\n", nrej, sigma_cutoff);
R1 = stat_r1_ignore(list1, list2, &scale_r1fi, config_unity);
STATUS("R1(F) = %5.4f %% (scale=%5.2e)"