aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/diffraction.c24
-rw-r--r--src/diffraction.h2
-rw-r--r--src/hdf5-file.c27
-rw-r--r--src/image.h16
-rw-r--r--src/indexamajig.c46
-rw-r--r--src/likelihood.c4
-rw-r--r--src/likelihood.h3
-rw-r--r--src/pattern_sim.c2
-rw-r--r--src/peaks.c116
-rw-r--r--src/peaks.h3
-rw-r--r--src/process_hkl.c25
11 files changed, 225 insertions, 43 deletions
diff --git a/src/diffraction.c b/src/diffraction.c
index 1960f786..b3e80e8f 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -240,8 +240,8 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
{
struct rvec q;
float twotheta, r, az;
- float rx = 0.0;
- float ry = 0.0;
+ float rx;
+ float ry;
struct panel *p;
const unsigned int x = xs / sampling;
@@ -267,6 +267,26 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
}
+double get_tt(struct image *image, unsigned int xs, unsigned int ys)
+{
+ float r, rx, ry;
+ struct panel *p;
+
+ const unsigned int x = xs;
+ const unsigned int y = ys; /* Integer part only */
+
+ p = find_panel(&image->det, x, y);
+
+ rx = ((float)xs - p->cx) / p->res;
+ ry = ((float)ys - p->cy) / p->res;
+
+ /* Calculate q-vector for this sub-pixel */
+ r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
+
+ return atan2(r, p->clen);
+}
+
+
void get_diffraction(struct image *image, int na, int nb, int nc,
const double *intensities, const unsigned int *counts,
UnitCell *cell, int do_water, GradientMethod m)
diff --git a/src/diffraction.h b/src/diffraction.h
index 85d06cc2..ed583feb 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -32,4 +32,6 @@ extern void get_diffraction(struct image *image, int na, int nb, int nc,
extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
unsigned int sampling, float *ttp, float k);
+extern double get_tt(struct image *image, unsigned int xs, unsigned int ys);
+
#endif /* DIFFRACTION_H */
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index d76f2bd4..5573267b 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -204,6 +204,24 @@ static double get_wavelength(struct hdfile *f)
}
+static double get_f0(struct hdfile *f)
+{
+ herr_t r;
+ hid_t dh;
+ double f0;
+
+ dh = H5Dopen(f->fh, "/LCLS/f_11_ENRC", H5P_DEFAULT);
+ if ( dh < 0 ) return -1.0;
+
+ r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, &f0);
+ H5Dclose(dh);
+ if ( r < 0 ) return -1.0;
+
+ return f0;
+}
+
+
static void debodge_saturation(struct hdfile *f, struct image *image)
{
hid_t dh, sh;
@@ -321,6 +339,15 @@ int hdf5_read(struct hdfile *f, struct image *image)
image->lambda = ph_en_to_lambda(eV_to_J(2000.0));
}
+ image->f0 = get_f0(f);
+ if ( image->f0 < 0.0 ) {
+ ERROR("Couldn't read incident intensity - using 1.0.\n");
+ image->f0 = 1.0;
+ image->f0_available = 0;
+ } else {
+ image->f0_available = 1;
+ }
+
debodge_saturation(f, image);
return 0;
diff --git a/src/image.h b/src/image.h
index 4ba9e5f9..db1b491e 100644
--- a/src/image.h
+++ b/src/image.h
@@ -67,6 +67,16 @@ struct rvec
};
+struct reflhit {
+ signed int h;
+ signed int k;
+ signed int l;
+ double min_distance;
+ int x;
+ int y;
+};
+
+
/* Structure describing an image */
struct image {
@@ -77,6 +87,8 @@ struct image {
int ncells;
struct detector det;
const char *filename;
+ struct reflhit *hits;
+ int n_hits;
int id;
@@ -85,6 +97,10 @@ struct image {
/* Wavelength must always be given */
double lambda; /* Wavelength in m */
+ /* Incident intensity (if unknown, put 1.0) */
+ double f0;
+ int f0_available;
+
int width;
int height;
diff --git a/src/indexamajig.c b/src/indexamajig.c
index bd3b53e2..6cf9cd15 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -57,6 +57,7 @@ struct process_args
int config_gpu;
int config_simulate;
int config_nomatch;
+ int config_unpolar;
IndexingMethod indm;
const double *intensities;
const unsigned int *counts;
@@ -66,6 +67,7 @@ struct process_args
struct process_result
{
int hit;
+ int peaks_sane;
};
@@ -118,6 +120,7 @@ static void show_help(const char *s)
" --no-match Don't attempt to match the indexed cell to the\n"
" model, just proceed with the one generated by the\n"
" auto-indexing procedure.\n"
+" --unpolarized Don't correct for the polarisation of the X-rays.\n"
"\n\nOptions for greater performance or verbosity:\n\n"
" --verbose Be verbose about indexing.\n"
" --gpu Use the GPU to speed up the simulation.\n"
@@ -190,6 +193,11 @@ static struct image *get_simage(struct image *template, int alternate)
image->features = template->features;
image->filename = template->filename;
image->indexed_cell = template->indexed_cell;
+ image->f0 = template->f0;
+
+ /* Prevent muppetry */
+ image->hits = NULL;
+ image->n_hits = 0;
return image;
}
@@ -237,6 +245,7 @@ static void *process_image(void *pargsv)
int config_gpu = pargs->config_gpu;
int config_simulate = pargs->config_simulate;
int config_nomatch = pargs->config_nomatch;
+ int config_unpolar = pargs->config_unpolar;
IndexingMethod indm = pargs->indm;
const double *intensities = pargs->intensities;
const unsigned int *counts = pargs->counts;
@@ -248,6 +257,8 @@ static void *process_image(void *pargsv)
image.indexed_cell = NULL;
image.id = pargs->id;
image.filename = filename;
+ image.hits = NULL;
+ image.n_hits = 0;
STATUS("Processing '%s'\n", image.filename);
@@ -295,6 +306,12 @@ static void *process_image(void *pargsv)
/* Perform 'fine' peak search */
search_peaks(&image);
+
+ /* Get rid of noise-filtered version at this point
+ * - it was strictly for the purposes of peak detection. */
+ free(image.data);
+ image.data = data_for_measurement;
+
if ( image_feature_count(image.features) < 5 ) goto done;
if ( config_dumpfound ) dump_peaks(&image, pargs->output_mutex);
@@ -312,19 +329,26 @@ static void *process_image(void *pargsv)
}
/* No cell at this point? Then we're done. */
+ result->peaks_sane = 0;
if ( image.indexed_cell == NULL ) goto done;
- simage = get_simage(&image, config_alternate);
+ /* Sanity check */
+ if ( !peak_sanity_check(&image, image.indexed_cell) ) {
+ STATUS("Failed peak sanity check.\n");
+ result->peaks_sane = 0;
+ goto done;
+ } else {
+ result->peaks_sane = 1;
+ }
/* Measure intensities if requested */
if ( config_nearbragg ) {
- /* Use original data (temporarily) */
- simage->data = data_for_measurement;
- output_intensities(simage, image.indexed_cell,
- pargs->output_mutex);
- simage->data = NULL;
+ output_intensities(&image, image.indexed_cell,
+ pargs->output_mutex, config_unpolar);
}
+ simage = get_simage(&image, config_alternate);
+
/* Simulate if requested */
if ( config_simulate ) {
if ( config_gpu ) {
@@ -350,7 +374,7 @@ done:
free(image.data);
free(image.det.panels);
image_feature_list_free(image.features);
- free(data_for_measurement);
+ free(image.hits);
hdfile_close(hdfile);
if ( image.indexed_cell == NULL ) {
@@ -371,6 +395,7 @@ int main(int argc, char *argv[])
char *rval = NULL;
int n_images;
int n_hits;
+ int n_sane;
int config_noindex = 0;
int config_dumpfound = 0;
int config_nearbragg = 0;
@@ -382,6 +407,7 @@ int main(int argc, char *argv[])
int config_gpu = 0;
int config_verbose = 0;
int config_alternate = 0;
+ int config_unpolar = 0;
IndexingMethod indm;
char *indm_str = NULL;
UnitCell *cell;
@@ -417,6 +443,7 @@ int main(int argc, char *argv[])
{"intensities", 1, NULL, 'q'},
{"pdb", 1, NULL, 'p'},
{"prefix", 1, NULL, 'x'},
+ {"unpolarized", 0, &config_unpolar, 1},
{0, 0, NULL, 0}
};
@@ -535,6 +562,7 @@ int main(int argc, char *argv[])
gsl_set_error_handler_off();
n_images = 0;
n_hits = 0;
+ n_sane = 0;
for ( i=0; i<nthreads; i++ ) {
worker_args[i] = malloc(sizeof(struct process_args));
@@ -570,6 +598,7 @@ int main(int argc, char *argv[])
pargs->config_gpu = config_gpu;
pargs->config_simulate = config_simulate;
pargs->config_nomatch = config_nomatch;
+ pargs->config_unpolar = config_unpolar;
pargs->cell = cell;
pargs->indm = indm;
pargs->intensities = intensities;
@@ -616,6 +645,7 @@ int main(int argc, char *argv[])
if ( result != NULL ) {
n_hits += result->hit;
+ n_sane += result->peaks_sane;
free(result);
}
@@ -666,7 +696,7 @@ int main(int argc, char *argv[])
fclose(fh);
STATUS("There were %i images.\n", n_images);
- STATUS("%i hits were found.\n", n_hits);
+ STATUS("%i hits were found, of which %i were sane.\n", n_hits, n_sane);
if ( gctx != NULL ) {
cleanup_gpu(gctx);
diff --git a/src/likelihood.c b/src/likelihood.c
index e8c62010..95f005c7 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -26,13 +26,13 @@ void detwin_intensities(const double *model, double *new_pattern,
void scale_intensities(const double *model, double *new_pattern,
const unsigned int *model_counts,
- unsigned int *new_counts)
+ unsigned int *new_counts, double f0, int f0_valid)
{
double s;
unsigned int i;
s = stat_scale_intensity(model, model_counts, new_pattern, new_counts);
- printf("%f\n", s);
+ if ( f0_valid ) printf("%f %f\n", s, f0);
/* NaN -> abort */
if ( isnan(s) ) return;
diff --git a/src/likelihood.h b/src/likelihood.h
index 9a773829..0a2bb407 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -23,7 +23,8 @@ extern void detwin_intensities(const double *model, double *new_pattern,
extern void scale_intensities(const double *model, double *new_pattern,
const unsigned int *model_counts,
- unsigned int *new_counts);
+ unsigned int *new_counts, double f0,
+ int f0_valid);
#endif /* LIKELIHOOD_H */
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index a0d81a9d..e09ab72a 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -346,7 +346,7 @@ int main(int argc, char *argv[])
record_image(&image, !config_nonoise);
if ( config_nearbragg ) {
- output_intensities(&image, cell, NULL);
+ output_intensities(&image, cell, NULL, 1);
}
if ( config_powder ) {
diff --git a/src/peaks.c b/src/peaks.c
index cf4eacf3..22e12949 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -34,21 +34,15 @@
#define MAX_HITS (1024)
/* How close a peak must be to an indexed position to be considered "close"
- * for the purposes of double hit detection etc. */
+ * for the purposes of double hit detection and sanity checking. */
#define PEAK_CLOSE (30.0)
/* How close a peak must be to an indexed position to be considered "close"
* for the purposes of integration. */
#define PEAK_REALLY_CLOSE (10.0)
-struct reflhit {
- signed int h;
- signed int k;
- signed int l;
- double min_distance;
- int x;
- int y;
-};
+/* Degree of polarisation of X-ray beam */
+#define POL (1.0)
#define PEAK_WINDOW_SIZE (10)
@@ -140,7 +134,8 @@ static void cull_peaks(struct image *image)
static void integrate_peak(struct image *image, int xp, int yp,
- float *xc, float *yc, float *intensity)
+ float *xc, float *yc, float *intensity,
+ int do_polar)
{
signed int x, y;
const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
@@ -154,6 +149,7 @@ static void integrate_peak(struct image *image, int xp, int yp,
struct panel *p;
double val, sa, pix_area, Lsq, dsq, proj_area;
float tt;
+ double phi, pa, pb, pol;
/* Circular mask */
if ( x*x + y*y > lim ) continue;
@@ -168,7 +164,7 @@ static void integrate_peak(struct image *image, int xp, int yp,
Lsq = pow(p->clen, 2.0);
/* Area of pixel as seen from crystal (approximate) */
- get_q(image, x+xp, y+yp, 1, &tt, 1.0 / image->lambda);
+ tt = get_tt(image, x+xp, y+yp);
proj_area = pix_area * cos(tt);
/* Calculate distance from crystal to pixel */
@@ -178,7 +174,16 @@ static void integrate_peak(struct image *image, int xp, int yp,
/* Projected area of pixel divided by distance squared */
sa = 1.0e7 * proj_area / (dsq + Lsq);
- val = image->data[(x+xp)+image->width*(y+yp)] / sa;
+ if ( do_polar ) {
+ phi = atan2(y+yp, x+xp);
+ pa = pow(sin(phi)*sin(tt), 2.0);
+ pb = pow(cos(tt), 2.0);
+ pol = 1.0 - 2.0*POL*(1-pa) + POL*(1.0+pb);
+ } else {
+ pol = 1.0;
+ }
+
+ val = image->data[(x+xp)+image->width*(y+yp)] / (sa*pol);
total += val;
@@ -305,8 +310,10 @@ void search_peaks(struct image *image)
continue;
}
- /* Centroid peak and get better coordinates */
- integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity);
+ /* Centroid peak and get better coordinates.
+ * Don't bother doing polarisation correction, because the
+ * intensity of this peak is only an estimate at this stage. */
+ integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity, 0);
/* It is possible for the centroid to fall outside the image */
if ( (fx < 0.0) || (fx > image->width)
@@ -369,23 +376,17 @@ void dump_peaks(struct image *image, pthread_mutex_t *mutex)
}
-void output_intensities(struct image *image, UnitCell *cell,
- pthread_mutex_t *mutex)
+static int find_projected_peaks(struct image *image, UnitCell *cell)
{
int x, y;
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
- double a, b, c, al, be, ga;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- struct reflhit hits[MAX_HITS];
+ struct reflhit *hits;
int n_hits = 0;
- int i;
- int n_found;
- int n_indclose = 0;
- int n_foundclose = 0;
+
+ hits = malloc(sizeof(struct reflhit)*MAX_HITS);
+ if ( hits == NULL ) return 0;
cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
@@ -446,6 +447,59 @@ void output_intensities(struct image *image, UnitCell *cell,
}
STATUS("Found %i reflections\n", n_hits);
+ image->hits = hits;
+ image->n_hits = n_hits;
+
+ return n_hits;
+}
+
+
+int peak_sanity_check(struct image *image, UnitCell *cell)
+{
+ int i;
+ int n_sane = 0;
+
+ find_projected_peaks(image, cell);
+ if ( image->n_hits == 0 ) return 0; /* Failed sanity check: no peaks */
+
+ for ( i=0; i<image->n_hits; i++ ) {
+
+ double d;
+ int idx;
+ struct imagefeature *f;
+
+ f = image_feature_closest(image->features,
+ image->hits[i].x, image->hits[i].y,
+ &d, &idx);
+ if ( (f != NULL) && (d < PEAK_CLOSE) ) {
+ n_sane++;
+ }
+
+ }
+
+ STATUS("Sanity factor: %f / %f = %f\n", (float)n_sane, (float)image->n_hits,
+ (float)n_sane / (float)image->n_hits);
+ if ( (float)n_sane / (float)image->n_hits < 0.8 ) return 0;
+
+ return 1;
+}
+
+
+void output_intensities(struct image *image, UnitCell *cell,
+ pthread_mutex_t *mutex, int unpolar)
+{
+ int i;
+ int n_found;
+ int n_indclose = 0;
+ int n_foundclose = 0;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double a, b, c, al, be, ga;
+ int n_hits = image->n_hits;
+ struct reflhit *hits = image->hits;
+
+ if ( image->n_hits == 0 ) return;
/* Get exclusive access to the output stream if necessary */
if ( mutex != NULL ) pthread_mutex_lock(mutex);
@@ -470,6 +524,13 @@ void output_intensities(struct image *image, UnitCell *cell,
printf("cstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
csx/1e9, csy/1e9, csz/1e9);
+ if ( image->f0_available ) {
+ printf("f0 = %7.5f (arbitrary gas detector units)\n",
+ image->f0);
+ } else {
+ printf("f0 = invalid\n");
+ }
+
for ( i=0; i<n_hits; i++ ) {
float x, y, intensity;
@@ -485,13 +546,14 @@ void output_intensities(struct image *image, UnitCell *cell,
/* f->intensity was measured on the filtered pattern,
* so instead re-integrate using old coordinates.
* This will produce further revised coordinates. */
- integrate_peak(image, f->x, f->y, &x, &y, &intensity);
+ integrate_peak(image, f->x, f->y, &x, &y, &intensity,
+ !unpolar);
intensity = f->intensity;
} else {
integrate_peak(image, hits[i].x, hits[i].y,
- &x, &y, &intensity);
+ &x, &y, &intensity, !unpolar);
}
diff --git a/src/peaks.h b/src/peaks.h
index a5c771b7..10cf5b0a 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -22,6 +22,7 @@
extern void search_peaks(struct image *image);
extern void dump_peaks(struct image *image, pthread_mutex_t *mutex);
extern void output_intensities(struct image *image, UnitCell *cell,
- pthread_mutex_t *mutex);
+ pthread_mutex_t *mutex, int unpolar);
+extern int peak_sanity_check(struct image *image, UnitCell *cell);
#endif /* PEAKS_H */
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 53638e7b..7c42cd3e 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -245,6 +245,8 @@ int main(int argc, char *argv[])
unsigned int *new_counts = NULL;
unsigned int n_total_patterns;
unsigned int *truecounts = NULL;
+ float f0;
+ int f0_valid;
/* Long options */
const struct option longopts[] = {
@@ -359,6 +361,7 @@ int main(int argc, char *argv[])
STATUS("There are %i patterns to process\n", n_total_patterns);
n_patterns = 0;
+ f0_valid = 0;
do {
char line[1024];
@@ -382,10 +385,17 @@ int main(int argc, char *argv[])
model_counts, new_counts);
}
+ /* Assume a default I0 if we don't have one by now */
+ if ( config_scale && !f0_valid ) {
+ ERROR("No f0 value.\n");
+ f0 = 1.0;
+ }
+
/* Scale if requested */
if ( config_scale ) {
scale_intensities(model, new_pattern,
- model_counts, new_counts);
+ model_counts, new_counts, f0,
+ f0_valid);
}
/* Start of second or later pattern */
@@ -409,6 +419,19 @@ int main(int argc, char *argv[])
progress_bar(n_patterns, n_total_patterns, "Merging");
+ f0_valid = 0;
+
+ }
+
+ if ( strncmp(line, "f0 = ", 5) == 0 ) {
+ r = sscanf(line, "f0 = %f", &f0);
+ if ( r != 1 ) {
+ ERROR("Couldn't understand f0 line.\n");
+ f0 = 1.0;
+ f0_valid = 0;
+ continue;
+ }
+ f0_valid = 1;
}
r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-09 23:54:57 +0000