diff options
author | Thomas White <taw@physics.org> | 2023-09-18 16:43:58 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2023-09-18 16:43:58 +0200 |
commit | 1ae9a458265df2bae2a4d05089e83390ae30084a (patch) | |
tree | bf68604d451365e6b85766bd9dc2758dac306fd2 /libcrystfel | |
parent | 042f9ca307bda09e35769f9eaa87da44877e34e5 (diff) |
Build EXC_WEIGHT into r_dev
This avoids weird weighting factors everywhere and much confusion.
Since Millepede doesn't have an easy way of weighting measurements (only
via altering the ESD values), treating it as a units conversion seems to
be easier.
Diffstat (limited to 'libcrystfel')
-rw-r--r-- | libcrystfel/src/crystfel-mille.c | 6 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.c | 32 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.h | 4 |
3 files changed, 20 insertions, 22 deletions
diff --git a/libcrystfel/src/crystfel-mille.c b/libcrystfel/src/crystfel-mille.c index 6a80e323..1f19d3dc 100644 --- a/libcrystfel/src/crystfel-mille.c +++ b/libcrystfel/src/crystfel-mille.c @@ -200,8 +200,8 @@ void write_mille(Mille *mille, int n, UnitCell *cell, Minvs[rps[i].peak->pn], 0, 0, 0, &local_gradients_fs[j], &local_gradients_ss[j]); - local_gradients_r[j] = EXC_WEIGHT * r_gradient(rvl[j], rps[i].refl, - cell, image->lambda); + local_gradients_r[j] = r_gradient(rvl[j], rps[i].refl, + cell, image->lambda); } /* Global gradients for each hierarchy level, starting at the @@ -249,7 +249,7 @@ void write_mille(Mille *mille, int n, UnitCell *cell, /* Add excitation error "measurement" (local-only) */ mille_add_measurement(mille, nl, local_gradients_r, - 0, NULL, NULL, r_dev(&rps[i])*EXC_WEIGHT, 1.0); + 0, NULL, NULL, r_dev(&rps[i]), 1.0); } } diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 55ecfa9d..ae638025 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -45,10 +45,14 @@ /** \file predict-refine.h */ +/* Weighting of excitation error term (m^-1) compared to position term (pixels) */ +#define EXC_WEIGHT (0.5e-7) + + double r_dev(struct reflpeak *rp) { /* Excitation error term */ - return get_exerr(rp->refl); + return EXC_WEIGHT * get_exerr(rp->refl); } @@ -93,31 +97,31 @@ double r_gradient(int param, Reflection *refl, UnitCell *cell, double wavelength switch ( param ) { case GPARAM_ASX : - return - hs * sin(phi) * cos(azi); + return - hs * sin(phi) * cos(azi) * EXC_WEIGHT; case GPARAM_BSX : - return - ks * sin(phi) * cos(azi); + return - ks * sin(phi) * cos(azi) * EXC_WEIGHT; case GPARAM_CSX : - return - ls * sin(phi) * cos(azi); + return - ls * sin(phi) * cos(azi) * EXC_WEIGHT; case GPARAM_ASY : - return - hs * sin(phi) * sin(azi); + return - hs * sin(phi) * sin(azi) * EXC_WEIGHT; case GPARAM_BSY : - return - ks * sin(phi) * sin(azi); + return - ks * sin(phi) * sin(azi) * EXC_WEIGHT; case GPARAM_CSY : - return - ls * sin(phi) * sin(azi); + return - ls * sin(phi) * sin(azi) * EXC_WEIGHT; case GPARAM_ASZ : - return - hs * cos(phi); + return - hs * cos(phi) * EXC_WEIGHT; case GPARAM_BSZ : - return - ks * cos(phi); + return - ks * cos(phi) * EXC_WEIGHT; case GPARAM_CSZ : - return - ls * cos(phi); + return - ls * cos(phi) * EXC_WEIGHT; /* Detector movements don't affect excitation error */ case GPARAM_DET_TX : @@ -615,7 +619,6 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, float fs_gradients[num_params]; float ss_gradients[num_params]; float r_gradients[num_params]; - double w; /* Calculate all gradients for this parameter */ for ( k=0; k<num_params; k++ ) { @@ -626,7 +629,6 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, } /* Excitation error terms */ - w = EXC_WEIGHT * rps[i].Ih; for ( k=0; k<num_params; k++ ) { int g; @@ -639,14 +641,14 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, /* Matrix is symmetric */ if ( g > k ) continue; - M_c = w * r_gradients[g] * r_gradients[k]; + M_c = r_gradients[g] * r_gradients[k]; M_curr = gsl_matrix_get(M, k, g); gsl_matrix_set(M, k, g, M_curr + M_c); gsl_matrix_set(M, g, k, M_curr + M_c); } - v_c = w * r_dev(&rps[i]); + v_c = r_dev(&rps[i]); v_c *= -r_gradients[k]; v_curr = gsl_vector_get(v, k); gsl_vector_set(v, k, v_curr + v_c); @@ -768,7 +770,7 @@ static double pred_residual(struct reflpeak *rps, int n, struct detgeom *det, res_fs = 0.0; res_ss = 0.0; for ( i=0; i<n; i++ ) { - res_r += EXC_WEIGHT * rps[i].Ih * pow(r_dev(&rps[i]), 2.0); + res_r += pow(r_dev(&rps[i]), 2.0); res_fs += pow(fs_dev(&rps[i], det), 2.0); res_ss += pow(ss_dev(&rps[i], det), 2.0); } diff --git a/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h index 1e72a048..604799c0 100644 --- a/libcrystfel/src/predict-refine.h +++ b/libcrystfel/src/predict-refine.h @@ -53,10 +53,6 @@ enum gparam { }; -/* Weighting of excitation error term (m^-1) compared to position term (pixels) */ -#define EXC_WEIGHT (0.5e-7) - - #include "crystal.h" #include "crystfel-mille.h" |