Commit e4a0b2c4 authored by Thomas White's avatar Thomas White
Browse files

Polarisation correction for extracted intensities

parent ecd7ed54
......@@ -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;
......@@ -118,6 +119,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"
......@@ -237,6 +239,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;
......@@ -321,7 +324,7 @@ static void *process_image(void *pargsv)
/* Use original data (temporarily) */
simage->data = data_for_measurement;
output_intensities(simage, image.indexed_cell,
pargs->output_mutex);
pargs->output_mutex, config_unpolar);
simage->data = NULL;
}
......@@ -382,6 +385,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 +421,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}
};
......@@ -570,6 +575,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;
......
......@@ -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 ) {
......
......@@ -41,6 +41,9 @@
* for the purposes of integration. */
#define PEAK_REALLY_CLOSE (10.0)
/* Degree of polarisation of X-ray beam */
#define POL (1.0)
struct reflhit {
signed int h;
signed int k;
......@@ -140,7 +143,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 +158,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;
......@@ -178,7 +183,12 @@ 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;
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);
val = image->data[(x+xp)+image->width*(y+yp)] / (sa*pol);
total += val;
......@@ -305,8 +315,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)
......@@ -370,7 +382,7 @@ void dump_peaks(struct image *image, pthread_mutex_t *mutex)
void output_intensities(struct image *image, UnitCell *cell,
pthread_mutex_t *mutex)
pthread_mutex_t *mutex, int unpolar)
{
int x, y;
double ax, ay, az;
......@@ -485,13 +497,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);
}
......
......@@ -22,6 +22,6 @@
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);
#endif /* PEAKS_H */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment