aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-03-16 11:11:11 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:20 +0100
commit4d654fcf91e742edd6940c4eda04113aa41ff8c8 (patch)
tree41a3aa9fe92f5a86cdb4ce3cd3ce16082a11126c /src
parente39c7d93bf878be9355e074475784119ab2e257c (diff)
Remove 'estimate_background'
We should be estimating proper SNRs with sigma(I) estimates
Diffstat (limited to 'src')
-rw-r--r--src/estimate_background.c231
1 files changed, 0 insertions, 231 deletions
diff --git a/src/estimate_background.c b/src/estimate_background.c
deleted file mode 100644
index 1f956377..00000000
--- a/src/estimate_background.c
+++ /dev/null
@@ -1,231 +0,0 @@
-/*
- * estimate-background.c
- *
- * Like 'process_hkl', but for signal/noise ratios
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdarg.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <unistd.h>
-#include <getopt.h>
-
-#include "utils.h"
-#include "statistics.h"
-#include "sfac.h"
-#include "reflections.h"
-#include "symmetry.h"
-#include "stream.h"
-
-
-/* Number of bins for plot of resolution shells */
-#define NBINS (10)
-
-
-
-static void show_help(const char *s)
-{
- printf("Syntax: %s [options]\n\n", s);
- printf(
-"Estimate peak SNR from FEL Bragg intensities.\n"
-"\n"
-" -h, --help Display this help message.\n"
-" -i, --input=<filename> Specify input filename (\"-\" for stdin).\n"
-);
-}
-
-
-
-int main(int argc, char *argv[])
-{
- int c;
- char *filename = NULL;
- FILE *fh;
- int rval;
- double total_vol, vol_per_shell;
- double rmin, rmax;
- double rmins[NBINS];
- double rmaxs[NBINS];
- double snrs[NBINS];
- unsigned int counts[NBINS];
- UnitCell *real_cell;
- double overall_snr;
- unsigned int overall_counts;
- int i;
-
- /* Long options */
- const struct option longopts[] = {
- {"help", 0, NULL, 'h'},
- {"input", 1, NULL, 'i'},
- {0, 0, NULL, 0}
- };
-
- /* Short options */
- while ((c = getopt_long(argc, argv, "hi:",
- longopts, NULL)) != -1) {
-
- switch (c) {
- case 'h' :
- show_help(argv[0]);
- return 0;
-
- case 'i' :
- filename = strdup(optarg);
- break;
-
- case 0 :
- break;
-
- default :
- return 1;
- }
-
- }
-
- if ( filename == NULL ) {
- ERROR("Please specify filename using the -i option\n");
- return 1;
- }
-
- /* Open the data stream */
- if ( strcmp(filename, "-") == 0 ) {
- fh = stdin;
- } else {
- fh = fopen(filename, "r");
- }
- free(filename);
- if ( fh == NULL ) {
- ERROR("Failed to open input file\n");
- return 1;
- }
-
- /* FIXME: Fixed resolution shells */
- rmin = 0.120e9;
- rmax = 1.172e9;
-
- for ( i=0; i<NBINS; i++ ) {
- counts[i] = 0;
- snrs[i] = 0.0;
- }
-
- total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
- vol_per_shell = total_vol / NBINS;
- rmins[0] = rmin;
- for ( i=1; i<NBINS; i++ ) {
-
- double r;
-
- r = vol_per_shell + pow(rmins[i-1], 3.0);
- r = pow(r, 1.0/3.0);
-
- /* Shells of constant volume */
- rmaxs[i-1] = r;
- rmins[i] = r;
-
- /* Shells of constant thickness */
- //rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS;
- //rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS;
-
- STATUS("Shell %i: %f to %f\n", i-1,
- rmins[i-1]/1e9, rmaxs[i-1]/1e9);
-
- }
- rmaxs[NBINS-1] = rmax;
- STATUS("Shell %i: %f to %f\n", NBINS-1,
- rmins[NBINS-1]/1e9, rmaxs[NBINS-1]/1e9);
-
- real_cell = load_cell_from_pdb("../../1JB0.pdb");
-
- do {
-
- char *rval2;
- UnitCell *cell;
- char *filename;
- char line[1024];
- int done = 0;
- double ph_ev;
-
- rval = find_chunk(fh, &cell, &filename, &ph_ev);
- if ( rval != 0 ) break;
-
- do {
-
- signed int h, k, l;
- float intensity, x, y, max, bg;
- int r, j;
- double d, snr;
- int bin;
-
- rval2 = fgets(line, 1024, fh);
-
- if ( strncmp(line, "Peak statistics:", 16) == 0 ) {
- done = 1;
- }
-
- r = sscanf(line, "%i %i %i %f (at %f,%f) max=%f bg=%f",
- &h, &k, &l, &intensity, &x, &y, &max, &bg);
-
- if ( r != 8 ) {
- continue;
- }
-
- d = resolution(real_cell, h, k, l) * 2.0;
- //STATUS("'%s'", line);
- //STATUS("-> %i %i %i %f\n", h, k, l, d/1e9);
-
- bin = -1;
- for ( j=0; j<NBINS; j++ ) {
- if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
- bin = j;
- break;
- }
- }
- if ( bin == -1 ) {
- ERROR("Warnung! %i %i %i %f\n", h, k, l, d/1e9);
- continue;
- }
-
- snr = max/bg;
- snrs[bin] += snr;
- counts[bin]++;
-
- } while ( !done );
-
- } while ( !rval );
-
- /* Print out results */
- overall_snr = 0.0;
- overall_counts = 0;
- for ( i=0; i<NBINS; i++ ) {
-
- double snr, cen;
-
- cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
- snr = snrs[i] / (double)counts[i];
- printf("%f %f (%f / %i)\n", cen*1.0e-9, snr, snrs[i], counts[i]);
- if ( i>0 ) {
- overall_snr += snrs[i];
- overall_counts += counts[i];
- }
-
- }
-
- printf("Overall: %f (%f / %i)\n", overall_snr/(double)overall_counts,
- overall_snr, overall_counts);
-
-
- fclose(fh);
-
- return 0;
-}