From 94b0050cc7735c3e1635cbc89c13c6b2c49c69c8 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 17 May 2019 17:03:48 +0200 Subject: Use Spectrum API for simulation --- src/diffraction-gpu.c | 26 ++-- src/diffraction-gpu.h | 4 +- src/diffraction.c | 329 +++----------------------------------------------- src/diffraction.h | 2 +- src/indexamajig.c | 63 +--------- src/partial_sim.c | 3 +- src/pattern_sim.c | 73 +++++++++-- src/process_image.h | 2 +- 8 files changed, 99 insertions(+), 403 deletions(-) (limited to 'src') diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c index 9cbfdf33..e2bce0cc 100644 --- a/src/diffraction-gpu.c +++ b/src/diffraction-gpu.c @@ -283,7 +283,7 @@ static int do_panels(struct gpu_context *gctx, struct image *image, int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, int na, int nb, int nc, UnitCell *ucell, - int no_fringes, int flat) + int no_fringes, int flat, int n_samples) { double ax, ay, az; double bx, by, bz; @@ -294,6 +294,8 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, int n_neg = 0; int n_nan = 0; int i; + double kmin, kmax, step; + double tot = 0.0; if ( gctx == NULL ) { ERROR("GPU setup failed.\n"); @@ -338,23 +340,25 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, } } - double tot = 0.0; - for ( i=0; insamples; i++ ) { + spectrum_get_range(image->spectrum, &kmin, &kmax); + step = (kmax-kmin)/n_samples; + for ( i=0; i<=n_samples; i++ ) { + + double k = kmin + i*step; + double prob; - printf("%.3f eV, weight = %.5f\n", - ph_lambda_to_eV(1.0/image->spectrum0[i].k), - image->spectrum0[i].weight); + /* Probability = p.d.f. times step width */ + prob = step * spectrum_get_density_at_k(image->spectrum, k); - err = do_panels(gctx, image, image->spectrum0[i].k, - image->spectrum0[i].weight, - &n_inf, &n_neg, &n_nan); + STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob); + err = do_panels(gctx, image, k, prob, &n_inf, &n_neg, &n_nan); if ( err ) return 1; - tot += image->spectrum0[i].weight; + tot += prob; } - printf("total weight = %f\n", tot); + STATUS("Total weight = %f\n", tot); if ( n_neg + n_inf + n_nan ) { ERROR("WARNING: The GPU calculation produced %i negative" diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h index bef12a8e..a4d16583 100644 --- a/src/diffraction-gpu.h +++ b/src/diffraction-gpu.h @@ -42,7 +42,7 @@ struct gpu_context; extern int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, int na, int nb, int nc, UnitCell *ucell, - int no_fringes, int flat); + int no_fringes, int flat, int n_samples); extern struct gpu_context *setup_gpu(int no_sfac, const double *intensities, const unsigned char *flags, @@ -53,7 +53,7 @@ extern void cleanup_gpu(struct gpu_context *gctx); static int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, int na, int nb, int nc, UnitCell *ucell, - int no_fringes, int flat) + int no_fringes, int flat, int n_samples) { /* Do nothing */ ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n"); diff --git a/src/diffraction.c b/src/diffraction.c index 4dc21c59..bcbc62e8 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -443,317 +443,11 @@ static void diffraction_at_k(struct image *image, const double *intensities, } -static void normalise_sampled_spectrum(struct sample *spec, int n, int nsamp) -{ - int i; - double total_w = 0.0; - double samp_w = 0.0; - - for ( i=0; iweight < sample2->weight ) { - return 1; - } - return -1; -} - - -static struct sample *get_gaussian_spectrum(double eV_cen, double eV_step, - double sigma, int spec_size, - double eV_start) -{ - struct sample *spectrum; - int i; - double eV; - - spectrum = malloc(spec_size * sizeof(struct sample)); - if ( spectrum == NULL ) return NULL; - - if (eV_start == 0) { /* eV starts at 3 sigma below the mean*/ - eV = eV_cen - (spec_size/2)*eV_step; - } else { - eV = eV_start; - } - - for ( i=0; i nsteps ) shiftLim = nsteps; - - noise = malloc(nsteps * sizeof(struct sample)); - if ( noise == NULL ) return 1; - - gaussianNoise = malloc(3 * nsteps * sizeof(double)); - if ( gaussianNoise == NULL ) { - free(noise); - return 1; - } - - /* Generate Gaussian noise of length of spectrum - * (replicate on both ends for circshift below) */ - for ( i=0; ibw * image->lambda / 2.0; /* m */ - double mink = 1.0/(image->lambda + halfwidth); - double maxk = 1.0/(image->lambda - halfwidth); - - spectrum = malloc(image->nsamples * sizeof(struct sample)); - if ( spectrum == NULL ) return NULL; - - k = mink; - k_step = (maxk-mink)/(image->nsamples-1); - for ( i=0; insamples; i++ ) { - spectrum[i].k = k; - spectrum[i].weight = 1.0/(double)image->nsamples; - k += k_step; - } - - image->spectrum_size = image->nsamples; - /* No need to call normalise_sampled_spectrum() in this case */ - - return spectrum; -} - - -struct sample *generate_spectrum_fromfile(struct image *image, - char *spectrum_fn) -{ - struct sample *spectrum; - int i; - double k, w; - double w_sum = 0; - - spectrum = malloc(image->nsamples * sizeof(struct sample)); - if ( spectrum == NULL ) return NULL; - - FILE *f; - f = fopen(spectrum_fn, "r"); - - int nsamples = 0; - for ( i=0; insamples; i++ ) { - if (fscanf(f, "%lf %lf", &k, &w) != EOF) { - spectrum[i].k = ph_eV_to_k(k); - spectrum[i].weight = w; - w_sum += w; - nsamples += 1; - } else break; - } - - for ( i=0; ispectrum_size = nsamples; - - return spectrum; -} - - -struct sample *generate_SASE(struct image *image, gsl_rng *rng) -{ - struct sample *spectrum; - const int spec_size = 1024; - double eV_cen; /* Central photon energy for this spectrum */ - const double jitter_sigma_eV = 8.0; - - /* Central wavelength jitters with Gaussian distribution */ - eV_cen = gaussian_noise(rng, ph_lambda_to_eV(image->lambda), - jitter_sigma_eV); - - /* Convert FWHM to standard deviation. Note that bandwidth is taken - * here to be "delta E over E" (E = photon energy), not the bandwidth in - * terms of wavelength (as it is everywhere else), but the difference - * should be very small */ - double sigma = (image->bw*eV_cen) / (2.0*sqrt(2.0*log(2.0))); - - /* The spectrum will be calculated to a resolution which spreads six - * sigmas of the original (no SASE noise) Gaussian pulse over spec_size - * points */ - double eV_step = 6.0*sigma/(spec_size-1); - - spectrum = get_gaussian_spectrum(eV_cen, eV_step, sigma, spec_size,0); - - /* Add SASE-type noise to Gaussian spectrum */ - add_sase_noise(spectrum, spec_size, rng); - - /* Sort samples in spectrum by weight. Diffraction calculation will - * take the requested number, starting from the brightest */ - qsort(spectrum, spec_size, sizeof(struct sample), compare_samples); - - normalise_sampled_spectrum(spectrum, spec_size, image->nsamples); - - image->spectrum_size = spec_size; - return spectrum; -} - - -struct sample *generate_twocolour(struct image *image) -{ - struct sample *spectrum; - struct sample *spectrum1; - struct sample *spectrum2; - int i; - double eV_cen1; /* Central photon energy for first colour */ - double eV_cen2; /* Central photon energy for second colour */ - double eV_cen; /* Central photon energy for this spectrum */ - const int spec_size = 1024; - - eV_cen = ph_lambda_to_eV(image->lambda); - - /* Convert FWHM to standard deviation. Note that bandwidth is taken - * here to be "delta E over E" (E = photon energy), not the bandwidth in - * terms of wavelength (as it is everywhere else), but the difference - * should be very small */ - double halfwidth = eV_cen*image->bw/2.0; /* eV */ - - eV_cen1 = eV_cen - halfwidth; - eV_cen2 = eV_cen + halfwidth; - - /* Hard-code sigma to be 1/5 of bandwidth */ - double sigma = eV_cen*image->bw/5.0; /* eV */ - - /* The spectrum will be calculated to a resolution which spreads six - * sigmas of the original (no SASE noise) Gaussian pulse over spec_size - * points */ - double eV_start = eV_cen1 - 3*sigma; - double eV_end = eV_cen2 + 3*sigma; - double eV_step = (eV_end - eV_start)/(spec_size-1); - - spectrum1 = get_gaussian_spectrum(eV_cen1, eV_step, sigma, spec_size, - eV_start); - spectrum2 = get_gaussian_spectrum(eV_cen2, eV_step, sigma, spec_size, - eV_start); - - spectrum = malloc(spec_size * sizeof(struct sample)); - if ( spectrum == NULL ) return NULL; - - for ( i=0; insamples); - - image->spectrum_size = spec_size; - return spectrum; -} - - void get_diffraction(struct image *image, int na, int nb, int nc, const double *intensities, const double *phases, const unsigned char *flags, UnitCell *cell, GradientMethod m, const SymOpList *sym, - int no_fringes, int flat) + int no_fringes, int flat, int n_samples) { double ax, ay, az; double bx, by, bz; @@ -762,6 +456,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc, double *lut_b; double *lut_c; int i; + double kmin, kmax, step; cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); @@ -769,17 +464,23 @@ void get_diffraction(struct image *image, int na, int nb, int nc, lut_b = get_sinc_lut(nb, no_fringes, flat); lut_c = get_sinc_lut(nc, no_fringes, flat); - for ( i=0; insamples; i++ ) { + spectrum_get_range(image->spectrum, &kmin, &kmax); + step = (kmax-kmin)/(n_samples+1); + STATUS("%i samples\n", n_samples); + for ( i=1; i<=n_samples; i++ ) { - printf("%.1f eV, weight = %.5f\n", - ph_lambda_to_eV(1.0/image->spectrum0[i].k), - image->spectrum0[i].weight); + double k = kmin + i*step; + double prob; + + /* Probability = p.d.f. times step width */ + prob = step * spectrum_get_density_at_k(image->spectrum, k); + + STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob); diffraction_at_k(image, intensities, phases, - flags, cell, m, sym, image->spectrum0[i].k, + flags, cell, m, sym, k, ax, ay, az, bx, by, bz, cx, cy, cz, - lut_a, lut_b, lut_c, image->spectrum0[i].weight); - + lut_a, lut_b, lut_c, prob); } diff --git a/src/diffraction.h b/src/diffraction.h index 0f95f83c..16b0386f 100644 --- a/src/diffraction.h +++ b/src/diffraction.h @@ -52,7 +52,7 @@ extern void get_diffraction(struct image *image, int na, int nb, int nc, const double *intensities, const double *phases, const unsigned char *flags, UnitCell *cell, GradientMethod m, const SymOpList *sym, - int no_fringes, int flat); + int no_fringes, int flat, int n_samples); extern struct sample *generate_tophat(struct image *image); diff --git a/src/indexamajig.c b/src/indexamajig.c index 681d5d31..e4568966 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -250,66 +250,6 @@ static void add_geom_beam_stuff_to_field_list(struct imagefile_field_list *copym } -static struct spectrum *read_spectrum_fromfile(char *fn) -{ - FILE *f; - struct spectrum *s; - int i; - double k, w; - double w_sum = 0; - - f = fopen(fn, "r"); - if ( f == NULL ) { - ERROR("Couldn't open '%s'\n", fn); - return NULL; - } - - s = malloc(sizeof(struct spectrum)); - if ( s == NULL ) return NULL; - - if ( fscanf(f, "%d", &s->n) == EOF ) { - return NULL; - } - - if ( s->n <= 0 ) { - return NULL; - } - - s->ks = malloc(s->n * sizeof(double)); - if ( s->ks == NULL ) { - ERROR("Failed to allocate spectrum!\n"); - return NULL; - } - - s->weights = malloc(s->n * sizeof(double)); - if ( s->weights == NULL ) { - ERROR("Failed to allocate spectrum!\n"); - return NULL; - } - - for ( i=0; in; i++ ) { - if ( fscanf(f, "%lf %lf", &k, &w) != EOF ) { - s->ks[i] = ph_eV_to_k(k); - s->weights[i] = w; - w_sum += w; - } else { - break; - } - } - - if ( i < s->n - 1 ) { - ERROR("Failed to read %d lines from %s\n", s->n, fn); - return NULL; - } - - for ( i=0; in; i++ ) { - s->weights[i] /= w_sum; - } - - return s; -} - - int main(int argc, char *argv[]) { int c; @@ -1164,13 +1104,12 @@ int main(int argc, char *argv[]) /* Load spectrum from file if given */ if ( spectrum_fn != NULL ) { - iargs.spectrum = read_spectrum_fromfile(spectrum_fn); + iargs.spectrum = spectrum_load(spectrum_fn); if ( iargs.spectrum == NULL ) { ERROR("Couldn't read spectrum (from %s)\n", spectrum_fn); return 1; } free(spectrum_fn); - STATUS("Read %d lines from %s\n", iargs.spectrum->n, spectrum_fn); } else { iargs.spectrum = NULL; } diff --git a/src/partial_sim.c b/src/partial_sim.c index 92f9703c..0ea45ad4 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -841,8 +841,7 @@ int main(int argc, char *argv[]) image.crystals = NULL; image.n_crystals = 0; image.indexed_by = INDEXING_SIMULATION; - image.spectrum_size = 0; - image.spectrum0 = NULL; + image.spectrum = NULL; image.serial = 0; image.event = NULL; diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 5ee53514..d13baff8 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -60,6 +60,7 @@ enum spectrum_type { SPECTRUM_TOPHAT, /**< A top hat distribution */ + SPECTRUM_GAUSSIAN, /**< A Gaussian distribution */ SPECTRUM_SASE, /**< A simulated SASE spectrum */ SPECTRUM_TWOCOLOUR, /**< A spectrum containing two peaks */ SPECTRUM_FROMFILE /**< An arbitrary spectrum read from a file */ @@ -102,7 +103,12 @@ static void show_help(const char *s) " --template= Take orientations from stream .\n" " --no-fringes Exclude the side maxima of Bragg peaks.\n" " --flat Make Bragg peaks flat.\n" -" --beam-bandwidth Beam bandwidth as a fraction. Default 1%%.\n" +" --beam-bandwidth Beam bandwidth (standard deviation of wavelength as\n" +" a fraction of wavelenth). Default 0.001 (1%%)\n" +" --sase-spike-width SASE spike width (standard deviation in m^-1)\n" +" Default 2e6 m^-1\n" +" --twocol-separation Separation between peaks in two-colour mode in m^-1\n" +" Default 8e6 m^-1\n" " --photon-energy Photon energy in eV. Default 9000.\n" " --nphotons Number of photons per X-ray pulse. Default 1e12.\n" " --beam-radius Radius of beam in metres (default 1e-6).\n" @@ -401,6 +407,8 @@ int main(int argc, char *argv[]) double beam_radius = 1e-6; /* metres */ double bandwidth = 0.01; double photon_energy = 9000.0; + double sase_spike_width = 2e6; /* m^-1 */ + double twocol_sep = 8e6; /* m^-1 */ struct beam_params beam; int i; @@ -438,6 +446,8 @@ int main(int argc, char *argv[]) {"nphotons", 1, NULL, 10}, {"beam-radius", 1, NULL, 11}, {"spectrum-file", 1, NULL, 12}, + {"sase-spike-width", 1, NULL, 13}, + {"twocol-separation", 1, NULL, 14}, {0, 0, NULL, 0} }; @@ -603,6 +613,30 @@ int main(int argc, char *argv[]) spectrum_fn = strdup(optarg); break; + case 13 : + sase_spike_width = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid SASE spike width.\n"); + return 1; + } + if ( beam_radius < 0.0 ) { + ERROR("SASE spike width must be positive.\n"); + return 1; + } + break; + + case 14 : + twocol_sep = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid two colour separation.\n"); + return 1; + } + if ( beam_radius < 0.0 ) { + ERROR("Two-colour separation must be positive.\n"); + return 1; + } + break; + case 0 : break; @@ -700,6 +734,8 @@ int main(int argc, char *argv[]) spectrum_type = SPECTRUM_TOPHAT; } else if ( strcasecmp(spectrum_str, "tophat") == 0) { spectrum_type = SPECTRUM_TOPHAT; + } else if ( strcasecmp(spectrum_str, "gaussian") == 0) { + spectrum_type = SPECTRUM_GAUSSIAN; } else if ( strcasecmp(spectrum_str, "sase") == 0) { spectrum_type = SPECTRUM_SASE; } else if ( strcasecmp(spectrum_str, "twocolour") == 0 || @@ -789,7 +825,6 @@ int main(int argc, char *argv[]) image.lambda = ph_en_to_lambda(eV_to_J(photon_energy)); image.bw = bandwidth; - image.nsamples = nsamples; /* Initialise stuff */ image.filename = NULL; @@ -826,7 +861,7 @@ int main(int argc, char *argv[]) powder.det = image.det; powder.beam = NULL; powder.lambda = 0.0; - powder.spectrum0 = NULL; + powder.spectrum = NULL; powder.dp = malloc(image.det->n_panels*sizeof(float *)); if ( powder.dp == NULL ) { ERROR("Failed to allocate powder data\n"); @@ -860,6 +895,11 @@ int main(int argc, char *argv[]) "width %.5f %%\n", image.bw*100.0); break; + case SPECTRUM_GAUSSIAN: + STATUS(" X-ray spectrum: Gaussian, " + "bandwidth %.5f %%\n", image.bw*100.0); + break; + case SPECTRUM_SASE: STATUS(" X-ray spectrum: SASE, " "bandwidth %.5f %%\n", image.bw*100.0); @@ -985,20 +1025,31 @@ int main(int argc, char *argv[]) switch ( spectrum_type ) { case SPECTRUM_TOPHAT : - image.spectrum0 = generate_tophat(&image); + image.spectrum = spectrum_generate_tophat(image.lambda, + image.bw); break; + case SPECTRUM_GAUSSIAN : + image.spectrum = spectrum_generate_gaussian(image.lambda, + image.bw); + break; + + case SPECTRUM_SASE : - image.spectrum0 = generate_SASE(&image, rng); + image.spectrum = spectrum_generate_sase(image.lambda, + image.bw, + sase_spike_width, + rng); break; case SPECTRUM_TWOCOLOUR : - image.spectrum0 = generate_twocolour(&image); + image.spectrum = spectrum_generate_twocolour(image.lambda, + image.bw, + twocol_sep); break; case SPECTRUM_FROMFILE : - image.spectrum0 = generate_spectrum_fromfile(&image, - spectrum_fn); + image.spectrum = spectrum_load(spectrum_fn); break; } @@ -1016,11 +1067,13 @@ int main(int argc, char *argv[]) gpu_dev); } err = get_diffraction_gpu(gctx, &image, na, nb, nc, - cell, no_fringes, flat); + cell, no_fringes, flat, + nsamples); } else { get_diffraction(&image, na, nb, nc, intensities, phases, - flags, cell, grad, sym, no_fringes, flat); + flags, cell, grad, sym, no_fringes, flat, + nsamples); } if ( err ) { ERROR("Diffraction calculation failed.\n"); diff --git a/src/process_image.h b/src/process_image.h index 2a43d11d..5939ae4d 100644 --- a/src/process_image.h +++ b/src/process_image.h @@ -112,7 +112,7 @@ struct index_args struct taketwo_options taketwo_opts; struct xgandalf_options xgandalf_opts; struct felix_options felix_opts; - struct spectrum *spectrum; + Spectrum *spectrum; signed int wait_for_file; /* -1 means wait forever */ }; -- cgit v1.2.3