diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/check_hkl.c | 6 | ||||
-rw-r--r-- | src/compare_hkl.c | 73 | ||||
-rw-r--r-- | src/partial_sim.c | 9 | ||||
-rw-r--r-- | src/process_hkl.c | 44 |
4 files changed, 89 insertions, 43 deletions
diff --git a/src/check_hkl.c b/src/check_hkl.c index 4927e60e..32f41cfd 100644 --- a/src/check_hkl.c +++ b/src/check_hkl.c @@ -282,13 +282,13 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, fprintf(fh, "1/d centre # refs Possible Compl " - "Meas Red SNR Std dev Mean\n"); + "Meas Red SNR Std dev Mean d(A)\n"); for ( i=0; i<NBINS; i++ ) { double cen; cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0; fprintf(fh, "%10.3f %8i %8i %6.2f %10i %5.1f" - " %5.2f %10.2f %10.2f\n", + " %5.2f %10.2f %10.2f %8.2f\n", cen*1.0e-9, measured[i], possible[i], @@ -297,7 +297,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, (double)measurements[i]/measured[i], snr[i]/(double)measured[i], sqrt(var[i]/measured[i]), - mean[i]); + mean[i], (1.0/cen)*1e10); } diff --git a/src/compare_hkl.c b/src/compare_hkl.c index b926c3ae..2f846087 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -85,6 +85,7 @@ static void show_help(const char *s) " Choose from: 'Rsplit', 'R1f' and 'R1i'.\n" " --rmin=<res> Fix lower resolution limit for --shells (m^-1).\n" " --rmax=<res> Fix upper resolution limit for --shells (m^-1).\n" +" --sigma-cutoff=<n> Discard reflections with I/sigma(I) < n.\n" "\n"); } @@ -106,7 +107,7 @@ static void plot_shells(RefList *list1, RefList *list2, double scale, Reflection *refl1; RefListIterator *iter; FILE *fh; - double den; + double den[NBINS]; int ctot, nout; if ( cell == NULL ) { @@ -116,6 +117,7 @@ static void plot_shells(RefList *list1, RefList *list2, double scale, for ( i=0; i<NBINS; i++ ) { num[i] = 0.0; + den[i] = 0.0; cts[i] = 0; measured[i] = 0; measurements[i] = 0; @@ -161,7 +163,7 @@ static void plot_shells(RefList *list1, RefList *list2, double scale, STATUS("Shell %i: %f to %f\n", NBINS-1, rmins[NBINS-1]/1e9, rmaxs[NBINS-1]/1e9); - den = 0.0; ctot = 0; nout = 0; + ctot = 0; nout = 0; for ( refl1 = first_refl(list1, &iter); refl1 != NULL; refl1 = next_refl(refl1, iter) ) @@ -201,18 +203,18 @@ static void plot_shells(RefList *list1, RefList *list2, double scale, switch ( config_shells ) { case R_SHELL_RSPLIT : - num[bin] += fabs(i1 - scale*i2); - den += i1 + scale*i2; + num[bin] += fabs(i1 - i2); + den[bin] += i1 + i2; break; case R_SHELL_R1I : num[bin] += fabs(i1 - scale*i2); - den += i1; + den[bin] += i1; break; case R_SHELL_R1F : num[bin] += fabs(f1 - scale*f2); - den += f1; + den[bin] += f1; break; default : break; @@ -237,19 +239,19 @@ static void plot_shells(RefList *list1, RefList *list2, double scale, switch ( config_shells ) { case R_SHELL_RSPLIT : - fprintf(fh, "1/d centre Rsplit / %%\n"); + fprintf(fh, "1/d centre Rsplit / %% nref d (A)\n"); break; case R_SHELL_R1I : - fprintf(fh, "1/d centre R1(I) / %%\n"); + fprintf(fh, "1/d centre R1(I) / %% nref d (A)\n"); break; case R_SHELL_R1F : - fprintf(fh, "1/d centre R1(F) ignoring -ves / %%\n"); + fprintf(fh, "1/d centre R1(F) ign -/%% nref d (A)\n"); break; default : - fprintf(fh, "1/d centre 0.0\n"); + fprintf(fh, "1/d centre 0.0 nref d (A)\n"); break; } @@ -262,12 +264,12 @@ static void plot_shells(RefList *list1, RefList *list2, double scale, switch ( config_shells ) { case R_SHELL_RSPLIT : - r = (2.0*(num[i]/den)*((double)ctot/cts[i]))/sqrt(2.0); + r = 2.0*(num[i]/den[i]) / sqrt(2.0); break; case R_SHELL_R1I : case R_SHELL_R1F : - r = (num[i]/den) * ((double)ctot/cts[i]); + r = num[i]/den[i]; break; default : @@ -276,8 +278,8 @@ static void plot_shells(RefList *list1, RefList *list2, double scale, } - fprintf(fh, "%10.3f %10.2f %10i\n", - cen*1.0e-9, r*100.0, cts[i]); + fprintf(fh, "%10.3f %10.2f %10i %10.2f\n", + cen*1.0e-9, r*100.0, cts[i], (1.0/cen)*1e10); } @@ -298,7 +300,7 @@ int main(int argc, char *argv[]) SymOpList *sym; double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson; double scale_rintint, scale_r1i, scale_r1, scale_r1fi; - int ncom; + int ncom, nrej; RefList *list1; RefList *list2; RefList *list1_raw; @@ -312,16 +314,18 @@ int main(int argc, char *argv[]) RefListIterator *iter; int config_unity = 0; double scale_for_shells; + float sigma_cutoff = -INFINITY; /* Long options */ const struct option longopts[] = { {"help", 0, NULL, 'h'}, {"ratio" , 1, NULL, 'o'}, {"symmetry", 1, NULL, 'y'}, - {"shells", 1, NULL, 4}, {"pdb", 1, NULL, 'p'}, - {"rmin", 1, NULL, 2}, - {"rmax", 1, NULL, 3}, + {"rmin", 1, NULL, 2}, + {"rmax", 1, NULL, 3}, + {"shells", 1, NULL, 4}, + {"sigma-cutoff", 1, NULL, 5}, {0, 0, NULL, 0} }; @@ -373,6 +377,14 @@ int main(int argc, char *argv[]) config_shells = get_r_shell(optarg); break; + case 5 : + if ( sscanf(optarg, "%f", &sigma_cutoff) != 1 ) { + ERROR("Invalid value for --sigma-cutoff\n"); + return 1; + } + break; + + default : return 1; @@ -431,12 +443,14 @@ int main(int argc, char *argv[]) /* Find common reflections and calculate ratio */ ratio = reflist_new(); ncom = 0; + nrej = 0; for ( refl1 = first_refl(list1, &iter); refl1 != NULL; refl1 = next_refl(refl1, iter) ) { signed int h, k, l; double val1, val2; + double esd1, esd2; Reflection *refl2; Reflection *tr; @@ -445,11 +459,21 @@ int main(int argc, char *argv[]) refl2 = find_refl(list2, h, k, l); if ( refl2 == NULL ) continue; - ncom++; - val1 = get_intensity(refl1); val2 = get_intensity(refl2); + esd1 = get_esd_intensity(refl1); + esd2 = get_esd_intensity(refl2); + + if ( (val1 < sigma_cutoff * esd1) + || (val2 < sigma_cutoff * esd2) ) + { + nrej++; + continue; + } + + ncom++; + /* Add divided version to 'output' list */ tr = add_refl(ratio, h, k, l); set_intensity(tr, val1/val2); @@ -463,8 +487,13 @@ int main(int argc, char *argv[]) gsl_set_error_handler_off(); - STATUS("%i,%i reflections: %i in common\n", - num_reflections(list1), num_reflections(list2), ncom); + STATUS("%i,%i reflections: %i in common where both versions have " + "I/sigma(I) >= %f.\n", + num_reflections(list1), num_reflections(list2), ncom, + sigma_cutoff); + + STATUS("Discarded %i reflections because either or both versions " + " had I/sigma(I) < %f\n", nrej, sigma_cutoff); R1 = stat_r1_ignore(list1, list2, &scale_r1fi, config_unity); STATUS("R1(F) = %5.4f %% (scale=%5.2e)" diff --git a/src/partial_sim.c b/src/partial_sim.c index 9e054675..3cfb94a2 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -189,6 +189,7 @@ struct queue_args pthread_mutex_t full_lock; int n_done; + int n_started; int n_to_do; SymOpList *sym; @@ -225,11 +226,16 @@ static void *create_job(void *vqargs) struct worker_args *wargs; struct queue_args *qargs = vqargs; + /* All done already? */ + if ( qargs->n_started == qargs->n_to_do ) return NULL; + wargs = malloc(sizeof(struct worker_args)); wargs->qargs = qargs; wargs->image = *qargs->template_image; + qargs->n_started++; + return wargs; } @@ -482,7 +488,7 @@ int main(int argc, char *argv[]) } if ( output_file == NULL ) { - ERROR("You must pgive a filename for the output.\n"); + ERROR("You must give a filename for the output.\n"); return 1; } ofh = fopen(output_file, "w"); @@ -512,6 +518,7 @@ int main(int argc, char *argv[]) pthread_mutex_init(&qargs.full_lock, NULL); qargs.n_to_do = n; qargs.n_done = 0; + qargs.n_started = 0; qargs.sym = sym; qargs.random_intensities = random_intensities; qargs.cell = cell; diff --git a/src/process_hkl.c b/src/process_hkl.c index 26abfd4f..f9e9af99 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -73,6 +73,7 @@ static void show_help(const char *s) " model.\n" " --reference=<file> Compare against intensities from <file> when\n" " scaling. \n" +" --no-polarisation Disable polarisation correction.\n" ); } @@ -170,7 +171,8 @@ static double scale_intensities(RefList *reference, RefList *new, static int merge_pattern(RefList *model, struct image *new, RefList *reference, const SymOpList *sym, double *hist_vals, signed int hist_h, - signed int hist_k, signed int hist_l, int *hist_n) + signed int hist_k, signed int hist_l, int *hist_n, + int config_nopolar) { Reflection *refl; RefListIterator *iter; @@ -214,19 +216,23 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, intensity = scale * get_intensity(refl); - /* Polarisation correction assuming 100% polarisation along the - * x direction */ - xl = h*asx + k*bsx + l*csx; - yl = h*asy + k*bsy + l*csy; - zl = h*asz + k*bsz + l*csz; + if ( !config_nopolar ) { - ool = 1.0 / new->lambda; - tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool); - phi = atan2(yl, xl); - pa = pow(sin(phi)*sin(tt), 2.0); - pb = pow(cos(tt), 2.0); - pol = 1.0 - 2.0*(1.0-pa) + (1.0+pb); - intensity /= pol; + /* Polarisation correction assuming 100% polarisation + * along the x direction */ + xl = h*asx + k*bsx + l*csx; + yl = h*asy + k*bsy + l*csy; + zl = h*asz + k*bsz + l*csz; + + ool = 1.0 / new->lambda; + tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool); + phi = atan2(yl, xl); + pa = pow(sin(phi)*sin(tt), 2.0); + pb = pow(cos(tt), 2.0); + pol = 1.0 - 2.0*(1.0-pa) + (1.0+pb); + intensity /= pol; + + } cur_intensity = get_intensity(model_version); set_intensity(model_version, cur_intensity + intensity); @@ -258,7 +264,7 @@ static void merge_all(FILE *fh, RefList *model, RefList *reference, int n_total_patterns, double *hist_vals, signed int hist_h, signed int hist_k, signed int hist_l, - int *hist_i) + int *hist_i, int config_nopolar) { int rval; int n_patterns = 0; @@ -289,7 +295,7 @@ static void merge_all(FILE *fh, RefList *model, RefList *reference, r = merge_pattern(model, &image, reference, sym, hist_vals, hist_h, hist_k, hist_l, - hist_i); + hist_i, config_nopolar); if ( r == 0 ) n_used++; @@ -355,6 +361,7 @@ int main(int argc, char *argv[]) int hist_i; int space_for_hist = 0; char *histo_params = NULL; + int config_nopolar = 0; /* Long options */ const struct option longopts[] = { @@ -367,6 +374,8 @@ int main(int argc, char *argv[]) {"start-after", 1, NULL, 'f'}, {"sum", 0, &config_sum, 1}, {"scale", 0, &config_scale, 1}, + {"no-polarisation", 0, &config_nopolar, 1}, + {"no-polarization", 0, &config_nopolar, 1}, {"symmetry", 1, NULL, 'y'}, {"histogram", 1, NULL, 'g'}, {"hist-parameters", 1, NULL, 'z'}, @@ -503,7 +512,7 @@ int main(int argc, char *argv[]) hist_i = 0; merge_all(fh, model, NULL, config_startafter, config_stopafter, sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l, - &hist_i); + &hist_i, config_nopolar); if ( ferror(fh) ) { ERROR("Stream read error.\n"); return 1; @@ -526,7 +535,8 @@ int main(int argc, char *argv[]) merge_all(fh, model, reference, config_startafter, config_stopafter, sym, n_total_patterns, - hist_vals, hist_h, hist_k, hist_l, &hist_i); + hist_vals, hist_h, hist_k, hist_l, &hist_i, + config_nopolar); if ( ferror(fh) ) { ERROR("Stream read error.\n"); |