aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/geometry.c
diff options
context:
space:
mode:
authorWolfgang Brehm <wolfgang.brehm@gmail.com>2019-08-01 18:14:32 +0200
committerThomas White <taw@physics.org>2019-08-26 14:04:06 +0200
commitcc9f746053a86fc88915f8a45b916ab77487a581 (patch)
tree941cdfec72e2831b96c744e9a8b6e3dbb68fad0b /libcrystfel/src/geometry.c
parent63542ec9e052b147da54a656c9eed6ec35efe431 (diff)
loop prediction for gaussian sum spectrum
Diffstat (limited to 'libcrystfel/src/geometry.c')
-rw-r--r--libcrystfel/src/geometry.c88
1 files changed, 69 insertions, 19 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index e5f33d64..366464fa 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -238,6 +238,13 @@ static double random_partiality(signed int h, signed int k, signed int l,
return p;
}
+static inline double safe_khalf(const double xl,
+ const double yl,
+ const double zl)
+{
+ if (zl>0) return 0.0/0.0;
+ return -(xl*xl+yl*yl+zl*zl) / (2.0*zl);
+}
static Reflection *check_reflection(struct image *image, Crystal *cryst,
signed int h, signed int k, signed int l,
@@ -245,32 +252,74 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
Reflection *updateme)
{
Reflection *refl;
- double R, top;
- double kmin, kmax, k0, knom, k1, khalf;
+ double R;
+ double knom=0, khalf=0, kpred=0;
double dcs, exerr;
+ double partiality = 0, mean_kpred=0, M2_kpred=0;
+ double sumw_k = 0, mean_k = 0, M2_k=0;
+ struct gaussian g;
+ /* this arbitrary value is there to mimic previous behavoir */
+ const double min_partiality = exp(-0.5*1.7*1.7);
/* Don't predict 000 */
if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL;
- /* Calculate the limiting wavelengths, lambda0 and lambda1
- * = 1/k0 and 1/k1 respectively */
R = fabs(crystal_get_profile_radius(cryst));
- top = R*R - xl*xl - yl*yl - zl*zl;
- k0 = top/(2.0*(zl+R));
- khalf = (- xl*xl - yl*yl - zl*zl) / (2.0*zl);
- k1 = top/(2.0*(zl-R));
-
- /* The reflection is excited if any of the reflection is within 2sigma
- * of the nominal * wavelength of the X-ray beam
- * (NB image->bw is full width) */
- kmin = 1.0/(image->lambda + image->lambda*image->bw);
- knom = 1.0/image->lambda;
- kmax = 1.0/(image->lambda - image->lambda*image->bw);
- if ( (updateme == NULL) && ((k1>kmax) || (k0<kmin)) ) return NULL;
+ assert(spectrum_get_num_gaussians(image->spectrum));
+ for ( int i = 0 ; i!=spectrum_get_num_gaussians(image->spectrum) ; ++i ) {
+ g = spectrum_get_gaussian(image->spectrum, i);
+
+ double exerr2 = 0;
+ /* project lattice point onto ewald sphere */
+ double x = xl, y = yl, z = zl + g.kcen;
+ double norm = 1.0/sqrt(x*x+y*y+z*z);
+ x*=norm; y*=norm; z*=norm;
+ z-=1;
+ /* width of ewald sphere in the direction of the projection */
+ const double sigma_proj = sqrt(x*x+y*y+z*z)*g.sigma;
+ mean_variance(g.kcen,g.area,&sumw_k,&mean_k,&M2_k);
+ M2_k+=g.area*g.sigma*g.sigma;
+ const double w0 = 1.0/(R*R);
+ const double w1 = 1.0/(sigma_proj*sigma_proj);
+ x*=g.kcen;y*=g.kcen;z*=g.kcen;
+ /* three because the general case fails in extreme cases */
+ if ( w0 / w1 <= DBL_MIN ) { /* "Laue" corner case */
+ kpred = g.kcen;
+ exerr2 = g.kcen - safe_khalf(xl,yl,zl);
+ exerr2*= exerr2;
+ } else if ( w1 / w0 <= DBL_MIN ) { /* monochromatic corner case */
+ kpred = safe_khalf(xl,yl,zl);
+ exerr2 = g.kcen - kpred;
+ exerr2*= exerr2;
+ } else { /* general case */
+ /* closest point on ewald sphere */
+ const double zlp0 = zl>0?zl:0; /* project zl to 0, bit of a hack... */
+ exerr2 = (x-xl)*(x-xl) + (y-yl)*(y-yl) + (z-zl)*(z-zl);
+ /* weighted average between projected lattice point and ewald sphere */
+ x = ( xl *w0 + x*w1 ) / ( w0 + w1 );
+ y = ( yl *w0 + y*w1 ) / ( w0 + w1 );
+ z = ( zlp0*w0 + z*w1 ) / ( w0 + w1 );
+ kpred = safe_khalf(x,y,z);
+ }
+ const double sigma2 = R*R + sigma_proj*sigma_proj;
+ const double exponent = - 0.5 * exerr2 / sigma2;
+ const double overlap_integral = exponent < -700.0 ? 0.0 :
+ exp( exponent )
+ * sqrt( 2 * M_PI * R * R )
+ / sqrt( 2 * M_PI * sigma2 );
+ mean_variance(kpred,g.area*overlap_integral,&partiality,&mean_kpred,&M2_kpred);
+ }
+ partiality *= sqrt( ( R*R + M2_k/sumw_k) / ( R*R ) ); /* reverting the lorentz factor */
+ if ( (updateme == NULL) && ( partiality < min_partiality ) ) return NULL;
+ kpred = mean_kpred;
/* Calculate excitation error */
+ knom = 1.0/image->lambda;
dcs = distance3d(0.0, 0.0, -knom, xl, yl, zl);
exerr = 1.0/image->lambda - dcs; /* Loss of precision */
+ /* to estimate excitation error from partiality but loosing sign */
+ /* that seems to be a problem */
+ /* exerr = R * sqrt( - 2 * log( partiality ) ) ; */
if ( updateme == NULL ) {
refl = reflection_new(h, k, l);
@@ -284,7 +333,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
if ( (image->det != NULL) && (updateme != NULL) ) {
double fs, ss;
- locate_peak_on_panel(xl, yl, zl, knom,
+ locate_peak_on_panel(xl, yl, zl, kpred,
get_panel(updateme), &fs, &ss);
set_detector_pos(refl, fs, ss);
@@ -296,7 +345,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
double fs, ss; /* Position on detector */
signed int p; /* Panel number */
- p = locate_peak(xl, yl, zl, knom,
+ p = locate_peak(xl, yl, zl, kpred,
image->det, &fs, &ss);
if ( p == -1 ) {
reflection_free(refl);
@@ -307,7 +356,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
}
- set_kpred(refl, knom);
+ khalf = (- xl*xl - yl*yl - zl*zl) / (2.0*zl);
+ set_kpred(refl, kpred);
set_khalf(refl, khalf);
set_exerr(refl, exerr);
set_lorentz(refl, 1.0);