From cc9f746053a86fc88915f8a45b916ab77487a581 Mon Sep 17 00:00:00 2001 From: Wolfgang Brehm Date: Thu, 1 Aug 2019 18:14:32 +0200 Subject: loop prediction for gaussian sum spectrum --- libcrystfel/src/geometry.c | 88 ++++++++++++++++++++++++++++++++++++---------- 1 file changed, 69 insertions(+), 19 deletions(-) (limited to 'libcrystfel/src/geometry.c') 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) || (k0spectrum)); + 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); -- cgit v1.2.3