aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/compare_hkl.c86
1 files changed, 77 insertions, 9 deletions
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 687d3a31..47b52419 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -8,6 +8,7 @@
*
* Authors:
* 2010-2012 Thomas White <taw@physics.org>
+ * 2013 Lorenzo Galli <lorenzo.galli@desy.de>
*
* This file is part of CrystFEL.
*
@@ -56,7 +57,8 @@ enum fom
FOM_RSPLIT,
FOM_CC,
FOM_CCSTAR,
- FOM_CCANO
+ FOM_CCANO,
+ FOM_CRDANO
};
@@ -69,6 +71,7 @@ static enum fom get_fom(const char *s)
if ( strcasecmp(s, "cc") == 0 ) return FOM_CC;
if ( strcasecmp(s, "ccstar") == 0 ) return FOM_CCSTAR;
if ( strcasecmp(s, "ccano") == 0 ) return FOM_CCANO;
+ if ( strcasecmp(s, "crdano") == 0 ) return FOM_CRDANO;
ERROR("Unknown figure of merit '%s'.\n", s);
exit(1);
@@ -84,7 +87,8 @@ static void show_help(const char *s)
" -y, --symmetry=<sym> The symmetry of both the input files.\n"
" -p, --pdb=<filename> PDB file to use.\n"
" --fom=<FoM> Calculate this figure of merit Choose from:\n"
-" R1I, R1F, R2, Rsplit, CC, CCstar, CCano.\n"
+" R1I, R1F, R2, Rsplit, CC, CCstar,\n"
+" CCano, CRDano.\n"
" --nshells=<n> Use <n> resolution shells.\n"
" -u Force scale factor to 1.\n"
" --shell-file=<file> Write resolution shells to <file>.\n"
@@ -149,6 +153,7 @@ static struct fom_context *init_fom(enum fom fom, int nmax, int nshells)
case FOM_CC :
case FOM_CCSTAR :
case FOM_CCANO :
+ case FOM_CRDANO :
fctx->vec1 = malloc(nshells*sizeof(double *));
fctx->vec2 = malloc(nshells*sizeof(double *));
if ( (fctx->vec1 == NULL) || (fctx->vec2 == NULL) ) return NULL;
@@ -213,6 +218,7 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2,
break;
case FOM_CCANO :
+ case FOM_CRDANO :
assert(fctx->n[bin] < fctx->nmax);
fctx->vec1[bin][fctx->n[bin]] = i1 - i1bij;
fctx->vec2[bin][fctx->n[bin]] = i2 - i2bij;
@@ -232,6 +238,10 @@ static double fom_overall(struct fom_context *fctx)
double *overall_vec1;
double *overall_vec2;
int overall_n;
+ double *overall_along_diagonal;
+ double *overall_perpend_diagonal;
+ double variance_signal;
+ double variance_error;
double cc = INFINITY;
switch ( fctx->fom ) {
@@ -268,6 +278,32 @@ static double fom_overall(struct fom_context *fctx)
free(overall_vec2);
break;
+ case FOM_CRDANO :
+ overall_along_diagonal = malloc(fctx->nmax*sizeof(double));
+ overall_perpend_diagonal = malloc(fctx->nmax*sizeof(double));
+ overall_n = 0;
+ for ( i=0; i<fctx->nshells; i++ ) {
+ int j;
+ for ( j=0; j<fctx->n[i]; j++ ) {
+ overall_along_diagonal[overall_n] =
+ ( fctx->vec1[i][j] + fctx->vec2[i][j] )
+ / sqrt(2.0);
+ overall_perpend_diagonal[overall_n] =
+ ( fctx->vec1[i][j] - fctx->vec2[i][j] )
+ / sqrt(2.0);
+ overall_n++;
+ }
+ }
+ variance_signal = gsl_stats_variance(overall_along_diagonal, 1,
+ overall_n);
+ variance_error = gsl_stats_variance(overall_perpend_diagonal, 1,
+ overall_n);
+ cc = sqrt(variance_signal / variance_error );
+
+ free(overall_along_diagonal);
+ free(overall_perpend_diagonal);
+ break;
+
}
switch ( fctx->fom ) {
@@ -284,6 +320,7 @@ static double fom_overall(struct fom_context *fctx)
case FOM_CC :
case FOM_CCANO :
+ case FOM_CRDANO :
return cc;
case FOM_CCSTAR :
@@ -299,6 +336,11 @@ static double fom_overall(struct fom_context *fctx)
static double fom_shell(struct fom_context *fctx, int i)
{
double cc;
+ int j;
+ double variance_signal;
+ double variance_error;
+ double along_diagonal[fctx->n[i]];
+ double perpend_diagonal[fctx->n[i]];
switch ( fctx->fom ) {
@@ -322,6 +364,20 @@ static double fom_shell(struct fom_context *fctx, int i)
fctx->n[i]);
return sqrt((2.0*cc)/(1.0+cc));
+
+ case FOM_CRDANO :
+ for ( j=0; j<fctx->n[i]; j++ ) {
+ along_diagonal[j] =( fctx->vec1[i][j] +
+ fctx->vec2[i][j] ) / sqrt(2.0);
+ perpend_diagonal[j] =( fctx->vec1[i][j] -
+ fctx->vec2[i][j] ) / sqrt(2.0);
+ }
+ variance_signal = gsl_stats_variance(along_diagonal, 1,
+ fctx->n[i]);
+ variance_error = gsl_stats_variance(perpend_diagonal, 1,
+ fctx->n[i]);
+ return sqrt(variance_signal / variance_error);
+
}
ERROR("This point is never reached.\n");
@@ -558,7 +614,7 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
i1 = get_intensity(refl1);
i2 = scale * get_intensity(refl2);
- if ( fom == FOM_CCANO ) {
+ if ( fom == FOM_CCANO || fom == FOM_CRDANO ) {
Reflection *refl1_bij = NULL;
Reflection *refl2_bij = NULL;
@@ -627,6 +683,10 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
STATUS("Overall CCano = %.7f\n", fom_overall(fctx));
break;
+ case FOM_CRDANO :
+ STATUS("Overall CRDano = %.7f\n", fom_overall(fctx));
+ break;
+
}
fh = fopen(filename, "w");
@@ -673,6 +733,10 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
fprintf(fh, "%s CCano nref%s\n", t1, t2);
break;
+ case FOM_CRDANO :
+ fprintf(fh, "%s CRDano nref%s\n", t1, t2);
+ break;
+
}
for ( i=0; i<nshells; i++ ) {
@@ -700,6 +764,7 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
case FOM_CC :
case FOM_CCSTAR :
case FOM_CCANO :
+ case FOM_CRDANO :
if ( config_intshells ) {
fprintf(fh, "%10.3f %10.7f %10i\n",
cen, r, cts[i]);
@@ -862,6 +927,7 @@ int main(int argc, char *argv[])
case FOM_CC :
case FOM_CCSTAR :
case FOM_CCANO :
+ case FOM_CRDANO :
break;
}
}
@@ -872,14 +938,16 @@ int main(int argc, char *argv[])
sym = get_pointgroup(sym_str);
free(sym_str);
- if ( (fom == FOM_CCANO) && is_centrosymmetric(sym) ) {
+ if ( (fom == FOM_CCANO || fom == FOM_CRDANO ) &&
+ is_centrosymmetric(sym) ) {
- /* Refuse to calculate CCano with a centrosymmetric point group.
+ /* Refuse to calculate CCano or CRDano with a
+ * centrosymmetric point group.
* It'll "work", in the sense of producing a number, but the
* number won't have anything to do with anomalous signal. */
- ERROR("You are trying to calculate CCano with a centrosymmetric"
- " point group.\n");
+ ERROR("You are trying to calculate CCano or CRDano with a"
+ "centrosymmetric point group.\n");
ERROR("This is a silly thing to do, and I'm refusing to help"
" you do it.\n");
ERROR("Please review your earlier processing steps and try"
@@ -1014,7 +1082,7 @@ int main(int argc, char *argv[])
}
}
- if ( fom == FOM_CCANO ) {
+ if ( fom == FOM_CCANO || fom == FOM_CRDANO ) {
Reflection *refl1_bij = NULL;
Reflection *refl2_bij = NULL;
@@ -1083,7 +1151,7 @@ int main(int argc, char *argv[])
if ( nbij > 0 ) {
STATUS("%i reflection pairs rejected because either or both"
- " versions did not have Bijvoet partners.\n", nbij);
+ " versions did not have Bijvoet partners.\n", nres);
}
if ( ncen > 0 ) {
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-07 15:43:14 +0000