aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--doc/man/compare_hkl.15
-rw-r--r--doc/man/process_hkl.14
-rw-r--r--libcrystfel/src/detector.h2
-rw-r--r--libcrystfel/src/reflist-utils.c4
-rw-r--r--src/check_hkl.c6
-rw-r--r--src/compare_hkl.c73
-rw-r--r--src/partial_sim.c9
-rw-r--r--src/process_hkl.c44
8 files changed, 101 insertions, 46 deletions
diff --git a/doc/man/compare_hkl.1 b/doc/man/compare_hkl.1
index a02cbd9c..b44f16df 100644
--- a/doc/man/compare_hkl.1
+++ b/doc/man/compare_hkl.1
@@ -42,6 +42,11 @@ Fix the lower resolution limit for the resolutions shells, as 1/d in m^-1. Refl
Fix the upper resolution limit for the resolutions shells, as 1/d in m^-1 Reflections outside the specified resolution range will still be included in the calculation of overall figures of merit.
.PD 0
+.IP \fB--sigma-cutoff=\fR\fIn\fR
+.PD
+Discard reflections with I/sigma(I) < \fIn\fR. Default: -infinity (no cutoff).
+
+.PD 0
.IP \fB--shells=\fR\fIFoM\fR
.PD
Calculate \fIFoM\fR in resolution shells. Possible figures of merit for \fIFoM\fR are:
diff --git a/doc/man/process_hkl.1 b/doc/man/process_hkl.1
index db191d4b..76c7505e 100644
--- a/doc/man/process_hkl.1
+++ b/doc/man/process_hkl.1
@@ -99,6 +99,10 @@ calculating 'powder' patterns.
Output the maximum intensity for each reflection rather than the mean value of
all intensities. This is usually not useful.
+.PD 0
+.IP \fB--no-polarisation\fR
+.PD
+Disable the polarisation correction.
.SH CHOICE OF POINT GROUP FOR MERGING
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index b740965c..1cd64716 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -149,6 +149,8 @@ extern double smallest_q(struct image *image);
extern struct panel *find_panel_by_name(struct detector *det, const char *name);
+extern int write_detector_geometry(const char *filename, struct detector *det);
+
#ifdef __cplusplus
}
#endif
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index 1404b545..2e6715f2 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -292,9 +292,7 @@ RefList *read_reflections_from_file(FILE *fh)
set_redundancy(refl, cts);
ph = strtod(phs, &v);
- if ( v != NULL ) set_phase(refl, deg2rad(ph));
-
- /* The 1/d value is actually ignored. */
+ if ( v != phs ) set_phase(refl, deg2rad(ph));
}
diff --git a/src/check_hkl.c b/src/check_hkl.c
index 4927e60e..32f41cfd 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -282,13 +282,13 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
fprintf(fh, "1/d centre # refs Possible Compl "
- "Meas Red SNR Std dev Mean\n");
+ "Meas Red SNR Std dev Mean d(A)\n");
for ( i=0; i<NBINS; i++ ) {
double cen;
cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
fprintf(fh, "%10.3f %8i %8i %6.2f %10i %5.1f"
- " %5.2f %10.2f %10.2f\n",
+ " %5.2f %10.2f %10.2f %8.2f\n",
cen*1.0e-9,
measured[i],
possible[i],
@@ -297,7 +297,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
(double)measurements[i]/measured[i],
snr[i]/(double)measured[i],
sqrt(var[i]/measured[i]),
- mean[i]);
+ mean[i], (1.0/cen)*1e10);
}
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index b926c3ae..2f846087 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");
}
@@ -106,7 +107,7 @@ static void plot_shells(RefList *list1, RefList *list2, double scale,
Reflection *refl1;
RefListIterator *iter;
FILE *fh;
- double den;
+ double den[NBINS];
int ctot, nout;
if ( cell == NULL ) {
@@ -116,6 +117,7 @@ static void plot_shells(RefList *list1, RefList *list2, double scale,
for ( i=0; i<NBINS; i++ ) {
num[i] = 0.0;
+ den[i] = 0.0;
cts[i] = 0;
measured[i] = 0;
measurements[i] = 0;
@@ -161,7 +163,7 @@ static void plot_shells(RefList *list1, RefList *list2, double scale,
STATUS("Shell %i: %f to %f\n", NBINS-1,
rmins[NBINS-1]/1e9, rmaxs[NBINS-1]/1e9);
- den = 0.0; ctot = 0; nout = 0;
+ ctot = 0; nout = 0;
for ( refl1 = first_refl(list1, &iter);
refl1 != NULL;
refl1 = next_refl(refl1, iter) )
@@ -201,18 +203,18 @@ static void plot_shells(RefList *list1, RefList *list2, double scale,
switch ( config_shells ) {
case R_SHELL_RSPLIT :
- num[bin] += fabs(i1 - scale*i2);
- den += i1 + scale*i2;
+ num[bin] += fabs(i1 - i2);
+ den[bin] += i1 + i2;
break;
case R_SHELL_R1I :
num[bin] += fabs(i1 - scale*i2);
- den += i1;
+ den[bin] += i1;
break;
case R_SHELL_R1F :
num[bin] += fabs(f1 - scale*f2);
- den += f1;
+ den[bin] += f1;
break;
default : break;
@@ -237,19 +239,19 @@ static void plot_shells(RefList *list1, RefList *list2, double scale,
switch ( config_shells ) {
case R_SHELL_RSPLIT :
- fprintf(fh, "1/d centre Rsplit / %%\n");
+ fprintf(fh, "1/d centre Rsplit / %% nref d (A)\n");
break;
case R_SHELL_R1I :
- fprintf(fh, "1/d centre R1(I) / %%\n");
+ fprintf(fh, "1/d centre R1(I) / %% nref d (A)\n");
break;
case R_SHELL_R1F :
- fprintf(fh, "1/d centre R1(F) ignoring -ves / %%\n");
+ fprintf(fh, "1/d centre R1(F) ign -/%% nref d (A)\n");
break;
default :
- fprintf(fh, "1/d centre 0.0\n");
+ fprintf(fh, "1/d centre 0.0 nref d (A)\n");
break;
}
@@ -262,12 +264,12 @@ static void plot_shells(RefList *list1, RefList *list2, double scale,
switch ( config_shells ) {
case R_SHELL_RSPLIT :
- r = (2.0*(num[i]/den)*((double)ctot/cts[i]))/sqrt(2.0);
+ r = 2.0*(num[i]/den[i]) / sqrt(2.0);
break;
case R_SHELL_R1I :
case R_SHELL_R1F :
- r = (num[i]/den) * ((double)ctot/cts[i]);
+ r = num[i]/den[i];
break;
default :
@@ -276,8 +278,8 @@ static void plot_shells(RefList *list1, RefList *list2, double scale,
}
- fprintf(fh, "%10.3f %10.2f %10i\n",
- cen*1.0e-9, r*100.0, cts[i]);
+ fprintf(fh, "%10.3f %10.2f %10i %10.2f\n",
+ cen*1.0e-9, r*100.0, cts[i], (1.0/cen)*1e10);
}
@@ -298,7 +300,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 +314,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 +377,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 +443,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 +459,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 +487,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)"
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 9e054675..3cfb94a2 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -189,6 +189,7 @@ struct queue_args
pthread_mutex_t full_lock;
int n_done;
+ int n_started;
int n_to_do;
SymOpList *sym;
@@ -225,11 +226,16 @@ static void *create_job(void *vqargs)
struct worker_args *wargs;
struct queue_args *qargs = vqargs;
+ /* All done already? */
+ if ( qargs->n_started == qargs->n_to_do ) return NULL;
+
wargs = malloc(sizeof(struct worker_args));
wargs->qargs = qargs;
wargs->image = *qargs->template_image;
+ qargs->n_started++;
+
return wargs;
}
@@ -482,7 +488,7 @@ int main(int argc, char *argv[])
}
if ( output_file == NULL ) {
- ERROR("You must pgive a filename for the output.\n");
+ ERROR("You must give a filename for the output.\n");
return 1;
}
ofh = fopen(output_file, "w");
@@ -512,6 +518,7 @@ int main(int argc, char *argv[])
pthread_mutex_init(&qargs.full_lock, NULL);
qargs.n_to_do = n;
qargs.n_done = 0;
+ qargs.n_started = 0;
qargs.sym = sym;
qargs.random_intensities = random_intensities;
qargs.cell = cell;
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 26abfd4f..f9e9af99 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -73,6 +73,7 @@ static void show_help(const char *s)
" model.\n"
" --reference=<file> Compare against intensities from <file> when\n"
" scaling. \n"
+" --no-polarisation Disable polarisation correction.\n"
);
}
@@ -170,7 +171,8 @@ static double scale_intensities(RefList *reference, RefList *new,
static int merge_pattern(RefList *model, struct image *new, RefList *reference,
const SymOpList *sym,
double *hist_vals, signed int hist_h,
- signed int hist_k, signed int hist_l, int *hist_n)
+ signed int hist_k, signed int hist_l, int *hist_n,
+ int config_nopolar)
{
Reflection *refl;
RefListIterator *iter;
@@ -214,19 +216,23 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
intensity = scale * get_intensity(refl);
- /* Polarisation correction assuming 100% polarisation along the
- * x direction */
- xl = h*asx + k*bsx + l*csx;
- yl = h*asy + k*bsy + l*csy;
- zl = h*asz + k*bsz + l*csz;
+ if ( !config_nopolar ) {
- ool = 1.0 / new->lambda;
- tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool);
- phi = atan2(yl, xl);
- pa = pow(sin(phi)*sin(tt), 2.0);
- pb = pow(cos(tt), 2.0);
- pol = 1.0 - 2.0*(1.0-pa) + (1.0+pb);
- intensity /= pol;
+ /* Polarisation correction assuming 100% polarisation
+ * along the x direction */
+ xl = h*asx + k*bsx + l*csx;
+ yl = h*asy + k*bsy + l*csy;
+ zl = h*asz + k*bsz + l*csz;
+
+ ool = 1.0 / new->lambda;
+ tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool);
+ phi = atan2(yl, xl);
+ pa = pow(sin(phi)*sin(tt), 2.0);
+ pb = pow(cos(tt), 2.0);
+ pol = 1.0 - 2.0*(1.0-pa) + (1.0+pb);
+ intensity /= pol;
+
+ }
cur_intensity = get_intensity(model_version);
set_intensity(model_version, cur_intensity + intensity);
@@ -258,7 +264,7 @@ static void merge_all(FILE *fh, RefList *model, RefList *reference,
int n_total_patterns,
double *hist_vals, signed int hist_h,
signed int hist_k, signed int hist_l,
- int *hist_i)
+ int *hist_i, int config_nopolar)
{
int rval;
int n_patterns = 0;
@@ -289,7 +295,7 @@ static void merge_all(FILE *fh, RefList *model, RefList *reference,
r = merge_pattern(model, &image, reference, sym,
hist_vals, hist_h, hist_k, hist_l,
- hist_i);
+ hist_i, config_nopolar);
if ( r == 0 ) n_used++;
@@ -355,6 +361,7 @@ int main(int argc, char *argv[])
int hist_i;
int space_for_hist = 0;
char *histo_params = NULL;
+ int config_nopolar = 0;
/* Long options */
const struct option longopts[] = {
@@ -367,6 +374,8 @@ int main(int argc, char *argv[])
{"start-after", 1, NULL, 'f'},
{"sum", 0, &config_sum, 1},
{"scale", 0, &config_scale, 1},
+ {"no-polarisation", 0, &config_nopolar, 1},
+ {"no-polarization", 0, &config_nopolar, 1},
{"symmetry", 1, NULL, 'y'},
{"histogram", 1, NULL, 'g'},
{"hist-parameters", 1, NULL, 'z'},
@@ -503,7 +512,7 @@ int main(int argc, char *argv[])
hist_i = 0;
merge_all(fh, model, NULL, config_startafter, config_stopafter,
sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l,
- &hist_i);
+ &hist_i, config_nopolar);
if ( ferror(fh) ) {
ERROR("Stream read error.\n");
return 1;
@@ -526,7 +535,8 @@ int main(int argc, char *argv[])
merge_all(fh, model, reference,
config_startafter, config_stopafter, sym,
n_total_patterns,
- hist_vals, hist_h, hist_k, hist_l, &hist_i);
+ hist_vals, hist_h, hist_k, hist_l, &hist_i,
+ config_nopolar);
if ( ferror(fh) ) {
ERROR("Stream read error.\n");
generated by cgit v1.2.3 (git 2.25.1) at 2024-08-13 21:35:52 +0000