diff options
62 files changed, 2308 insertions, 593 deletions
@@ -109,3 +109,4 @@ Contributors * Yaroslav Gevorkov <yaroslav.gevorkov@desy.de> peakfinder9 + xgandalf diff --git a/CMakeLists.txt b/CMakeLists.txt index 8f5f82f0..a3f15912 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,12 +12,51 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/CMake") find_package(HDF5 REQUIRED COMPONENTS C) find_package(Threads REQUIRED) find_package(GSL REQUIRED) -find_package(GTK2 2.0.0 COMPONENTS gtk) find_package(TIFF) find_package(OpenCL) find_package(PkgConfig) + +pkg_check_modules (GLIB glib-2.0) +include_directories (${GLIB_INCLUDE_DIRS}) +link_directories (${GLIB_LIBRARY_DIRS}) + +set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wall") + +pkg_search_module(GTK gtk+-3.0) +if (NOT GTK_FOUND) + + # No GTK3. Try GTK2? + pkg_search_module(GTK gtk+-2.0) + if (GTK_FOUND) + message(STATUS "Using GTK 2") + pkg_search_module(GDK gdk-2.0) + endif () + +else () + + message(STATUS "Using GTK 3") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wno-deprecated-declarations") + pkg_search_module(GDK gdk-3.0) + +endif () + +# If no version of GTK was found, try for GDK +if (NOT GTK_FOUND) + message(STATUS "GTK not found. hdfsee and cell_explorer will not be compiled") + message(STATUS "Looking separately for GDK") + pkg_search_module(GDK gdk-3.0) + if (NOT GDK_FOUND) + pkg_search_module(GDK gdk-2.0) + if (GDK_FOUND) + message(STATUS "Using GDK 2") + endif() + else () + message(STATUS "Using GDK 3") + endif() +endif() +link_directories(${GTK_LIBRARY_DIRS}) + pkg_search_module(CAIRO cairo) -pkg_search_module(GDK gdk-2.0) pkg_search_module(GDKPIXBUF gdk-pixbuf-2.0) set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fdiagnostics-color=always") @@ -53,7 +92,7 @@ include(CheckLibraryExists) set(HAVE_CAIRO ${CAIRO_FOUND}) set(HAVE_TIFF ${TIFF_FOUND}) -set(HAVE_GTK ${GTK2_FOUND}) +set(HAVE_GTK ${GTK_FOUND}) set(HAVE_OPENCL ${OpenCL_FOUND}) set(HAVE_GDKPIXBUF ${GDKPIXBUF_FOUND}) set(HAVE_GDK ${GDK_FOUND}) @@ -64,9 +103,16 @@ check_include_file(fcntl.h HAVE_FCNTL_H) check_include_file(stdlib.h HAVE_STDLIB_H) check_include_file(unistd.h HAVE_UNISTD_H) -# If we supposedly have OpenCL, check if the header is actually "CL/cl.h" +# Workarounds for weird OpenCL quirks.. if (OpenCL_FOUND) + + # If include path ends with "OpenCL.framework", add the missing "/Headers" + string(REGEX REPLACE "OpenCL.framework$" "OpenCL.framework/Headers" + OpenCL_INCLUDE_DIRS ${OpenCL_INCLUDE_DIRS}) + + # Check if the header is actually "CL/cl.h" check_include_file(CL/cl.h HAVE_CL_CL_H "-I${OpenCL_INCLUDE_DIRS}") + endif (OpenCL_FOUND) # Check for nice clock function @@ -136,13 +182,13 @@ list(APPEND CRYSTFEL_EXECUTABLES pattern_sim) # ---------------------------------------------------------------------- # hdfsee -if (GTK2_FOUND) +if (GTK_FOUND) set(HDFSEE_SOURCES src/hdfsee.c src/dw-hdfsee.c src/hdfsee-render.c) add_executable(hdfsee ${HDFSEE_SOURCES}) - target_include_directories(hdfsee PRIVATE ${COMMON_INCLUDES} ${GTK2_INCLUDE_DIRS}) - target_link_libraries(hdfsee ${COMMON_LIBRARIES} ${GTK2_LIBRARIES}) + target_include_directories(hdfsee PRIVATE ${COMMON_INCLUDES} ${GTK_INCLUDE_DIRS}) + target_link_libraries(hdfsee ${COMMON_LIBRARIES} ${GTK_LIBRARIES}) if (TIFF_FOUND) target_include_directories(hdfsee PRIVATE ${TIFF_INCLUDE_DIRS}) @@ -151,22 +197,22 @@ if (GTK2_FOUND) list(APPEND CRYSTFEL_EXECUTABLES hdfsee) -endif (GTK2_FOUND) +endif (GTK_FOUND) # ---------------------------------------------------------------------- # cell_explorer -if (GTK2_FOUND) +if (GTK_FOUND) set(CELL_EXPLORER_SOURCES src/cell_explorer.c src/multihistogram.c) add_executable(cell_explorer ${CELL_EXPLORER_SOURCES}) - target_include_directories(cell_explorer PRIVATE ${COMMON_INCLUDES} ${GTK2_INCLUDE_DIRS}) - target_link_libraries(cell_explorer ${COMMON_LIBRARIES} ${GTK2_LIBRARIES}) + target_include_directories(cell_explorer PRIVATE ${COMMON_INCLUDES} ${GTK_INCLUDE_DIRS}) + target_link_libraries(cell_explorer ${COMMON_LIBRARIES} ${GTK_LIBRARIES}) list(APPEND CRYSTFEL_EXECUTABLES cell_explorer) -endif (GTK2_FOUND) +endif (GTK_FOUND) # ---------------------------------------------------------------------- # render_hkl @@ -1,4 +1,25 @@ -CrystFEL version 0.7.0, 17th May 2018 +Changes in this development version up to commit 019220ed +--------------------------------------------------------- + +- Switch to CMake build system +- indexamajig: Warn if a frame takes more than 10 minutes to process +- Add support for peakfinder9 +- Add support for XGANDALF indexer +- check_{peak_detection,near_bragg}: Warn if stream is older than image +- Some fixes for peakfinder8 +- indexamajig: Clean up temporary folders on clean exit +- Use GTK3 if available +- partialator: Add deltaCChalf rejection test +- turbo-index-lsf: Match turbo-index-slurm as closely as possible +- sum-peaks: Accept "-" for stdin +- Handle 16bpp and gzipped CBF files +- Multiple fixes for XDS indexer interface +- Fix for indexing centered cells with Mosflm without prior cell parameters +- Handle numbers (e.g. "k-2h") in symmetry operator specifications (e.g. + ambigator) + + +CrystFEL version 0.7.0, 17th May 2018 -------------------------------------- - indexamajig: Auto-detection of indexing methods diff --git a/doc/man/crystfel_geometry.5 b/doc/man/crystfel_geometry.5 index cf64d76e..db444214 100644 --- a/doc/man/crystfel_geometry.5 +++ b/doc/man/crystfel_geometry.5 @@ -28,10 +28,6 @@ through. The coordinates are specified inclusively, meaning that a minimum of 0 and a maximum of 9 results in a width of ten pixels. Counting begins from zero. All pixels in the image must be assigned to a panel - gaps are not permitted. .PP -In the current version, panels are assumed to be perpendicular to the incident -beam and to have their edges parallel. Within these limitations, any geometry -can be constructed. - The job of the geometry file is to establish a relationship an array of pixel values in the HDF5 file, defined in terms only of the "fast scan" and "slow scan" directions, and the laboratory coordinate system defined as follows: @@ -166,7 +162,7 @@ The camera length offset (in metres) for this panel. This number will be added .IP \fBfs\fR .IP \fBss\fR .PD -For this panel, the fast and slow scan directions correspond to the given directions in the lab coordinate system described above, measured in pixels. Example: "panel0/fs = 0.5x+0.5y-0.0001z". Including a component in the z direction means that the panel normal is not parallel to the X-ray beam. +For this panel, the fast and slow scan directions correspond to the given directions in the lab coordinate system described above, measured in pixels. Example: "panel0/fs = 0.5x+0.5y-0.0001z". Including a component in the z direction means that the panel is not perpendicular to the X-ray beam. .PD 0 .IP \fBcorner_x\fR diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1 index 37e80801..ab581120 100644 --- a/doc/man/indexamajig.1 +++ b/doc/man/indexamajig.1 @@ -86,6 +86,10 @@ Invoke XDS, and use its REFIDX procedure to attempt to index the pattern. .PD Use the TakeTwo algorithm. See Ginn et al., Acta Cryst. (2016). D72, 956-965. +.IP \fBxgandalf\fR +.PD +Invoke XGANDALF - eXtended GrAdient Descent Algorithm for Lattice Finding. Xgandalf must be installed in order to be able to use it. + .PP You can add one or more of the following to the above indexing methods, to control what information should be provided to them. Note that indexamajig performs a series of checks on the indexing results, including checking that the result is consistent with the target unit cell parameters. To get completely "raw" indexing, you need to disable these checks (see below) \fBand\fR not provide prior information. @@ -397,6 +401,26 @@ The defaults are: \fB--taketwo-member-threshold=20\fR, \fB--taketwo-len-tolernac .PD 0 These set low-level parameters for the Felix indexing algorithm. +.PD 0 +.IP \fB--xgandalf-sampling-pitch=\fIn\fR +.IP \fB--xgandalf-grad-desc-iterations=\fIn\fR +.IP \fB--xgandalf-tolerance=\fIn\fR +.IP \fB--xgandalf-no-deviation-from-provided-cell\fR +.IP \fB--xgandalf-max-lattice-vector-length=\fIn\fR +.IP \fB--xgandalf-min-lattice-vector-length=\fIn\fR +.PD +These set low-level parameters for the XGANDALF indexing algorithm. +.IP +\fB--xgandalf-sampling-pitch\fR selects how dense the reciprocal space is sampled. [0-4]: extremelyLoose to extremelyDense. [5-7]: standardWithSeondaryMillerIndices to extremelyDenseWithSeondaryMillerIndices. Default is 6 (denseWithSeondaryMillerIndices). +.IP +\fB--xgandalf-grad-desc-iterations\fR selects how many gradient descent iterations are performed. [0-5]: veryFew to extremelyMany. Default is 4 (manyMany). +.IP +\fB--xgandalf-tolerance\fR relative tolerance of the lattice vectors. Default is 0.02. +.IP +\fB--xgandalf-no-deviation-from-provided-cell\fR if a prior unit cell was provided, and this flag is set, the found unit cell will have exactly the same size as the provided one. +.IP +\fB--xgandalf-min-lattice-vector-length\fR and \fB--xgandalf-min-lattice-vector-length\fR minimum and maximum possible lattice vector lengths (unit is A). Used for fitting without prior lattice as starting point for gradient descent, so the final minimum lattice vector length can be smaller/highier as min/max. Note: This is valid for the uncentered cell, i.e. the P-cell! Default is 30A and 250A respectively. + .SH INTEGRATION OPTIONS .PD 0 .IP \fB--integration=\fR\fImethod\fR diff --git a/doc/man/partialator.1 b/doc/man/partialator.1 index aa25d2a7..7a4b3310 100644 --- a/doc/man/partialator.1 +++ b/doc/man/partialator.1 @@ -67,11 +67,21 @@ Run \fIn\fR cycles of scaling and post refinement. Disable the scaling part of the refinement calculation. .PD 0 +.IP \fB--no-Bscale\fR +.PD +Disable the Debye-Waller part of the scaling calculation. + +.PD 0 .IP \fB--no-pr\fR .PD Disable the orientation/physics model part of the refinement calculation. .PD 0 +.IP \fB--no-deltacchalf\fR +.PD +Disable rejection based on deltaCChalf. + +.PD 0 .IP "\fB-m\fR \fImodel\fR" .IP \fB--model=\fR\fImodel\fR .PD @@ -93,6 +103,11 @@ Disable the polarisation correction. Include reflections only if their peak values were less than \fIn\fR. That means, \fIn\fR is the saturation value of the detector. The default is infinity, i.e. no cutoff. .PD 0 +.IP \fB--min-res=\fR\fIn\fR +.PD +Merge crystals only if they diffract to beyond \fIn\fR Angstroms resolution. The default is infinity, i.e. all crystals are included. The resolution is taken from the diffraction_resolution_limit line in the stream. + +.PD 0 .IP \fB--min-measurements=\fR\fIn\fR .PD Include a reflection in the output only if it appears at least least \fIn\fR times. The default is \fB--min-measurements=2\fR. diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt index f6f25ec9..3d7bec38 100644 --- a/libcrystfel/CMakeLists.txt +++ b/libcrystfel/CMakeLists.txt @@ -6,6 +6,7 @@ find_package(XGANDALF) find_package(PINKINDEXER) find_package(NBP) find_package(FDIP) +find_package(ZLIB REQUIRED) pkg_search_module(FFTW fftw3) set(HAVE_CURSES ${CURSES_FOUND}) @@ -68,6 +69,7 @@ set(LIBCRYSTFEL_SOURCES src/felix.c src/peakfinder8.c src/taketwo.c + src/xgandalf.c ) if (HAVE_FFTW) @@ -106,6 +108,7 @@ set(LIBCRYSTFEL_HEADERS src/felix.h src/peakfinder8.h src/taketwo.h + src/xgandalf.h ) add_library(${PROJECT_NAME} SHARED @@ -123,8 +126,9 @@ target_include_directories(${PROJECT_NAME} INTERFACE ${PROJECT_SOURCE_DIR}/src) include_directories(${CMAKE_CURRENT_BINARY_DIR}) add_definitions(-DHAVE_CONFIG_H) -target_include_directories(${PROJECT_NAME} PRIVATE ${HDF5_INCLUDE_DIRS}) -target_link_libraries(${PROJECT_NAME} PRIVATE util ${HDF5_C_LIBRARIES} Threads::Threads GSL::gsl m) +target_include_directories(${PROJECT_NAME} PRIVATE ${HDF5_INCLUDE_DIRS} ${ZLIB_INCLUDE_DIRS}) +target_link_libraries(${PROJECT_NAME} PRIVATE util ${HDF5_C_LIBRARIES} ${ZLIB_LIBRARIES} + Threads::Threads GSL::gsl m) if (XGANDALF_FOUND) target_include_directories(${PROJECT_NAME} PRIVATE ${XGANDALF_INCLUDES}) @@ -162,9 +166,11 @@ if (CURSES_FOUND) endif (CURSES_FOUND) target_compile_options(${PROJECT_NAME} PRIVATE -Wall) +set_target_properties(${PROJECT_NAME} PROPERTIES PUBLIC_HEADER "${LIBCRYSTFEL_HEADERS}") -install (TARGETS libcrystfel - LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}) +install(TARGETS libcrystfel + PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/crystfel + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}) # crystfel.pc configure_file(crystfel.pc.in crystfel.pc) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 5465052c..7b1984bb 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -208,7 +208,10 @@ int right_handed(UnitCell *cell) rh_direct = aCb_dot_c > 0.0; - assert(rh_reciprocal == rh_direct); + if ( rh_reciprocal != rh_direct ) { + ERROR("Whoops, reciprocal and real space handedness are " + "not the same!\n"); + } return rh_direct; } diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h index f4c7477b..cc5754ca 100644 --- a/libcrystfel/src/crystal.h +++ b/libcrystfel/src/crystal.h @@ -36,7 +36,6 @@ #include "cell.h" -#include "reflist.h" /** @@ -47,6 +46,8 @@ **/ typedef struct _crystal Crystal; +#include "reflist.h" + #ifdef __cplusplus extern "C" { #endif diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index 8d33731a..d4d8b768 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -1429,6 +1429,8 @@ struct detector *get_detector_geometry_2(const char *filename, if ( det->panels[i].dim_structure->dims[di] == HYSL_UNDEFINED ) { dim_dim_reject = 1; + ERROR("Dimension %i for panel %s is undefined.\n", + di, det->panels[i].name); } if ( det->panels[i].dim_structure->dims[di] == HYSL_PLACEHOLDER ) { @@ -1446,18 +1448,23 @@ struct detector *get_detector_geometry_2(const char *filename, } if ( found_ss != 1 ) { - ERROR("Only one slow scan dim coordinate is allowed\n"); + ERROR("Exactly one slow scan dim coordinate is needed " + "(found %i for panel %s)\n", found_ss, + det->panels[i].name); dim_dim_reject = 1; } if ( found_fs != 1 ) { - ERROR("Only one fast scan dim coordinate is allowed\n"); + ERROR("Exactly one fast scan dim coordinate is needed " + "(found %i for panel %s)\n", found_fs, + det->panels[i].name); dim_dim_reject = 1; } if ( panel_dim_dim > 1 ) { ERROR("Maximum one placeholder dim coordinate is " - "allowed\n"); + "allowed (found %i for panel %s)\n", + panel_dim_dim, det->panels[i].name); dim_dim_reject = 1; } diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index 7e3d6b11..f9c2aa97 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -1047,9 +1047,9 @@ int hdf5_write_image(const char *filename, const struct image *image, write_photon_energy(fh, ph_lambda_to_eV(image->lambda), ph_en_loc); - if ( image->spectrum != NULL && image->spectrum_size > 0 ) { + if ( image->spectrum0 != NULL && image->spectrum_size > 0 ) { - write_spectrum(fh, image->spectrum, image->spectrum_size, + write_spectrum(fh, image->spectrum0, image->spectrum_size, image->nsamples); } diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 0dda66ed..26d48007 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -34,6 +34,7 @@ #include <math.h> #include <stdio.h> #include <hdf5.h> +#include <zlib.h> #ifdef HAVE_CBFLIB #ifdef HAVE_CBFLIB_CBF_H @@ -439,7 +440,7 @@ static char *cbf_strerr(int e) } -static int unpack_panels(struct image *image, signed int *data, int data_width, +static int unpack_panels(struct image *image, float *data, int data_width, int data_height) { int pi; @@ -587,7 +588,39 @@ static void cbf_fill_in_clen(struct detector *det, struct imagefile *f) } -static int read_cbf(struct imagefile *f, struct image *image) +static float *convert_sint32_float(int32_t *data, int w, int h) +{ + float *df; + long int i; + + df = malloc(sizeof(float)*w*h); + if ( df == NULL ) return NULL; + + for ( i=0; i<w*h; i++ ) { + df[i] = data[i]; + } + + return df; +} + + +static float *convert_sint16_float(int16_t *data, int w, int h) +{ + float *df; + long int i; + + df = malloc(sizeof(float)*w*h); + if ( df == NULL ) return NULL; + + for ( i=0; i<w*h; i++ ) { + df[i] = data[i]; + } + + return df; +} + + +static float *read_cbf_data(struct imagefile *f, int *w, int *h, cbf_handle *pcbfh) { cbf_handle cbfh; FILE *fh; @@ -596,29 +629,70 @@ static int read_cbf(struct imagefile *f, struct image *image) int binary_id, minelement, maxelement, elsigned, elunsigned; size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding; const char *byteorder; - signed int *data; + void *data; + float *dataf; + void *buf = NULL; - if ( image->det == NULL ) { - ERROR("read_cbf() needs a geometry\n"); - return 1; + if ( f->type == IMAGEFILE_CBF ) { + + fh = fopen(f->filename, "rb"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", f->filename); + return NULL; + } + + } else if ( f->type == IMAGEFILE_CBFGZ ) { + + gzFile gzfh; + size_t len, len_read; + const size_t bufinc = 8*1024*1024; /* Allocate buffer in 8Mb chunks */ + size_t bufsz = bufinc; + + gzfh = gzopen(f->filename, "rb"); + if ( gzfh == NULL ) return NULL; + + /* Set larger buffer size for hopefully faster uncompression */ + gzbuffer(gzfh, 128*1024); + + buf = malloc(bufsz); + if ( buf == NULL ) return NULL; + + len = 0; + do { + + len_read = gzread(gzfh, buf+len, bufinc); + if ( len_read == -1 ) return NULL; + len += len_read; + + if ( len_read == bufinc ) { + bufsz += bufinc; + buf = realloc(buf, bufsz); + if ( buf == NULL ) return NULL; + } + + } while ( len_read == bufinc ); + + fh = fmemopen(buf, len, "rb"); + if ( fh == NULL ) return NULL; + + gzclose_r(gzfh); + + } else { + /* Don't know how we ended up here */ + return NULL; } if ( cbf_make_handle(&cbfh) ) { ERROR("Failed to allocate CBF handle\n"); - return 1; + return NULL; } - fh = fopen(f->filename, "rb"); - if ( fh == NULL ) { - ERROR("Failed to open '%s'\n", f->filename); - return 1; - } - /* CBFlib calls fclose(fh) when it's ready */ + /* CBFlib will call fclose(fh) when it's ready */ if ( cbf_read_widefile(cbfh, fh, 0) ) { ERROR("Failed to read CBF file '%s'\n", f->filename); cbf_free_handle(cbfh); - return 1; + return NULL; } /* Select row 0 in data column inside array_data */ @@ -639,40 +713,46 @@ static int read_cbf(struct imagefile *f, struct image *image) ERROR("Failed to read CBF array parameters: %s\n", err); free(err); cbf_free_handle(cbfh); - return 1; + return NULL; } if ( dimslow != 0 ) { ERROR("CBF data array is 3D - don't know what to do with it\n"); cbf_free_handle(cbfh); - return 1; + return NULL; } if ( dimfast*dimmid*elsize > 10e9 ) { ERROR("CBF data is far too big (%i x %i x %i bytes).\n", (int)dimfast, (int)dimmid, (int)elsize); cbf_free_handle(cbfh); - return 1; + return NULL; } - if ( elsize != 4 ) { + if ( (elsize != 4) && (elsize != 2) ) { STATUS("Don't know what to do with element size %i\n", (int)elsize); cbf_free_handle(cbfh); - return 1; + return NULL; + } + + if ( !elsigned ) { + STATUS("Don't know what to do with unsigned data (yet)\n"); + cbf_free_handle(cbfh); + return NULL; } if ( strcmp(byteorder, "little_endian") != 0 ) { STATUS("Don't know what to do with non-little-endian datan\n"); cbf_free_handle(cbfh); - return 1; + return NULL; } data = malloc(elsize*dimfast*dimmid); if ( data == NULL ) { ERROR("Failed to allocate memory for CBF data\n"); cbf_free_handle(cbfh); - return 1; + return NULL; } r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1, @@ -682,10 +762,44 @@ static int read_cbf(struct imagefile *f, struct image *image) ERROR("Failed to read CBF array: %s\n", err); free(err); cbf_free_handle(cbfh); + return NULL; + } + + if ( elsize == 4 ) { + dataf = convert_sint32_float(data, dimfast, dimmid); + } else if ( elsize == 2 ) { + dataf = convert_sint16_float(data, dimfast, dimmid); + } else { + ERROR("Don't know how to convert element size %i\n", + (int)elsize); + cbf_free_handle(cbfh); + return NULL; + } + + free(data); + + free(buf); /* Might be NULL */ + + *w = dimfast; + *h = dimmid; + *pcbfh = cbfh; + return dataf; +} + + +static int read_cbf(struct imagefile *f, struct image *image) +{ + float *data; + int w, h; + cbf_handle cbfh; + + data = read_cbf_data(f, &w, &h, &cbfh); + if ( data == NULL ) { + ERROR("Failed to read CBF data\n"); return 1; } - unpack_panels(image, data, dimfast, dimmid); + unpack_panels(image, data, w, h); free(data); if ( image->beam != NULL ) { @@ -704,126 +818,26 @@ static int read_cbf(struct imagefile *f, struct image *image) } -static float *convert_float(signed int *data, int w, int h) -{ - float *df; - long int i; - - df = malloc(sizeof(float)*w*h); - if ( df == NULL ) return NULL; - - for ( i=0; i<w*h; i++ ) { - df[i] = data[i]; - } - - return df; -} - - static int read_cbf_simple(struct imagefile *f, struct image *image) { + float *data; + int w, h; cbf_handle cbfh; - FILE *fh; - int r; - unsigned int compression; - int binary_id, minelement, maxelement, elsigned, elunsigned; - size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding; - const char *byteorder; - signed int *data; - - if ( cbf_make_handle(&cbfh) ) { - ERROR("Failed to allocate CBF handle\n"); - return 1; - } - - fh = fopen(f->filename, "rb"); - if ( fh == NULL ) { - ERROR("Failed to open '%s'\n", f->filename); - return 1; - } - /* CBFlib calls fclose(fh) when it's ready */ - - if ( cbf_read_widefile(cbfh, fh, 0) ) { - ERROR("Failed to read CBF file '%s'\n", f->filename); - cbf_free_handle(cbfh); - return 1; - } - - /* Select row 0 in data column inside array_data */ - cbf_find_category(cbfh, "array_data"); - cbf_find_column(cbfh, "data"); - cbf_select_row(cbfh, 0); - - /* Get parameters for array read */ - r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id, - &elsize, &elsigned, &elunsigned, - &elements, - &minelement, &maxelement, - &byteorder, - &dimfast, &dimmid, &dimslow, - &padding); - if ( r ) { - char *err = cbf_strerr(r); - ERROR("Failed to read CBF array parameters: %s\n", err); - free(err); - cbf_free_handle(cbfh); - return 1; - } - if ( dimslow != 0 ) { - ERROR("CBF data array is 3D - don't know what to do with it\n"); - cbf_free_handle(cbfh); - return 1; - } - - if ( dimfast*dimmid*elsize > 10e9 ) { - ERROR("CBF data is far too big (%i x %i x %i bytes).\n", - (int)dimfast, (int)dimmid, (int)elsize); - cbf_free_handle(cbfh); - return 1; - } - - if ( elsize != 4 ) { - STATUS("Don't know what to do with element size %i\n", - (int)elsize); - cbf_free_handle(cbfh); - return 1; - } - - if ( strcmp(byteorder, "little_endian") != 0 ) { - STATUS("Don't know what to do with non-little-endian datan\n"); - cbf_free_handle(cbfh); - return 1; - } - - data = malloc(elsize*dimfast*dimmid); + data = read_cbf_data(f, &w, &h, &cbfh); if ( data == NULL ) { - ERROR("Failed to allocate memory for CBF data\n"); - cbf_free_handle(cbfh); + ERROR("Failed to read CBF data\n"); return 1; } - r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1, - elsize*dimfast*dimmid, &elread); - if ( r ) { - char *err = cbf_strerr(r); - ERROR("Failed to read CBF array: %s\n", err); - free(err); - cbf_free_handle(cbfh); - return 1; - } - - image->det = simple_geometry(image, dimfast, dimmid); + image->det = simple_geometry(image, w, h); image->dp = malloc(sizeof(float *)); if ( image->dp == NULL ) { ERROR("Failed to allocate dp array\n"); return 1; } - image->dp[0] = convert_float(data, dimfast, dimmid); - if ( image->dp[0] == NULL ) { - ERROR("Failed to allocate dp array\n"); - return 1; - } + + image->dp[0] = data; if ( image->beam != NULL ) { cbf_fill_in_beam_parameters(image->beam, f, image); @@ -840,6 +854,7 @@ static int read_cbf_simple(struct imagefile *f, struct image *image) return 0; } + #else /* HAVE_CBFLIB */ static int read_cbf_simple(struct imagefile *f, struct image *image) @@ -881,6 +896,24 @@ signed int is_cbf_file(const char *filename) } +signed int is_cbfgz_file(const char *filename) +{ + gzFile gzfh; + char line[1024]; + + gzfh = gzopen(filename, "rb"); + if ( gzfh == NULL ) return -1; + if ( gzgets(gzfh, line, 1024) == NULL ) return -1; + gzclose_r(gzfh); + + if ( strstr(line, "CBF") == NULL ) { + return 0; + } + + return 1; +} + + struct imagefile *imagefile_open(const char *filename) { struct imagefile *f; @@ -903,6 +936,15 @@ struct imagefile *imagefile_open(const char *filename) f->type = IMAGEFILE_CBF; + } else if ( is_cbfgz_file(filename) ) { + + f->type = IMAGEFILE_CBFGZ; + + } else { + + ERROR("Unrecognised file type: %s\n", filename); + return NULL; + } f->filename = strdup(filename); @@ -917,6 +959,8 @@ int imagefile_read(struct imagefile *f, struct image *image, return hdf5_read2(f->hdfile, image, event, 0); } else if ( f->type == IMAGEFILE_CBF ) { return read_cbf(f, image); + } else if ( f->type == IMAGEFILE_CBFGZ ) { + return read_cbf(f, image); } else { ERROR("Unknown file type %i\n", f->type); return 1; @@ -930,8 +974,13 @@ int imagefile_read_simple(struct imagefile *f, struct image *image) { if ( f->type == IMAGEFILE_HDF5 ) { return hdf5_read(f->hdfile, image, NULL, 0); - } else { + } else if ( f->type == IMAGEFILE_CBF ) { + return read_cbf_simple(f, image); + } else if ( f->type == IMAGEFILE_CBFGZ ) { return read_cbf_simple(f, image); + } else { + ERROR("Unknown file type %i\n", f->type); + return 1; } } diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 40717ab3..cc1bcc5d 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -60,6 +60,7 @@ struct imagefile_field_list; * @SPECTRUM_TOPHAT: A top hat distribution of wavelengths * @SPECTRUM_SASE: A simulated SASE spectrum * @SPECTRUM_TWOCOLOUR: A spectrum containing two peaks + * @SPECTRUM_FROMFILE: An arbitrary spectrum read from input file * * A %SpectrumType represents a type of X-ray energy spectrum to use for * generating simulated data. @@ -67,7 +68,8 @@ struct imagefile_field_list; typedef enum { SPECTRUM_TOPHAT, SPECTRUM_SASE, - SPECTRUM_TWOCOLOUR + SPECTRUM_TWOCOLOUR, + SPECTRUM_FROMFILE } SpectrumType; @@ -114,17 +116,27 @@ struct imagefeature { enum imagefile_type { IMAGEFILE_HDF5, - IMAGEFILE_CBF + IMAGEFILE_CBF, + IMAGEFILE_CBFGZ }; /* An opaque type representing a list of image features */ typedef struct _imagefeaturelist ImageFeatureList; + +struct spectrum +{ + int n; + double *ks; /* 1/m */ + double *weights; +}; + + /* Structure describing a wavelength sample from a spectrum */ struct sample { - double k; + double k; /* 1/m */ double weight; }; @@ -145,9 +157,11 @@ struct beam_params /** * image: + * @hit: Non-zero if the frame was determined to be a "hit" * @crystals: Array of crystals in the image * @n_crystals: The number of crystals in the image * @indexed_by: Indexing method which indexed this pattern + * @n_indexing_tries: Number of times the indexer was tried before indexing * @det: Detector structure * @beam: Beam parameters structure * @filename: Filename for the image file @@ -168,7 +182,7 @@ struct beam_params * @avg_clen: Mean of camera length values for all panels * @spectrum: Spectrum information * @nsamples: Number of spectrum samples - * @spectrum_size: SIze of spectrum array + * @spectrum_size: Size of spectrum array * * The field <structfield>data</structfield> contains the raw image data, if it * is currently available. The data might be available throughout the @@ -190,9 +204,11 @@ struct image { int **bad; /* Bad pixels by panel */ float **sat; /* Per-pixel saturation values */ + int hit; Crystal **crystals; int n_crystals; IndexingMethod indexed_by; + int n_indexing_tries; struct detector *det; struct beam_params *beam; /* The nominal beam parameters */ @@ -209,7 +225,11 @@ struct image { int serial; /* Monotonically ascending serial * number for this image */ - struct sample *spectrum; + struct spectrum *spectrum; /* Beam spectrum for pink beam data */ + + /* FIXME: These are only used in pattern_sim, which should be changed to + * use the "struct spectrum" from above later */ + struct sample *spectrum0; int nsamples; /* Number of wavelengths */ int spectrum_size; /* Size of "spectrum" */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index ff84c629..1476304d 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -57,6 +57,7 @@ #include "felix.h" #include "predict-refine.h" #include "taketwo.h" +#include "xgandalf.h" struct _indexingprivate @@ -66,6 +67,7 @@ struct _indexingprivate float tolerance[4]; struct taketwo_options *ttopts; + struct xgandalf_options *xgandalf_opts; int n_methods; IndexingMethod *methods; @@ -182,6 +184,10 @@ static char *base_indexer_str(IndexingMethod indm) strcpy(str, "taketwo"); break; + case INDEXING_XGANDALF: + strcpy(str, "xgandalf"); + break; + case INDEXING_SIMULATION : strcpy(str, "simulation"); break; @@ -219,6 +225,7 @@ static char *friendly_indexer_name(IndexingMethod m) static void *prepare_method(IndexingMethod *m, UnitCell *cell, + struct xgandalf_options *xgandalf_opts, struct felix_options *felix_opts) { char *str; @@ -259,6 +266,10 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell, priv = taketwo_prepare(m, cell); break; + case INDEXING_XGANDALF : + priv = xgandalf_prepare(m, cell, xgandalf_opts); + break; + default : ERROR("Don't know how to prepare indexing method %i\n", *m); break; @@ -290,6 +301,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, struct detector *det, float *ltl, IndexingFlags flags, struct taketwo_options *ttopts, + struct xgandalf_options *xgandalf_opts, struct felix_options *felix_opts) { int i, n; @@ -386,6 +398,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, int j; ipriv->engine_private[i] = prepare_method(&methods[i], cell, + xgandalf_opts, felix_opts); if ( ipriv->engine_private[i] == NULL ) return NULL; @@ -411,6 +424,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i]; ipriv->ttopts = ttopts; + ipriv->xgandalf_opts = xgandalf_opts; STATUS("List of indexing methods:\n"); for ( i=0; i<n; i++ ) { @@ -467,6 +481,10 @@ void cleanup_indexing(IndexingPrivate *ipriv) taketwo_cleanup(ipriv->engine_private[n]); break; + case INDEXING_XGANDALF : + xgandalf_cleanup(ipriv->engine_private[n]); + break; + default : ERROR("Don't know how to clean up indexing method %i\n", ipriv->methods[n]); @@ -577,6 +595,10 @@ static int try_indexer(struct image *image, IndexingMethod indm, r = taketwo_index(image, ipriv->ttopts, mpriv); break; + case INDEXING_XGANDALF : + r = run_xgandalf(image, mpriv); + break; + default : ERROR("Unrecognised indexing method: %i\n", indm); return 0; @@ -584,7 +606,11 @@ static int try_indexer(struct image *image, IndexingMethod indm, } /* Stop a really difficult to debug situation in its tracks */ - assert(image->n_crystals - n_before == r); + if ( image->n_crystals - n_before != r ) { + ERROR("Whoops, indexer didn't return the right number " + "of crystals!\n"); + exit(1); + } /* For all the crystals found this time ... */ for ( i=0; i<r; i++ ) { @@ -649,7 +675,7 @@ static int try_indexer(struct image *image, IndexingMethod indm, if ( compare_cells(crystal_get_cell(cr), crystal_get_cell(that_cr), - 0.1, deg2rad(5.0), NULL) ) + 0.1, deg2rad(0.5), NULL) ) { crystal_set_user_flag(cr, 1); } @@ -829,7 +855,10 @@ void index_pattern_2(struct image *image, IndexingPrivate *ipriv, int *ping) /* Stop now if the pattern is indexed (don't try again for more * crystals with a different indexing method) */ - if ( success ) break; + if ( success ) { + image->n_indexing_tries = ntry; + break; + } } @@ -945,6 +974,11 @@ IndexingMethod get_indm_from_string_2(const char *str, int *err) method = INDEXING_DEFAULTS_TAKETWO; have_method = 1; + } else if ( strcmp(bits[i], "xgandalf") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_DEFAULTS_XGANDALF; + have_method = 1; + } else if ( strcmp(bits[i], "none") == 0) { if ( have_method ) return warn_method(str); method = INDEXING_NONE; @@ -1042,6 +1076,7 @@ char *detect_indexing_methods(UnitCell *cell) do_probe(dirax_probe, cell, methods); do_probe(asdf_probe, cell, methods); do_probe(xds_probe, cell, methods); + do_probe(xgandalf_probe, cell, methods); /* Don't automatically use TakeTwo or Felix (yet) */ //do_probe(taketwo_probe, cell, methods); //do_probe(felix_probe, cell, methods); diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 8ece9227..2099b4d7 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -56,6 +56,8 @@ #define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE \ | INDEXING_USE_CELL_PARAMETERS) +#define INDEXING_DEFAULTS_XGANDALF (INDEXING_XGANDALF | INDEXING_USE_CELL_PARAMETERS) + /** * IndexingMethod: * @INDEXING_NONE: No indexing to be performed @@ -67,6 +69,7 @@ * @INDEXING_DEBUG: Results injector for debugging * @INDEXING_ASDF: Use in-built "asdf" indexer * @INDEXING_TAKETWO: Use in-built "taketwo" indexer + * @INDEXING_XGANDALF: Invoke XGANDALF * @INDEXING_ERROR: Special value for unrecognised indexing engine name * @INDEXING_USE_LATTICE_TYPE: Use lattice type and centering information to * guide the indexing process. @@ -90,6 +93,7 @@ typedef enum { INDEXING_DEBUG = 7, INDEXING_ASDF = 8, INDEXING_TAKETWO = 9, + INDEXING_XGANDALF = 10, INDEXING_ERROR = 255, /* Unrecognised indexing engine */ @@ -136,6 +140,7 @@ extern IndexingMethod get_indm_from_string_2(const char *method, int *err); #include "cell.h" #include "image.h" #include "taketwo.h" +#include "xgandalf.h" #include "felix.h" @@ -143,6 +148,7 @@ extern IndexingPrivate *setup_indexing(const char *methods, UnitCell *cell, struct detector *det, float *ltl, IndexingFlags flags, struct taketwo_options *ttopts, + struct xgandalf_options *xgandalf_opts, struct felix_options *felix_opts); extern char *detect_indexing_methods(UnitCell *cell); diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index 9dc99162..d48dcef2 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -570,18 +570,20 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) cell_get_parameters(mosflm->mp->template, &a, &b, &c, &alpha, &beta, &gamma); cen = cell_get_centering(mosflm->mp->template); + /* Old MOSFLM simply ignores CELL and CENTERING subkeywords. + * So this doesn't do any harm. + * TODO: but still better to show WARNING if MOSFLM is old. */ + snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 " + "MAXCELL 1000 REFINE " + "CELL %.1f %.1f %.1f %.1f %.1f %.1f " + "CENTERING %c\n", + mosflm->sptfile, a*1e10, b*1e10, c*1e10, + rad2deg(alpha), rad2deg(beta), rad2deg(gamma), + cen); } else { - cen = 'P'; - a = 0; /* Disables prior-cell algorithm in MOSFLM */ + snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 " + "MAXCELL 1000 REFINE\n", mosflm->sptfile); } - /* Old MOSFLM simply ignores CELL and CENTERING subkeywords. - * So this doesn't do any harm. - * TODO: but still better to show WARNING if MOSFLM is old. */ - snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 MAXCELL 1000 " - "REFINE " - "CELL %.1f %.1f %.1f %.1f %.1f %.1f CENTERING %c\n", - mosflm->sptfile, a*1e10, b*1e10, c*1e10, - rad2deg(alpha), rad2deg(beta), rad2deg(gamma), cen); mosflm_sendline(tmp, mosflm); break; @@ -844,6 +846,15 @@ void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell) *indm &= INDEXING_METHOD_MASK | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS; + if ( (cell != NULL) && (cell_get_centering(cell) == 'I') + && (cell_get_lattice_type(cell) == L_MONOCLINIC) ) + { + ERROR("WARNING: Mosflm always gives the monoclinic C cell in " + "preference to the monoclinic I cell choice.\n"); + ERROR("To get a higher indexing rate, convert your cell to the " + "monoclinic C cell choice.\n"); + } + mp = malloc(sizeof(struct mosflm_private)); if ( mp == NULL ) return NULL; @@ -875,20 +886,6 @@ static void chop_word(char *s) } -static int file_exists(const char *filename) -{ - struct stat s; - - if ( stat(filename, &s) != 0 ) { - if ( errno == ENOENT ) return 0; - ERROR("Failed to check for %s.\n", filename); - exit(1); - } - - return 1; -} - - const char *mosflm_probe(UnitCell *cell) { pid_t pid; @@ -899,15 +896,6 @@ const char *mosflm_probe(UnitCell *cell) int ok = 0; int l; - /* Mosflm will write mosflm.lp and SUMMARY when we test it, which we are - * are going to delete afterwards. Better check they don't exist first, - * in case they were important. */ - if ( file_exists("mosflm.lp") || file_exists("SUMMARY") ) { - ERROR("Please move or delete mosflm.lp and SUMMARY from the " - "working directory first.\n"); - exit(1); - } - pid = forkpty(&pty, NULL, NULL, NULL); if ( pid == -1 ) { return NULL; @@ -947,9 +935,6 @@ const char *mosflm_probe(UnitCell *cell) close(pty); waitpid(pid, &status, 0); - unlink("mosflm.lp"); - unlink("SUMMARY"); - if ( !ok ) return NULL; if ( cell_has_parameters(cell) ) return "mosflm-cell-nolatt,mosflm-latt-nocell"; diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 39fe8a79..9090d48f 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -298,35 +298,35 @@ static struct radial_stats* allocate_radial_stats(int num_rad_bins) rstats->rthreshold = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->rthreshold == NULL ) { - free(rstats); free(rstats->roffset); + free(rstats); return NULL; } rstats->lthreshold = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->lthreshold == NULL ) { - free(rstats); free(rstats->rthreshold); free(rstats->roffset); + free(rstats); return NULL; } rstats->rsigma = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->rsigma == NULL ) { - free(rstats); free(rstats->roffset); free(rstats->rthreshold); free(rstats->lthreshold); + free(rstats); return NULL; } rstats->rcount = (int *)malloc(num_rad_bins*sizeof(int)); if ( rstats->rcount == NULL ) { - free(rstats); free(rstats->roffset); free(rstats->rthreshold); free(rstats->lthreshold); free(rstats->rsigma); + free(rstats); return NULL; } @@ -434,8 +434,8 @@ struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) pkdata->npix = (int *)malloc(max_num_peaks*sizeof(int)); if ( pkdata->npix == NULL ) { - free(pkdata); free(pkdata->npix); + free(pkdata); return NULL; } diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 61a94e05..cccd5d0b 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -319,7 +319,7 @@ static int pair_peaks(struct image *image, Crystal *cr, } -void refine_radius(Crystal *cr, struct image *image) +int refine_radius(Crystal *cr, struct image *image) { int n, n_acc; struct reflpeak *rps; @@ -328,13 +328,13 @@ void refine_radius(Crystal *cr, struct image *image) /* Maximum possible size */ rps = malloc(image_feature_count(image->features) * sizeof(struct reflpeak)); - if ( rps == NULL ) return; + if ( rps == NULL ) return 1; reflist = reflist_new(); n_acc = pair_peaks(image, cr, reflist, rps); if ( n_acc < 3 ) { free(rps); - return; + return 1; } crystal_set_reflections(cr, reflist); update_predictions(cr); @@ -347,6 +347,8 @@ void refine_radius(Crystal *cr, struct image *image) reflist_free(reflist); free(rps); + + return 0; } diff --git a/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h index c763d2ca..fe700f47 100644 --- a/libcrystfel/src/predict-refine.h +++ b/libcrystfel/src/predict-refine.h @@ -39,7 +39,7 @@ struct image; extern int refine_prediction(struct image *image, Crystal *cr); -extern void refine_radius(Crystal *cr, struct image *image); +extern int refine_radius(Crystal *cr, struct image *image); #endif /* PREDICT_REFINE_H */ diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c index 54c467b3..70548994 100644 --- a/libcrystfel/src/reflist-utils.c +++ b/libcrystfel/src/reflist-utils.c @@ -3,11 +3,11 @@ * * Utilities to complement the core reflist.c * - * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2011-2017 Thomas White <taw@physics.org> + * 2011-2018 Thomas White <taw@physics.org> * 2014 Valerio Mariani * * This file is part of CrystFEL. @@ -648,6 +648,30 @@ RefList *copy_reflist(RefList *list) } +/** + * free_contribs: + * @list: A %RefList + * + * Goes through @list and frees all the reflection contribution structures. + **/ +void free_contribs(RefList *list) +{ + Reflection *refl; + RefListIterator *iter; + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + struct reflection_contributions *c; + c = get_contributions(refl); + free(c->contribs); + free(c->contrib_crystals); + free(c); + } +} + + static char *full_command_line(int argc, char *argv[]) { int i; diff --git a/libcrystfel/src/reflist-utils.h b/libcrystfel/src/reflist-utils.h index f64e9f51..c955491a 100644 --- a/libcrystfel/src/reflist-utils.h +++ b/libcrystfel/src/reflist-utils.h @@ -3,11 +3,11 @@ * * Utilities to complement the core reflist.c * - * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2011-2017 Thomas White <taw@physics.org> + * 2011-2018 Thomas White <taw@physics.org> * 2014 Valerio Mariani * * This file is part of CrystFEL. @@ -69,6 +69,8 @@ extern RefList *res_cutoff(RefList *list, UnitCell *cell, extern RefList *copy_reflist(RefList *list); +extern void free_contribs(RefList *list); + extern void reflist_add_command_and_version(RefList *list, int argcv, char *argv[]); diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c index a35aa575..f1d759ea 100644 --- a/libcrystfel/src/reflist.c +++ b/libcrystfel/src/reflist.c @@ -3,11 +3,11 @@ * * Fast reflection/peak list * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2011-2016 Thomas White <taw@physics.org> + * 2011-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -100,6 +100,9 @@ struct _refldata { double peak; double mean_bg; + /* Contributions */ + struct reflection_contributions *contribs; + /* User-specified temporary values */ double temp1; double temp2; @@ -587,6 +590,18 @@ int get_flag(const Reflection *refl) } +/** + * get_contributions: + * @refl: A %Reflection + * + * Returns: the reflection's contribution list + * + **/ +struct reflection_contributions *get_contributions(const Reflection *refl) +{ + return refl->data.contribs; +} + /********************************** Setters ***********************************/ /** @@ -634,6 +649,7 @@ void set_panel(Reflection *refl, struct panel *p) refl->data.panel = p; } + /** * set_khalf: * @refl: A %Reflection @@ -647,7 +663,6 @@ void set_khalf(Reflection *refl, double khalf) } - /** * set_kpred: * @refl: A %Reflection @@ -848,6 +863,21 @@ void set_flag(Reflection *refl, int flag) } +/** + * set_contributions: + * @refl: A %Reflection + * @contribs: Pointer to the contribution list + * + * Note that the pointer will be stored, not the contents of the structure. + * + **/ +void set_contributions(Reflection *refl, + struct reflection_contributions *contribs) +{ + refl->data.contribs = contribs; +} + + /********************************* Insertion **********************************/ static Reflection *rotate_once(Reflection *refl, int dir) diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h index 9f494474..30bbaa28 100644 --- a/libcrystfel/src/reflist.h +++ b/libcrystfel/src/reflist.h @@ -3,11 +3,11 @@ * * Fast reflection/peak list * - * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2011-2017 Thomas White <taw@physics.org> + * 2011-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -38,6 +38,7 @@ #define GET_K(serial) ((((serial) & 0x000ffc00)>>10)-512) #define GET_L(serial) (((serial) & 0x000003ff)-512) + /** * RefList: * @@ -69,10 +70,21 @@ typedef struct _reflection Reflection; **/ typedef struct _reflistiterator RefListIterator; +#include "crystal.h" + #ifdef __cplusplus extern "C" { #endif +/* Structure representing the contributions to a merged reflection */ +struct reflection_contributions +{ + int n_contrib; + int max_contrib; + Reflection **contribs; + Crystal **contrib_crystals; +}; + /* Creation/deletion */ extern RefList *reflist_new(void); extern void reflist_free(RefList *list); @@ -105,6 +117,7 @@ extern double get_phase(const Reflection *refl, int *have_phase); extern double get_peak(const Reflection *refl); extern double get_mean_bg(const Reflection *refl); extern int get_flag(const Reflection *refl); +extern struct reflection_contributions *get_contributions(const Reflection *refl); /* Set */ extern void copy_data(Reflection *to, const Reflection *from); @@ -126,6 +139,8 @@ extern void set_mean_bg(Reflection *refl, double mean_bg); extern void set_symmetric_indices(Reflection *refl, signed int hs, signed int ks, signed int ls); extern void set_flag(Reflection *refl, int flag); +extern void set_contributions(Reflection *refl, + struct reflection_contributions *contribs); /* Insertion */ extern Reflection *add_refl(RefList *list, diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index c11e38df..bcd2627e 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -832,9 +832,13 @@ int write_chunk(Stream *st, struct image *i, struct imagefile *imfile, fprintf(st->fh, "Image serial number: %i\n", i->serial); + fprintf(st->fh, "hit = %i\n", i->hit); indexer = indexer_str(i->indexed_by); fprintf(st->fh, "indexed_by = %s\n", indexer); free(indexer); + if ( i->indexed_by != INDEXING_NONE ) { + fprintf(st->fh, "n_indexing_tries = %i\n", i->n_indexing_tries); + } fprintf(st->fh, "photon_energy_eV = %f\n", J_to_eV(ph_lambda_to_en(i->lambda))); @@ -1355,6 +1359,7 @@ static void read_audit_lines(Stream *st) ERROR("Failed to allocate memory for audit information\n"); return; } + st->audit_info[0] = '\0'; /* Read lines from stream until one of them starts with "-----", * then rewind to the start of that line */ diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 5f3b67f2..06cc368c 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -36,6 +36,7 @@ #include <string.h> #include <math.h> #include <assert.h> +#include <ctype.h> #include "symmetry.h" #include "utils.h" @@ -1656,19 +1657,45 @@ static IntegerMatrix *parse_symmetry_operation(const char *s) int c; size_t cl; - int sign = +1; - int nh = 0; - int nk = 0; - int nl = 0; - + signed int nh = 0; + signed int nk = 0; + signed int nl = 0; + signed int mult = 1; + int ndigit = 0; + signed int sign = +1; + + /* We have one expression something like "-2h+k" */ cl = strlen(els[i]); - for ( c=0; c<cl; c++ ) { - if ( els[i][c] == '-' ) sign = -1; - if ( els[i][c] == '+' ) sign = +1; - if ( els[i][c] == 'h' ) nh += sign; - if ( els[i][c] == 'k' ) nk += sign; - if ( els[i][c] == 'l' ) nl += sign; + + if ( els[i][c] == '-' ) sign *= -1; + if ( els[i][c] == 'h' ) { + nh = mult*sign; + mult = 1; + ndigit = 0; + sign = +1; + } + if ( els[i][c] == 'k' ) { + nk = mult*sign; + mult = 1; + ndigit = 0; + sign = +1; + } + if ( els[i][c] == 'l' ) { + nl = mult*sign; + mult = 1; + ndigit = 0; + sign = +1; + } + if ( isdigit(els[i][c]) ) { + if ( ndigit > 0 ) { + mult *= 10; + mult += els[i][c] - '0'; + } else { + mult *= els[i][c] - '0'; + } + ndigit++; + } } intmat_set(m, i, 0, nh); diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index a759ff15..2c756ad7 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -185,6 +185,8 @@ static inline int within_tolerance(double a, double b, double percent) /* Photon energy (eV) to wavelength (m) */ #define ph_eV_to_lambda(a) ph_en_to_lambda(eV_to_J(a)) +/* Photon energy (eV) to k (1/m) */ +#define ph_eV_to_k(a) ((a)*ELECTRON_CHARGE/PLANCK/C_VACUO) /* ------------------------------ Message macros ---------------------------- */ diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c index 01c12dca..cef619c2 100644 --- a/libcrystfel/src/xds.c +++ b/libcrystfel/src/xds.c @@ -73,18 +73,43 @@ struct xds_private }; +/* Essentially the reverse of spacegroup_for_lattice(), below */ +static int convert_spacegroup_number(int spg, LatticeType *lt, char *cen, + char *ua) +{ + switch ( spg ) { + + case 1: *lt = L_TRICLINIC; *cen = 'P'; *ua = '*'; return 0; + case 3: *lt = L_MONOCLINIC; *cen = 'P'; *ua = 'b'; return 0; + case 5: *lt = L_MONOCLINIC; *cen = 'C'; *ua = 'b'; return 0; + case 16: *lt = L_ORTHORHOMBIC; *cen = 'P'; *ua = '*'; return 0; + case 21: *lt = L_ORTHORHOMBIC; *cen = 'C'; *ua = '*'; return 0; + case 22: *lt = L_ORTHORHOMBIC; *cen = 'F'; *ua = '*'; return 0; + case 23: *lt = L_ORTHORHOMBIC; *cen = 'I'; *ua = '*'; return 0; + case 75: *lt = L_TETRAGONAL; *cen = 'P'; *ua = 'c'; return 0; + case 79: *lt = L_TETRAGONAL; *cen = 'I'; *ua = 'c'; return 0; + case 143: *lt = L_HEXAGONAL; *cen = 'P'; *ua = 'c'; return 0; + case 146: *lt = L_HEXAGONAL; *cen = 'H'; *ua = 'c'; return 0; + case 195: *lt = L_CUBIC; *cen = 'P'; *ua = '*'; return 0; + case 196: *lt = L_CUBIC; *cen = 'F'; *ua = '*'; return 0; + case 197: *lt = L_CUBIC; *cen = 'I'; *ua = '*'; return 0; + default: return 1; + + } +} + + static int read_cell(struct image *image) { FILE * fh; - float axstar, aystar, azstar; - float bxstar, bystar, bzstar; - float cxstar, cystar, czstar; - char asx[11], asy[11], asz[11]; - char bsx[11], bsy[11], bsz[11]; - char csx[11], csy[11], csz[11]; + float ax, ay, az; + float bx, by, bz; + float cx, cy, cz; + int spg; char *rval, line[1024]; - int r; UnitCell *cell; + LatticeType latticetype; + char centering, ua; Crystal *cr; fh = fopen("IDXREF.LP", "r"); @@ -97,74 +122,66 @@ static int read_cell(struct image *image) return 0; } - } while ( strcmp(line, " # COORDINATES OF REC. BASIS VECTOR" - " LENGTH 1/LENGTH\n") != 0 ); + } while ( strcmp(line, " ***** DIFFRACTION PARAMETERS USED AT START OF " + "INTEGRATION *****\n") != 0 ); - /* Free line after chunk */ - rval = fgets(line, 1023, fh); - if ( rval == NULL ) { + /* Find and read space group number */ + do { + rval = fgets(line, 1023, fh); + if ( rval == NULL ) { + fclose(fh); + return 0; + } + } while ( strncmp(line, " SPACE GROUP NUMBER ", 20) != 0 ); + sscanf(line+20, "%i\n", &spg); + + /* Find and read a */ + do { + rval = fgets(line, 1023, fh); + if ( rval == NULL ) { + fclose(fh); + return 0; + } + } while ( strncmp(line, " COORDINATES OF UNIT CELL A-AXIS ", 33) != 0 ); + if ( sscanf(line+33, "%f %f %f\n", &ax, &ay, &az) < 3 ) { fclose(fh); return 0; } - /* Get first vector */ + /* Read b */ rval = fgets(line, 1023, fh); if ( rval == NULL ) { fclose(fh); return 0; } - if ( line[4] != '1' ) { - ERROR("No first vector from XDS.\n"); + if ( sscanf(line+33, "%f %f %f\n", &bx, &by, &bz) < 3 ) { + fclose(fh); return 0; } - memcpy(asx, line+7, 10); asx[10] = '\0'; - memcpy(asy, line+17, 10); asy[10] = '\0'; - memcpy(asz, line+27, 10); asz[10] = '\0'; - /* Get second vector */ + /* Read c */ rval = fgets(line, 1023, fh); if ( rval == NULL ) { fclose(fh); return 0; } - if ( line[4] != '2' ) { - ERROR("No second vector from XDS.\n"); + if ( sscanf(line+33, "%f %f %f\n", &cx, &cy, &cz) < 3 ) { + fclose(fh); return 0; } - memcpy(bsx, line+7, 10); bsx[10] = '\0'; - memcpy(bsy, line+17, 10); bsy[10] = '\0'; - memcpy(bsz, line+27, 10); bsz[10] = '\0'; - /* Get third vector */ - rval = fgets(line, 1023, fh); - fclose(fh); - if ( rval == NULL ) return 0; - if ( line[4] != '3' ) return 0; /* No error message this time - * - happens a lot */ - memcpy(csx, line+7, 10); csx[10] = '\0'; - memcpy(csy, line+17, 10); csy[10] = '\0'; - memcpy(csz, line+27, 10); csz[10] = '\0'; - - r = sscanf(asx, "%f", &cxstar); - r += sscanf(asy, "%f", &cystar); - r += sscanf(asz, "%f", &czstar); - r += sscanf(bsx, "%f", &bxstar); - r += sscanf(bsy, "%f", &bystar); - r += sscanf(bsz, "%f", &bzstar); - r += sscanf(csx, "%f", &axstar); - r += sscanf(csy, "%f", &aystar); - r += sscanf(csz, "%f", &azstar); - - if ( r != 9 ) { - STATUS("Fewer than 9 parameters found in NEWMAT file.\n"); + cell = cell_new(); + cell_set_cartesian(cell, + ax*1e-10, ay*1e-10, az*1e-10, + bx*1e-10, by*1e-10, bz*1e-10, + -cx*1e-10, -cy*1e-10, -cz*1e-10); + if ( convert_spacegroup_number(spg, &latticetype, ¢ering, &ua) ) { + ERROR("Failed to convert XDS space group number (%i)\n", spg); return 0; } - - cell = cell_new(); - cell_set_reciprocal(cell, - axstar*10e9, aystar*10e9, azstar*10e9, - bxstar*10e9, bystar*10e9, bzstar*10e9, - -cxstar*10e9, -cystar*10e9, -czstar*10e9); + cell_set_lattice_type(cell, latticetype); + cell_set_centering(cell, centering); + cell_set_unique_axis(cell, ua); cr = crystal_new(); if ( cr == NULL ) { @@ -235,13 +252,15 @@ static const char *spacegroup_for_lattice(UnitCell *cell) switch ( latt ) { case L_TRICLINIC : - g = "1"; + if ( centering == 'P' ) { + g = "1"; + } break; case L_MONOCLINIC : if ( centering == 'P' ) { g = "3"; - } else { + } else if ( centering == 'C' ) { g = "5"; } break; @@ -250,10 +269,10 @@ static const char *spacegroup_for_lattice(UnitCell *cell) if ( centering == 'P' ) { g = "16"; } else if ( centering == 'C' ) { - g = "20"; + g = "21"; } else if ( centering == 'F' ) { g = "22"; - } else { + } else if ( centering == 'I' ) { g = "23"; } break; @@ -261,34 +280,34 @@ static const char *spacegroup_for_lattice(UnitCell *cell) case L_TETRAGONAL : if ( centering == 'P' ) { g = "75"; - } else { + } else if ( centering == 'I' ) { g = "79"; } break; + /* Unfortunately, XDS only does "hexagonal H" */ case L_RHOMBOHEDRAL : + return NULL; + + case L_HEXAGONAL : if ( centering == 'P' ) { g = "143"; - } else { + } + if ( centering == 'H' ) { g = "146"; } break; - case L_HEXAGONAL : - g = "168"; - break; - case L_CUBIC : if ( centering == 'P' ) { g = "195"; } else if ( centering == 'F' ) { g = "196"; - } else { + } else if ( centering == 'I' ) { g = "197"; } break; } - assert(g != NULL); return g; } @@ -348,7 +367,7 @@ static int write_inp(struct image *image, struct xds_private *xp) fprintf(fh, "DETECTOR= CSPAD\n"); fprintf(fh, "MINIMUM_VALID_PIXEL_VALUE= 1\n"); fprintf(fh, "OVERLOAD= 200000000\n"); - fprintf(fh, "INDEX_ERROR= 0.05\n"); + fprintf(fh, "INDEX_ERROR= 0.2\n"); //fprintf(fh, "INDEX_QUALITY= 0.5\n"); fprintf(fh, "REFINE(IDXREF)= CELL ORIENTATION\n"); //fprintf(fh, "MINIMUM_NUMBER_OF_PIXELS_IN_A_SPOT= 1\n"); @@ -442,6 +461,11 @@ void *xds_prepare(IndexingMethod *indm, UnitCell *cell) return NULL; } + if ( (cell != NULL) && (spacegroup_for_lattice(cell) == NULL) ) { + ERROR("Don't know how to ask XDS for your cell.\n"); + return NULL; + } + xp = calloc(1, sizeof(*xp)); if ( xp == NULL ) return NULL; diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c new file mode 100644 index 00000000..2434ce2e --- /dev/null +++ b/libcrystfel/src/xgandalf.c @@ -0,0 +1,335 @@ +/* + * xgandalf.c + * + * Interface to XGANDALF indexer + * + * Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#include "xgandalf.h" + +#ifdef HAVE_XGANDALF +#include <stdlib.h> + +#include "utils.h" +#include "cell-utils.h" +#include "peaks.h" + +#include "xgandalf/adaptions/crystfel/Lattice.h" +#include "xgandalf/adaptions/crystfel/ExperimentSettings.h" +#include "xgandalf/adaptions/crystfel/IndexerPlain.h" + +struct xgandalf_private_data { + IndexerPlain *indexer; + reciprocalPeaks_1_per_A_t reciprocalPeaks_1_per_A; + + IndexingMethod indm; + UnitCell *cellTemplate; + Lattice_t sampleRealLattice_A; //same as cellTemplate + UnitCellTransformation *uncenteringTransformation; + LatticeTransform_t latticeReductionTransform; +}; + +#define FAKE_DETECTOR_DISTANCE (0.1) +#define FAKE_DETECTOR_RADIUS (0.1) +#define FAKE_BEAM_ENERGY (1) +#define FAKE_DIVERGENCE_ANGLE_DEG (0.05) +#define FAKE_NON_MONOCHROMATICITY (0.005) +#define FAKE_REFLECTION_RADIUS (0.0001) + +#define MAX_ASSEMBLED_LATTICES_COUNT (10) + +static void reduceCell(UnitCell* cell, LatticeTransform_t* appliedReductionTransform); +static void restoreCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform); +static void makeRightHanded(UnitCell* cell); + +int run_xgandalf(struct image *image, void *ipriv) +{ + struct xgandalf_private_data *xgandalf_private_data = (struct xgandalf_private_data*) ipriv; + reciprocalPeaks_1_per_A_t *reciprocalPeaks_1_per_A = &(xgandalf_private_data->reciprocalPeaks_1_per_A); + + int peakCountMax = image_feature_count(image->features); + reciprocalPeaks_1_per_A->peakCount = 0; + for (int i = 0; i < peakCountMax && i < MAX_PEAK_COUNT_FOR_INDEXER; i++) { + struct imagefeature *f; + f = image_get_feature(image->features, i); + if (f == NULL) { + continue; + } + + reciprocalPeaks_1_per_A->coordinates_x[reciprocalPeaks_1_per_A->peakCount] = f->rx * 1e-10; + reciprocalPeaks_1_per_A->coordinates_y[reciprocalPeaks_1_per_A->peakCount] = f->ry * 1e-10; + reciprocalPeaks_1_per_A->coordinates_z[reciprocalPeaks_1_per_A->peakCount] = f->rz * 1e-10; + reciprocalPeaks_1_per_A->peakCount++; + } + + Lattice_t assembledLattices[MAX_ASSEMBLED_LATTICES_COUNT]; + int assembledLatticesCount; + IndexerPlain_index(xgandalf_private_data->indexer, + assembledLattices, + &assembledLatticesCount, + MAX_ASSEMBLED_LATTICES_COUNT, + *reciprocalPeaks_1_per_A, + NULL); + + if (assembledLatticesCount > 0) { //no multi-lattice at the moment + assembledLatticesCount = 1; + } + + int goodLatticesCount = assembledLatticesCount; + for (int i = 0; i < assembledLatticesCount && i < 1; i++) { + reorderLattice(&(xgandalf_private_data->sampleRealLattice_A), + &assembledLattices[i]); + + UnitCell *uc; + uc = cell_new(); + + Lattice_t *l = &assembledLattices[i]; + + cell_set_cartesian(uc, l->ax * 1e-10, l->ay * 1e-10, l->az * 1e-10, + l->bx * 1e-10, l->by * 1e-10, l->bz * 1e-10, + l->cx * 1e-10, l->cy * 1e-10, l->cz * 1e-10); + makeRightHanded(uc); + + if(xgandalf_private_data->cellTemplate != NULL){ + restoreCell(uc, &xgandalf_private_data->latticeReductionTransform); + + UnitCell *new_cell_trans = cell_transform_inverse(uc, xgandalf_private_data->uncenteringTransformation); + cell_free(uc); + uc = new_cell_trans; + + cell_set_lattice_type(new_cell_trans, cell_get_lattice_type(xgandalf_private_data->cellTemplate)); + cell_set_centering(new_cell_trans, cell_get_centering(xgandalf_private_data->cellTemplate)); + cell_set_unique_axis(new_cell_trans, cell_get_unique_axis(xgandalf_private_data->cellTemplate)); + } + + if (validate_cell(uc)) { + STATUS("Problem with returned cell!\n"); + } + + Crystal *cr = crystal_new(); + if (cr == NULL) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + crystal_set_cell(cr, uc); + image_add_crystal(image, cr); + + } + + return goodLatticesCount; +} + +void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, + struct xgandalf_options *xgandalf_opts) +{ + struct xgandalf_private_data *xgandalf_private_data = malloc(sizeof(struct xgandalf_private_data)); + allocReciprocalPeaks(&(xgandalf_private_data->reciprocalPeaks_1_per_A)); + xgandalf_private_data->indm = *indm; + xgandalf_private_data->cellTemplate = NULL; + xgandalf_private_data->uncenteringTransformation = NULL; + + float tolerance = xgandalf_opts->tolerance; + samplingPitch_t samplingPitch = xgandalf_opts->sampling_pitch; + gradientDescentIterationsCount_t gradientDescentIterationsCount = xgandalf_opts->grad_desc_iterations; + + if (*indm & INDEXING_USE_CELL_PARAMETERS) { + + xgandalf_private_data->cellTemplate = cell; + + UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->uncenteringTransformation); + + UnitCell *uc = cell_new_from_cell(primitiveCell); + reduceCell(primitiveCell, &xgandalf_private_data->latticeReductionTransform); + + cell_free(uc); + + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; + int ret = cell_get_reciprocal(primitiveCell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + if (ret != 0) { + ERROR("cell_get_reciprocal did not finish properly!"); + } + + Lattice_t sampleReciprocalLattice_1_per_A = { + .ax = asx * 1e-10, .ay = asy * 1e-10, .az = asz * 1e-10, + .bx = bsx * 1e-10, .by = bsy * 1e-10, .bz = bsz * 1e-10, + .cx = csx * 1e-10, .cy = csy * 1e-10, .cz = csz * 1e-10 }; + + double ax, ay, az, bx, by, bz, cx, cy, cz; + ret = cell_get_cartesian(primitiveCell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + if (ret != 0) { + ERROR("cell_get_cartesian did not finish properly!"); + } + Lattice_t sampleRealLattice_A = { + .ax = ax * 1e10, .ay = ay * 1e10, .az = az * 1e10, + .bx = bx * 1e10, .by = by * 1e10, .bz = bz * 1e10, + .cx = cx * 1e10, .cy = cy * 1e10, .cz = cz * 1e10 }; + xgandalf_private_data->sampleRealLattice_A = sampleRealLattice_A; + + ExperimentSettings *experimentSettings = + ExperimentSettings_new(FAKE_BEAM_ENERGY, + FAKE_DETECTOR_DISTANCE, + FAKE_DETECTOR_RADIUS, + FAKE_DIVERGENCE_ANGLE_DEG, + FAKE_NON_MONOCHROMATICITY, + sampleReciprocalLattice_1_per_A, + tolerance, + FAKE_REFLECTION_RADIUS); + + xgandalf_private_data->indexer = IndexerPlain_new(experimentSettings); + IndexerPlain_setSamplingPitch(xgandalf_private_data->indexer, + samplingPitch); + IndexerPlain_setGradientDescentIterationsCount(xgandalf_private_data->indexer, + gradientDescentIterationsCount); + + if (xgandalf_opts->no_deviation_from_provided_cell) { + IndexerPlain_setRefineWithExactLattice(xgandalf_private_data->indexer, 1); + } + + ExperimentSettings_delete(experimentSettings); + cell_free(primitiveCell); + + } else { + + Lattice_t sampleRealLattice_A = { 0, 0, 0, 0, 0, 0, 0, 0, 0 }; + xgandalf_private_data->sampleRealLattice_A = sampleRealLattice_A; + + ExperimentSettings *experimentSettings = + ExperimentSettings_new_nolatt(FAKE_BEAM_ENERGY, + FAKE_DETECTOR_DISTANCE, + FAKE_DETECTOR_RADIUS, + FAKE_DIVERGENCE_ANGLE_DEG, + FAKE_NON_MONOCHROMATICITY, + xgandalf_opts->minLatticeVectorLength_A, + xgandalf_opts->maxLatticeVectorLength_A, + FAKE_REFLECTION_RADIUS); + + xgandalf_private_data->indexer = IndexerPlain_new(experimentSettings); + IndexerPlain_setSamplingPitch(xgandalf_private_data->indexer, + samplingPitch); + IndexerPlain_setGradientDescentIterationsCount(xgandalf_private_data->indexer, + gradientDescentIterationsCount); + + ExperimentSettings_delete(experimentSettings); + } + + /* Flags that XGANDALF knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_USE_CELL_PARAMETERS; + + return xgandalf_private_data; +} + + +void xgandalf_cleanup(void *pp) +{ + struct xgandalf_private_data *xgandalf_private_data = pp; + + freeReciprocalPeaks(xgandalf_private_data->reciprocalPeaks_1_per_A); + IndexerPlain_delete(xgandalf_private_data->indexer); + if(xgandalf_private_data->uncenteringTransformation != NULL){ + tfn_free(xgandalf_private_data->uncenteringTransformation); + } +} + +static void reduceCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform) +{ + double ax, ay, az, bx, by, bz, cx, cy, cz; + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz }; + + reduceLattice(&l, appliedReductionTransform); + + cell_set_cartesian(cell, l.ax, l.ay, l.az, + l.bx, l.by, l.bz, + l.cx, l.cy, l.cz); + + makeRightHanded(cell); +} + +static void restoreCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform) +{ + + double ax, ay, az, bx, by, bz, cx, cy, cz; + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz }; + + restoreLattice(&l, appliedReductionTransform); + + cell_set_cartesian(cell, l.ax, l.ay, l.az, + l.bx, l.by, l.bz, + l.cx, l.cy, l.cz); + + makeRightHanded(cell); +} + +static void makeRightHanded(UnitCell *cell) +{ + double ax, ay, az, bx, by, bz, cx, cy, cz; + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + if ( !right_handed(cell) ) { + cell_set_cartesian(cell, -ax, -ay, -az, -bx, -by, -bz, -cx, -cy, -cz); + } +} + + +const char *xgandalf_probe(UnitCell *cell) +{ + return "xgandalf"; +} + +#else + +int run_xgandalf(struct image *image, void *ipriv) +{ + ERROR("This copy of CrystFEL was compiled without XGANDALF support.\n"); + return 0; +} + + +void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, + struct xgandalf_options *xgandalf_opts) +{ + ERROR("This copy of CrystFEL was compiled without XGANDALF support.\n"); + ERROR("To use XGANDALF indexing, recompile with XGANDALF.\n"); + return NULL; +} + + +void xgandalf_cleanup(void *pp) +{ +} + + +const char *xgandalf_probe(UnitCell *cell) +{ + return NULL; +} + +#endif // HAVE_XGANDALF diff --git a/libcrystfel/src/xgandalf.h b/libcrystfel/src/xgandalf.h new file mode 100644 index 00000000..1aced417 --- /dev/null +++ b/libcrystfel/src/xgandalf.h @@ -0,0 +1,58 @@ +/* + * xgandalf.h + * + * Interface to XGANDALF indexer + * + * Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifndef LIBCRYSTFEL_SRC_XGANDALF_H +#define LIBCRYSTFEL_SRC_XGANDALF_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stddef.h> + +struct xgandalf_options { + unsigned int sampling_pitch; + unsigned int grad_desc_iterations; + float tolerance; + unsigned int no_deviation_from_provided_cell; + float minLatticeVectorLength_A; + float maxLatticeVectorLength_A; +}; + +#include "index.h" + +extern int run_xgandalf(struct image *image, void *ipriv); + +extern void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, + struct xgandalf_options *xgandalf_opts); + +extern void xgandalf_cleanup(void *pp); +extern const char *xgandalf_probe(UnitCell *cell); + + +#endif /* LIBCRYSTFEL_SRC_XGANDALF_H */ diff --git a/scripts/ave-resolution b/scripts/ave-resolution index 69c21195..370b3889 100755 --- a/scripts/ave-resolution +++ b/scripts/ave-resolution @@ -14,7 +14,11 @@ import sys import numpy import matplotlib.pyplot as plt -f = open(sys.argv[1]) +if sys.argv[1] == '-': + f = sys.stdin +else: + f = open(sys.argv[1]) + a = [] while True: diff --git a/scripts/euxfel-train-analysis b/scripts/euxfel-train-analysis new file mode 100755 index 00000000..0d20decd --- /dev/null +++ b/scripts/euxfel-train-analysis @@ -0,0 +1,86 @@ +#!/usr/bin/env python + +# +# Analysis train and pulse IDs +# +# Copyright (c) 2017 Deutsches Elektronen-Synchrotron DESY, +# a research centre of the Helmholtz Association. +# +# Author: +# 2017 Thomas White <taw@physics.org> +# + +import sys +from collections import deque + +f = open(sys.argv[1], 'r') +n_crystals = 0 +n_frames = 0 +crystal = deque() +in_crystal = 0 +trainID = 0 +pulseID = 0 +trains = dict() +while True: + + fline = f.readline() + if not fline: + break + + fline = fline.rstrip("\r\n") + + if fline.find("Image filename") != -1: + filename = fline.split(": ")[1] + + if fline.find("Event") != -1: + event = fline.split(": ")[1] + + if fline.find("End chunk") != -1: + if n_crystals > 0: + if trainID not in trains: + trains[trainID] = [] + frame = dict() + frame['filename'] = filename+" "+event + frame['ncrystals'] = n_crystals + frame['crystals'] = deque(crystal) + frame['pulseID'] = pulseID + trains[trainID].append(frame) + + if fline.find("Begin chunk") != -1: + n_crystals = 0 + filename = "" + event = "" + pulseID = "" + trainID = "" + crystal.clear() + + if fline.find("Begin crystal") != -1: + n_crystals += 1 + + if fline.find("Cell param") != -1: + in_crystal = 1 + + if fline.find("hdf5/instrument/trainID") != -1: + trainID = fline.split(" = ")[1] + + if fline.find("hdf5/instrument/pulseID") != -1: + pulseID = fline.split(" = ")[1] + + if in_crystal: + crystal.append(fline) + + if fline.find("unique_axis") != -1: + in_crystal = 0 + +for trainID in trains: + print("Train: "+trainID) + print("There are "+str(len(trains[trainID]))+" indexed frames in this train") + if len(trains[trainID]) > 1: + for frame in trains[trainID]: + print("\nFrame: "+frame['filename']) + print("Pulse: "+frame['pulseID']) + print("There are "+str(frame['ncrystals'])+" crystals in this frame") + for line in frame['crystals']: + print(line) + print("\n") + print("\n") diff --git a/scripts/extract-geom b/scripts/extract-geom new file mode 100755 index 00000000..13f7ac2e --- /dev/null +++ b/scripts/extract-geom @@ -0,0 +1,13 @@ +#!/usr/bin/env python + +import sys + +with open(sys.argv[1]) as f: + s = "" + while "Begin geometry file" not in s: + s = f.readline() + s = f.readline() + while "End geometry file" not in s: + print(s), + s = f.readline() + diff --git a/scripts/find-multiples b/scripts/find-multiples new file mode 100755 index 00000000..6196e83c --- /dev/null +++ b/scripts/find-multiples @@ -0,0 +1,58 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Find multiple crystal hits +# +# Copyright © 2017 Deutsches Elektronen-Synchrotron DESY, +# a research centre of the Helmholtz Association. +# +# Author: +# 2017 Thomas White <taw@physics.org> +# + +import sys +from collections import deque + +f = open(sys.argv[1], 'r') +n_crystals = 0 +crystal = deque() +in_crystal = 0 +while True: + + fline = f.readline() + if not fline: + break + + fline = fline.rstrip("\r\n") + + if fline.find("Image filename") != -1: + filename = fline + + if fline.find("Event") != -1: + event = fline + + if fline.find("End chunk") != -1: + if n_crystals > 1: + print("\nFilename: "+filename) + print("Event: "+event) + for line in crystal: + print(line) + + if fline.find("Begin chunk") != -1: + n_crystals = 0 + filename = "" + event = "" + crystal.clear() + + if fline.find("Begin crystal") != -1: + n_crystals += 1 + + if fline.find("Cell param") != -1: + in_crystal = 1 + + if in_crystal: + crystal.append(fline) + + if fline.find("unique_axis") != -1: + in_crystal = 0 + diff --git a/scripts/sum-peaks b/scripts/sum-peaks index 5b528661..c239e113 100755 --- a/scripts/sum-peaks +++ b/scripts/sum-peaks @@ -15,37 +15,41 @@ import numpy as np import h5py from os.path import basename, splitext -with open(sys.argv[1], 'r') as stream: - reading_geometry = False - reading_chunks = False - reading_peaks = False - max_fs = -100500 - max_ss = -100500 - for line in stream: - if reading_chunks: - if line.startswith('End of peak list'): - reading_peaks = False - elif line.startswith(' fs/px ss/px (1/d)/nm^-1 Intensity Panel'): - reading_peaks = True - elif reading_peaks: - fs, ss, dump, intensity = [float(i) for i in line.split()[:4]] - powder[int(ss), int(fs)] += intensity - elif line.startswith('----- End geometry file -----'): - reading_geometry = False - powder = np.zeros((max_ss + 1, max_fs + 1)) - elif reading_geometry: - try: - par, val = line.split('=') - if par.split('/')[-1].strip() == 'max_fs' and int(val) > max_fs: - max_fs = int(val) - elif par.split('/')[-1].strip() == 'max_ss' and int(val) > max_ss: - max_ss = int(val) - except ValueError: - pass - elif line.startswith('----- Begin geometry file -----'): - reading_geometry = True - elif line.startswith('----- Begin chunk -----'): - reading_chunks = True +if sys.argv[1] == '-': + stream = sys.stdin +else: + stream = open(sys.argv[1], 'r') + +reading_geometry = False +reading_chunks = False +reading_peaks = False +max_fs = -100500 +max_ss = -100500 +for line in stream: + if reading_chunks: + if line.startswith('End of peak list'): + reading_peaks = False + elif line.startswith(' fs/px ss/px (1/d)/nm^-1 Intensity Panel'): + reading_peaks = True + elif reading_peaks: + fs, ss, dump, intensity = [float(i) for i in line.split()[:4]] + powder[int(ss), int(fs)] += intensity + elif line.startswith('----- End geometry file -----'): + reading_geometry = False + powder = np.zeros((max_ss + 1, max_fs + 1)) + elif reading_geometry: + try: + par, val = line.split('=') + if par.split('/')[-1].strip() == 'max_fs' and int(val) > max_fs: + max_fs = int(val) + elif par.split('/')[-1].strip() == 'max_ss' and int(val) > max_ss: + max_ss = int(val) + except ValueError: + pass + elif line.startswith('----- Begin geometry file -----'): + reading_geometry = True + elif line.startswith('----- Begin chunk -----'): + reading_chunks = True f = h5py.File(splitext(basename(sys.argv[1]))[0]+'-powder.h5', 'w') f.create_dataset('/data/data', data=powder) diff --git a/scripts/transfer-geom b/scripts/transfer-geom new file mode 100755 index 00000000..6109f5ed --- /dev/null +++ b/scripts/transfer-geom @@ -0,0 +1,91 @@ +#!/usr/bin/python3 + +import sys + +corner_x = dict() +corner_y = dict() +fs = dict() +ss = dict() +coffset = dict() + +# Go through "ingeom" and find all corner_{x,y}/fs/ss/coffset lines +with open(sys.argv[1]) as ingeom: + while True: + line = ingeom.readline() + if not line: + break + line = line.rstrip("\r\n") + slash = line.find("/") + if slash != -1: + panel = line[0:slash] + else: + panel = None + if line.find("/corner_x") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + corner_x[panel] = line.split(" = ")[1] + if line.find("/corner_y") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + corner_y[panel] = line.split(" = ")[1] + if line.find("/fs") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + fs[panel] = line.split(" = ")[1] + if line.find("/ss") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + ss[panel] = line.split(" = ")[1] + if line.find("/coffset") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + coffset[panel] = line.split(" = ")[1] + +# Go through "template", and replace all corner_{x,y}/fs/ss lines from memory +out = open(sys.argv[3], "w") +with open(sys.argv[2]) as template: + while True: + line = template.readline() + if not line: + break + line = line.rstrip("\r\n") + slash = line.find("/") + if slash != -1: + panel = line[0:slash] + else: + panel = None + if line.find("/corner_x") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + out.write("%s/corner_x = %s\n" % (panel, corner_x[panel])) + elif line.find("/corner_y") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + out.write("%s/corner_y = %s\n" % (panel, corner_y[panel])) + elif line.find("/fs") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + out.write("%s/fs = %s\n" % (panel, fs[panel])) + elif line.find("/ss") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + out.write("%s/ss = %s\n" % (panel, ss[panel])) + elif line.find("/coffset") != -1: + if not panel: + print("No panel for line: %s\n" % line) + sys.exit(1) + out.write("%s/coffset = %s\n" % (panel, coffset[panel])) + else: + out.write("%s\n" % line) + +out.close() + diff --git a/scripts/turbo-index-lsf b/scripts/turbo-index-lsf index 00e4ec15..1640ec3d 100644 --- a/scripts/turbo-index-lsf +++ b/scripts/turbo-index-lsf @@ -1,26 +1,76 @@ #!/bin/sh -RUN=$1 -NOSAMPLE=`echo $RUN | sed -e 's/\-.*$//'` +# Split a large indexing job into many small tasks and submit using LSF -GEOM=my.geom # Name of your geometry file +# ./turbo-index my-files.lst label my.geom /location/for/streams -find /path/to/CXI/files/$RUN -name '*.cxi' > files-${RUN}.lst # Set location of files -list_events -i files-${RUN}.lst -g $GEOM -o events-${RUN}.lst +# Copyright © 2016-2017 Deutsches Elektronen-Synchrotron DESY, +# a research centre of the Helmholtz Association. +# +# Authors: +# 2016 Steve Aplin <steve.aplin@desy.de> +# 2016-2017 Thomas White <taw@physics.org> + +SPLIT=1000 # Size of job chunks + +INPUT=$1 +RUN=$2 +GEOM=$3 +STREAMDIR=$4 + +# Set up environment here if necessary +#source /path/to/crystfel/setup.sh + +# Generate event list from file above +list_events -i $INPUT -g $GEOM -o events-${RUN}.lst +if [ $? != 0 ]; then + echo "list_events failed" + exit 1 +fi +# If you are using single-event files instead of multi-event ("CXI") ones, +# comment out the above lines and uncomment the following one: +#cp $INPUT events-${RUN}.lst + +# Count total number of events wc -l events-${RUN}.lst -rm -f split-events-${RUN}.lst files-${RUN}.lst -split -d -l 500 events-${RUN}.lst split-events-${RUN}.lst + +# Split the events up, will create files with $SPLIT lines +split -a 3 -d -l $SPLIT events-${RUN}.lst split-events-${RUN}.lst + +# Clean up rm -f events-${RUN}.lst +# Loop over the event list files, and submit a batch job for each of them for FILE in split-events-${RUN}.lst*; do - STREAM=`echo $FILE | sed -e "s/split-events-${RUN}.lst/${RUN}.stream/"` - NAME=`echo $FILE | sed -e "s/split-events-${RUN}.lst/${NOSAMPLE}-/"` - echo "$NAME: $FILE ---> $STREAM" + # Stream file is the output of crystfel + STREAM=`echo $FILE | sed -e "s/split-events-${RUN}.lst/${RUN}.stream/"` + + # Job name + NAME=`echo $FILE | sed -e "s/split-events-${RUN}.lst/${RUN}-/"` + + # Job number + NUMBER=${NAME##$RUN-} + POS=`expr $NUMBER \* $SPLIT + 1` + + echo "$NAME (serial start $POS): $FILE ---> $STREAM" + + SLURMFILE="${NAME}.sh" + + echo "#!/bin/sh" > $SLURMFILE + echo >> $SLURMFILE + + # Set up environment here (again) if necessary + echo "#source /path/to/crystfel/setup.sh" >> $SLURMFILE + echo >> $SLURMFILE + + command="indexamajig -i $FILE -o $STREAMDIR/$STREAM --serial-start=$POS" + command="$command -j \`nproc\` -g $GEOM" + #command="$command --peaks=zaef" # Indexing parameters here + + echo $command >> $SLURMFILE - # Set indexing parameters here - bsub -q psanaq -o $NAME.log -J $NAME -n 12 -R "span[hosts=1]" \ - indexamajig \ - -i $FILE -o $STREAM -j 32 -g $GEOM --peaks=cxi + # Set your queue after "-q" below + bsub -q myqueue -o $NAME.out -e $NAME.err -J $NAME -x sh $SLURMFILE done diff --git a/src/cell_explorer.c b/src/cell_explorer.c index 94e7e985..69dd7178 100644 --- a/src/cell_explorer.c +++ b/src/cell_explorer.c @@ -37,7 +37,7 @@ #include <string.h> #include <gtk/gtk.h> #include <math.h> -#include <gdk/gdkkeysyms.h> +#include <gdk/gdkkeysyms-compat.h> #include <gsl/gsl_multifit_nlin.h> #include "stream.h" @@ -189,6 +189,12 @@ static void set_col(cairo_t *cr, CellWindow *w, int cat) static gboolean destroy_sig(GtkWidget *da, CellWindow *w) { + multihistogram_free(w->hist_a->h); + multihistogram_free(w->hist_b->h); + multihistogram_free(w->hist_c->h); + multihistogram_free(w->hist_al->h); + multihistogram_free(w->hist_be->h); + multihistogram_free(w->hist_ga->h); gtk_main_quit(); return FALSE; } @@ -315,21 +321,20 @@ static void draw_label(cairo_t *cr, HistoBox *b, int width, int height) } -static gboolean draw_sig(GtkWidget *da, GdkEventExpose *event, HistoBox *b) +static gboolean draw_sig(GtkWidget *da, cairo_t *cr, HistoBox *b) { int width, height; int i, max; double h_height; - cairo_t *cr; double gstep; int *data_p, *data_a, *data_b, *data_c, *data_i, *data_f; int *data_r, *data_h, *data_excl; int start, stop; + GtkAllocation allocation; - cr = gdk_cairo_create(da->window); - - width = da->allocation.width; - height = da->allocation.height; + gtk_widget_get_allocation(da, &allocation); + width = allocation.width; + height = allocation.height; b->width = width; /* Store for later use when dragging */ /* Overall background */ @@ -478,8 +483,16 @@ static gboolean draw_sig(GtkWidget *da, GdkEventExpose *event, HistoBox *b) draw_axis(cr, b, width, height); draw_label(cr, b, width, height); - cairo_destroy(cr); + return FALSE; +} + +static gboolean expose_sig(GtkWidget *da, GdkEventExpose *event, HistoBox *b) +{ + cairo_t *cr; + cr = gdk_cairo_create(gtk_widget_get_window(da)); + draw_sig(da, cr, b); + cairo_destroy(cr); return FALSE; } @@ -510,11 +523,13 @@ static gint keyclick_sig(GtkWidget *widget, GdkEventButton *event, CellWindow *w) { int width, cat; + GtkAllocation alloc; /* Ignore extra events for double click */ if ( event->type != GDK_BUTTON_PRESS ) return FALSE; - width = widget->allocation.width; + gtk_widget_get_allocation(widget, &alloc); + width = alloc.width; cat = 8*event->x / width; @@ -534,16 +549,15 @@ static gint keyclick_sig(GtkWidget *widget, GdkEventButton *event, -static gboolean keydraw_sig(GtkWidget *da, GdkEventExpose *event, CellWindow *w) +static gboolean keydraw_sig(GtkWidget *da, cairo_t *cr, CellWindow *w) { int width, height; - cairo_t *cr; double x; + GtkAllocation allocation; - cr = gdk_cairo_create(da->window); - - width = da->allocation.width; - height = da->allocation.height; + gtk_widget_get_allocation(da, &allocation); + width = allocation.width; + height = allocation.height; /* Overall background */ cairo_rectangle(cr, 0.0, 0.0, width, height); @@ -602,6 +616,17 @@ static gboolean keydraw_sig(GtkWidget *da, GdkEventExpose *event, CellWindow *w) } +static gboolean keyexpose_sig(GtkWidget *da, GdkEventExpose *event, + CellWindow *w) +{ + cairo_t *cr; + cr = gdk_cairo_create(gtk_widget_get_window(da)); + keydraw_sig(da, cr, w); + cairo_destroy(cr); + return FALSE; +} + + static int check_exclude(HistoBox *h, double v) { double min, max; @@ -647,7 +672,7 @@ static void scan_cells(CellWindow *w) for ( i=0; i<n_cells; i++ ) { double a, b, c, al, be, ga; - int cat, j; + unsigned int cat, j; int ignore = 0; for ( j=0; j<w->n_unique_indms; j++ ) { @@ -687,7 +712,7 @@ static void scan_cells(CellWindow *w) check_exclude(w->hist_be, be) || check_exclude(w->hist_ga, ga) ) { - cat = 1<<CAT_EXCLUDE; + cat = (unsigned)1<<CAT_EXCLUDE; n_excl++; } @@ -1162,8 +1187,8 @@ static int write_cell_to_file(UnitCell *cell, const char *filename) static gint savecell_sig(GtkWidget *widget, CellWindow *w) { GtkWidget *d; - gchar *output_filename; UnitCell *cell; + gint r; cell = get_cell(w); if ( cell == NULL ) return FALSE; @@ -1176,19 +1201,19 @@ static gint savecell_sig(GtkWidget *widget, CellWindow *w) NULL); gtk_file_chooser_set_do_overwrite_confirmation(GTK_FILE_CHOOSER(d), TRUE); + r = gtk_dialog_run(GTK_DIALOG(d)); + if ( r == GTK_RESPONSE_ACCEPT ) { - if ( gtk_dialog_run(GTK_DIALOG(d)) == GTK_RESPONSE_CANCEL ) { - gtk_widget_destroy(d); - return FALSE; - } - output_filename = gtk_file_chooser_get_filename(GTK_FILE_CHOOSER(d)); + gchar *output_filename; - if ( write_cell_to_file(cell, output_filename) ) { - error_box(w, "Failed to save unit cell"); - } + output_filename = gtk_file_chooser_get_filename(GTK_FILE_CHOOSER(d)); + if ( write_cell_to_file(cell, output_filename) ) { + error_box(w, "Failed to save unit cell"); + } + g_free(output_filename); + } gtk_widget_destroy(d); - g_free(output_filename); return FALSE; } @@ -1205,7 +1230,7 @@ static gint about_sig(GtkWidget *widget, CellWindow *w) window = gtk_about_dialog_new(); gtk_window_set_transient_for(GTK_WINDOW(window), GTK_WINDOW(w->window)); - gtk_about_dialog_set_name(GTK_ABOUT_DIALOG(window), + gtk_about_dialog_set_program_name(GTK_ABOUT_DIALOG(window), "Unit Cell Explorer"); gtk_about_dialog_set_version(GTK_ABOUT_DIALOG(window), "0.0.1"); gtk_about_dialog_set_copyright(GTK_ABOUT_DIALOG(window), @@ -1336,8 +1361,11 @@ static gint motion_sig(GtkWidget *da, GdkEventMotion *event, HistoBox *h) gtk_widget_queue_draw(h->da); - if ( event->is_hint ) gdk_window_get_pointer(da->window, - NULL, NULL, NULL); + if ( event->is_hint ) { + gdk_window_get_pointer(gtk_widget_get_window(da), + NULL, NULL, NULL); + } + return TRUE; } @@ -1419,8 +1447,14 @@ static HistoBox *histobox_new(CellWindow *w, const char *units, const char *n) | GDK_SCROLL_MASK | GDK_KEY_PRESS_MASK); - g_signal_connect(G_OBJECT(h->da), "expose_event", G_CALLBACK(draw_sig), - h); + if ( g_signal_lookup("draw", GTK_TYPE_DRAWING_AREA) ) { + g_signal_connect(G_OBJECT(h->da), "draw", + G_CALLBACK(draw_sig), h); + } else { + g_signal_connect(G_OBJECT(h->da), "expose-event", + G_CALLBACK(expose_sig), h); + } + g_signal_connect(G_OBJECT(h->da), "button-press-event", G_CALLBACK(press_sig), h); g_signal_connect(G_OBJECT(h->da), "button-release-event", @@ -1466,8 +1500,15 @@ static void indexing_method_list(CellWindow *w, GtkWidget *vbox) key = gtk_drawing_area_new(); gtk_box_pack_end(GTK_BOX(w->indmlist), key, FALSE, FALSE, 5.0); gtk_widget_add_events(GTK_WIDGET(key), GDK_BUTTON_PRESS_MASK); - g_signal_connect(G_OBJECT(key), "expose_event", G_CALLBACK(keydraw_sig), - w); + + if ( g_signal_lookup("draw", GTK_TYPE_DRAWING_AREA) ) { + g_signal_connect(G_OBJECT(key), "draw", + G_CALLBACK(keydraw_sig), w); + } else { + g_signal_connect(G_OBJECT(key), "expose-event", + G_CALLBACK(keyexpose_sig), w); + } + g_signal_connect(G_OBJECT(key), "configure-event", G_CALLBACK(keyconf_sig), w); g_signal_connect(G_OBJECT(key), "button-press-event", diff --git a/src/compare_hkl.c b/src/compare_hkl.c index feb79879..5812cad2 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -1438,6 +1438,8 @@ int main(int argc, char *argv[]) list1 = asymmetric_indices(list1_raw, sym); list2 = asymmetric_indices(list2_raw, sym); + reflist_free(list1_raw); + reflist_free(list2_raw); /* Select reflections to be used */ ncom = 0; @@ -1521,12 +1523,58 @@ int main(int argc, char *argv[]) } } - if ( (fom == FOM_CCANO) || (fom == FOM_CRDANO) - || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) ) + refl1_acc = add_refl(list1_acc, h, k, l); + copy_data(refl1_acc, refl1); + set_intensity(refl1_acc, val1); + + refl2_acc = add_refl(list2_acc, h, k, l); + copy_data(refl2_acc, refl2); + set_intensity(refl2_acc, val2); + + if ( val1 > max_I ) max_I = val1; + if ( val1 < min_I ) min_I = val1; + + ncom++; + + } + + reflist_free(list1); + reflist_free(list2); + + /* For anomalous figures of merit, we additionally require that we have + * all the Bijvoet pairs after the above rejection tests */ + if ( (fom == FOM_CCANO) || (fom == FOM_CRDANO) + || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) ) + { + list1 = list1_acc; + list2 = list2_acc; + list1_acc = reflist_new(); + list2_acc = reflist_new(); + + min_I = +INFINITY; + max_I = -INFINITY; + ncom = 0; + + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { Reflection *refl1_bij = NULL; Reflection *refl2_bij = NULL; + signed int h, k, l; signed int hb, kb, lb; + Reflection *refl1_acc; + Reflection *refl2_acc; + Reflection *refl2; + double val1, val2; + + get_indices(refl1, &h, &k, &l); + + refl2 = find_refl(list2, h, k, l); + assert(refl2 != NULL); + + val1 = get_intensity(refl1); + val2 = get_intensity(refl2); if ( is_centric(h, k, l, sym) ) { ncen++; @@ -1549,21 +1597,20 @@ int main(int argc, char *argv[]) nbij++; continue; } - } - refl1_acc = add_refl(list1_acc, h, k, l); - copy_data(refl1_acc, refl1); - set_intensity(refl1_acc, val1); + refl1_acc = add_refl(list1_acc, h, k, l); + copy_data(refl1_acc, refl1); + set_intensity(refl1_acc, val1); - refl2_acc = add_refl(list2_acc, h, k, l); - copy_data(refl2_acc, refl2); - set_intensity(refl2_acc, val2); + refl2_acc = add_refl(list2_acc, h, k, l); + copy_data(refl2_acc, refl2); + set_intensity(refl2_acc, val2); - if ( val1 > max_I ) max_I = val1; - if ( val1 < min_I ) min_I = val1; - - ncom++; + if ( val1 > max_I ) max_I = val1; + if ( val1 < min_I ) min_I = val1; + ncom++; + } } gsl_set_error_handler_off(); @@ -1611,11 +1658,6 @@ int main(int argc, char *argv[]) " (%.2f to %.2f Angstroms).\n", rmin/1e9, rmax/1e9, 1e10/rmin, 1e10/rmax); - reflist_free(list1_raw); - reflist_free(list2_raw); - reflist_free(list1); - reflist_free(list2); - if ( rmin_fix >= 0.0 ) { rmin = rmin_fix; } diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c index b0ed52ea..79b5443d 100644 --- a/src/diffraction-gpu.c +++ b/src/diffraction-gpu.c @@ -342,16 +342,16 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, for ( i=0; i<image->nsamples; i++ ) { printf("%.3f eV, weight = %.5f\n", - ph_lambda_to_eV(1.0/image->spectrum[i].k), - image->spectrum[i].weight); + ph_lambda_to_eV(1.0/image->spectrum0[i].k), + image->spectrum0[i].weight); - err = do_panels(gctx, image, image->spectrum[i].k, - image->spectrum[i].weight, + err = do_panels(gctx, image, image->spectrum0[i].k, + image->spectrum0[i].weight, &n_inf, &n_neg, &n_nan); if ( err ) return 1; - tot += image->spectrum[i].weight; + tot += image->spectrum0[i].weight; } printf("total weight = %f\n", tot); diff --git a/src/diffraction.c b/src/diffraction.c index 4cceba0e..4dc21c59 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -608,6 +608,40 @@ struct sample *generate_tophat(struct image *image) } +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; i<image->nsamples; 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; i<nsamples; i++ ) { + spectrum[i].weight /= w_sum; + } + + image->spectrum_size = nsamples; + + return spectrum; +} + + struct sample *generate_SASE(struct image *image, gsl_rng *rng) { struct sample *spectrum; @@ -738,13 +772,13 @@ void get_diffraction(struct image *image, int na, int nb, int nc, for ( i=0; i<image->nsamples; i++ ) { printf("%.1f eV, weight = %.5f\n", - ph_lambda_to_eV(1.0/image->spectrum[i].k), - image->spectrum[i].weight); + ph_lambda_to_eV(1.0/image->spectrum0[i].k), + image->spectrum0[i].weight); diffraction_at_k(image, intensities, phases, - flags, cell, m, sym, image->spectrum[i].k, + flags, cell, m, sym, image->spectrum0[i].k, ax, ay, az, bx, by, bz, cx, cy, cz, - lut_a, lut_b, lut_c, image->spectrum[i].weight); + lut_a, lut_b, lut_c, image->spectrum0[i].weight); } diff --git a/src/diffraction.h b/src/diffraction.h index ca4e3bde..0f95f83c 100644 --- a/src/diffraction.h +++ b/src/diffraction.h @@ -60,4 +60,7 @@ extern struct sample *generate_SASE(struct image *image, gsl_rng *rng); extern struct sample *generate_twocolour(struct image *image); +extern struct sample *generate_spectrum_fromfile(struct image *image, + char *spectrum_fn); + #endif /* DIFFRACTION_H */ diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c index 7fa49bb2..7545e4ab 100644 --- a/src/dw-hdfsee.c +++ b/src/dw-hdfsee.c @@ -41,7 +41,7 @@ #include <string.h> #include <cairo.h> #include <gdk-pixbuf/gdk-pixbuf.h> -#include <gdk/gdkkeysyms.h> +#include <gdk/gdkkeysyms-compat.h> #include <assert.h> #include "dw-hdfsee.h" @@ -531,7 +531,8 @@ static void redraw_window(DisplayWindow *dw) width, dw->height); draw_stuff(dw->surf, dw); - gdk_window_invalidate_rect(dw->drawingarea->window, NULL, FALSE); + gdk_window_invalidate_rect(gtk_widget_get_window(dw->drawingarea), + NULL, FALSE); } @@ -612,20 +613,21 @@ static void displaywindow_update(DisplayWindow *dw) } +static gboolean displaywindow_draw(GtkWidget *da, cairo_t *cr, DisplayWindow *dw) +{ + cairo_set_source_surface(cr, dw->surf, 0.0, 0.0); + cairo_paint(cr); + return FALSE; +} + + static gboolean displaywindow_expose(GtkWidget *da, GdkEventExpose *event, DisplayWindow *dw) { cairo_t *cr; - - cr = gdk_cairo_create(da->window); - - cairo_set_source_surface(cr, dw->surf, 0.0, 0.0); - cairo_rectangle(cr, event->area.x, event->area.y, - event->area.width, event->area.height); - cairo_fill(cr); - + cr = gdk_cairo_create(gtk_widget_get_window(da)); + displaywindow_draw(da, cr, dw); cairo_destroy(cr); - return FALSE; } @@ -721,6 +723,7 @@ static gint displaywindow_set_binning(GtkWidget *widget, DisplayWindow *dw) GtkWidget *label; char tmp[64]; double minx, maxx, miny, maxy; + GtkWidget *dvbox; if ( dw->binning_dialog != NULL ) { return 0; @@ -741,10 +744,10 @@ static gint displaywindow_set_binning(GtkWidget *widget, DisplayWindow *dw) GTK_STOCK_OK, GTK_RESPONSE_OK, NULL); + dvbox = gtk_dialog_get_content_area((GTK_DIALOG(bd->window))); vbox = gtk_vbox_new(FALSE, 0); hbox = gtk_hbox_new(TRUE, 0); - gtk_box_pack_start(GTK_BOX(GTK_DIALOG(bd->window)->vbox), - GTK_WIDGET(hbox), FALSE, FALSE, 7); + gtk_box_pack_start(GTK_BOX(dvbox), GTK_WIDGET(hbox), FALSE, FALSE, 7); gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5); table = gtk_table_new(3, 2, FALSE); @@ -850,6 +853,7 @@ static gint displaywindow_set_boostint(GtkWidget *widget, DisplayWindow *dw) GtkWidget *table; GtkWidget *label; char tmp[64]; + GtkWidget *dvbox; if ( dw->boostint_dialog != NULL ) { return 0; @@ -869,10 +873,10 @@ static gint displaywindow_set_boostint(GtkWidget *widget, DisplayWindow *dw) GTK_STOCK_CANCEL, GTK_RESPONSE_CLOSE, GTK_STOCK_OK, GTK_RESPONSE_OK, NULL); + dvbox = gtk_dialog_get_content_area(GTK_DIALOG(bd->window)); vbox = gtk_vbox_new(FALSE, 0); hbox = gtk_hbox_new(TRUE, 0); - gtk_box_pack_start(GTK_BOX(GTK_DIALOG(bd->window)->vbox), - GTK_WIDGET(hbox), FALSE, FALSE, 7); + gtk_box_pack_start(GTK_BOX(dvbox), GTK_WIDGET(hbox), FALSE, FALSE, 7); gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5); table = gtk_table_new(3, 2, FALSE); @@ -1043,6 +1047,7 @@ static gint displaywindow_set_newevent(GtkWidget *widget, DisplayWindow *dw) GtkWidget *hbox; GtkWidget *table; GtkWidget *label; + GtkWidget *dvbox; char tmp[1024]; if ( dw->event_dialog != NULL ) { @@ -1063,10 +1068,10 @@ static gint displaywindow_set_newevent(GtkWidget *widget, DisplayWindow *dw) GTK_STOCK_CANCEL, GTK_RESPONSE_CLOSE, GTK_STOCK_OK, GTK_RESPONSE_OK, NULL); + dvbox = gtk_dialog_get_content_area(GTK_DIALOG(ed->window)); vbox = gtk_vbox_new(FALSE, 0); hbox = gtk_hbox_new(TRUE, 0); - gtk_box_pack_start(GTK_BOX(GTK_DIALOG(ed->window)->vbox), - GTK_WIDGET(hbox), FALSE, FALSE, 7); + gtk_box_pack_start(GTK_BOX(dvbox), GTK_WIDGET(hbox), FALSE, FALSE, 7); gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5); table = gtk_table_new(3, 2, FALSE); @@ -1159,6 +1164,7 @@ static gint displaywindow_set_ringradius(GtkWidget *widget, DisplayWindow *dw) GtkWidget *hbox; GtkWidget *table; GtkWidget *label; + GtkWidget *dvbox; char tmp[64]; if ( dw->ringradius_dialog != NULL ) { @@ -1179,10 +1185,10 @@ static gint displaywindow_set_ringradius(GtkWidget *widget, DisplayWindow *dw) GTK_STOCK_CANCEL, GTK_RESPONSE_CLOSE, GTK_STOCK_OK, GTK_RESPONSE_OK, NULL); + dvbox = gtk_dialog_get_content_area(GTK_DIALOG(rd->window)); vbox = gtk_vbox_new(FALSE, 0); hbox = gtk_hbox_new(TRUE, 0); - gtk_box_pack_start(GTK_BOX(GTK_DIALOG(rd->window)->vbox), - GTK_WIDGET(hbox), FALSE, FALSE, 7); + gtk_box_pack_start(GTK_BOX(dvbox), GTK_WIDGET(hbox), FALSE, FALSE, 7); gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5); table = gtk_table_new(3, 2, FALSE); @@ -1414,7 +1420,7 @@ static gint displaywindow_about(GtkWidget *widget, DisplayWindow *dw) gtk_window_set_transient_for(GTK_WINDOW(window), GTK_WINDOW(dw->window)); - gtk_about_dialog_set_name(GTK_ABOUT_DIALOG(window), "hdfsee"); + gtk_about_dialog_set_program_name(GTK_ABOUT_DIALOG(window), "hdfsee"); gtk_about_dialog_set_version(GTK_ABOUT_DIALOG(window), PACKAGE_VERSION); gtk_about_dialog_set_copyright(GTK_ABOUT_DIALOG(window), "© 2012-2015 Deutsches Elektronen-Synchrotron DESY," @@ -1542,7 +1548,7 @@ static gint displaywindow_set_calibmode(GtkWidget *d, DisplayWindow *dw) w = gtk_ui_manager_get_widget(dw->ui, "/ui/displaywindow/tools/calibmode"); if ( dw->simple ) { - gtk_check_menu_item_set_state(GTK_CHECK_MENU_ITEM(w), 0); + gtk_check_menu_item_set_active(GTK_CHECK_MENU_ITEM(w), 0); return 0; } @@ -1550,8 +1556,7 @@ static gint displaywindow_set_calibmode(GtkWidget *d, DisplayWindow *dw) dw->rg_coll = dw->image->det->rigid_group_collections[0]; if (dw->rg_coll == NULL) { ERROR("Cannot find asuitable rigid group collection.\n"); - gtk_check_menu_item_set_state(GTK_CHECK_MENU_ITEM(w), - 0); + gtk_check_menu_item_set_active(GTK_CHECK_MENU_ITEM(w), 0); return 0; } } else { @@ -1560,8 +1565,7 @@ static gint displaywindow_set_calibmode(GtkWidget *d, DisplayWindow *dw) if (dw->rg_coll == NULL) { ERROR("Cannot find rigid group collection: %s\n", dw->rg_coll_name); - gtk_check_menu_item_set_state(GTK_CHECK_MENU_ITEM(w), - 0); + gtk_check_menu_item_set_active(GTK_CHECK_MENU_ITEM(w), 0); return 0; } } @@ -1807,18 +1811,18 @@ static gint displaywindow_save(GtkWidget *widget, DisplayWindow *dw) hbox = gtk_hbox_new(FALSE, 0); gtk_file_chooser_set_extra_widget(GTK_FILE_CHOOSER(d), hbox); - cb = gtk_combo_box_new_text(); + cb = gtk_combo_box_text_new(); gtk_box_pack_end(GTK_BOX(hbox), GTK_WIDGET(cb), TRUE, TRUE, 5); l = gtk_label_new("Save as type:"); gtk_box_pack_end(GTK_BOX(hbox), GTK_WIDGET(l), FALSE, FALSE, 5); - gtk_combo_box_append_text(GTK_COMBO_BOX(cb), + gtk_combo_box_text_append_text(GTK_COMBO_BOX_TEXT(cb), "PNG (colour)"); - gtk_combo_box_append_text(GTK_COMBO_BOX(cb), + gtk_combo_box_text_append_text(GTK_COMBO_BOX_TEXT(cb), "TIFF (floating point)"); - gtk_combo_box_append_text(GTK_COMBO_BOX(cb), + gtk_combo_box_text_append_text(GTK_COMBO_BOX_TEXT(cb), "TIFF (16 bit signed integer)"); - gtk_combo_box_append_text(GTK_COMBO_BOX(cb), + gtk_combo_box_text_append_text(GTK_COMBO_BOX_TEXT(cb), "ADSC (16 bit unsigned integer)"); gtk_combo_box_set_active(GTK_COMBO_BOX(cb), 0); @@ -1885,6 +1889,7 @@ static gint displaywindow_show_numbers(GtkWidget *widget, DisplayWindow *dw) GtkWidget *hbox2; GtkWidget *table; GtkWidget *label; + GtkWidget *dvbox; unsigned int x, y; if ( dw->numbers_window != NULL ) { @@ -1905,10 +1910,10 @@ static gint displaywindow_show_numbers(GtkWidget *widget, DisplayWindow *dw) GTK_STOCK_CLOSE, GTK_RESPONSE_CLOSE, NULL); + dvbox = gtk_dialog_get_content_area(GTK_DIALOG(nw->window)); vbox = gtk_vbox_new(FALSE, 0); hbox = gtk_hbox_new(TRUE, 0); - gtk_box_pack_start(GTK_BOX(GTK_DIALOG(nw->window)->vbox), - GTK_WIDGET(hbox), FALSE, FALSE, 7); + gtk_box_pack_start(GTK_BOX(dvbox), GTK_WIDGET(hbox), FALSE, FALSE, 7); gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5); table = gtk_table_new(17, 17, FALSE); @@ -2624,7 +2629,7 @@ static gint displaywindow_release(GtkWidget *widget, GdkEventButton *event, { if ( (event->type == GDK_BUTTON_RELEASE) && (event->button == 1) ) { - g_signal_handler_disconnect(GTK_OBJECT(dw->drawingarea), + g_signal_handler_disconnect(G_OBJECT(dw->drawingarea), dw->motion_callback); dw->motion_callback = 0; @@ -2699,7 +2704,7 @@ static gint displaywindow_press(GtkWidget *widget, GdkEventButton *event, if ( (event->type == GDK_BUTTON_PRESS) && (event->button == 1) ) { dw->motion_callback = g_signal_connect( - GTK_OBJECT(dw->drawingarea), + G_OBJECT(dw->drawingarea), "motion-notify-event", G_CALLBACK(displaywindow_motion), dw); @@ -3132,8 +3137,13 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, dw->drawingarea); gtk_box_pack_start(GTK_BOX(vbox), dw->scrollarea, TRUE, TRUE, 0); - g_signal_connect(GTK_OBJECT(dw->drawingarea), "expose-event", - G_CALLBACK(displaywindow_expose), dw); + if ( g_signal_lookup("draw", GTK_TYPE_DRAWING_AREA) ) { + g_signal_connect(G_OBJECT(dw->drawingarea), "draw", + G_CALLBACK(displaywindow_draw), dw); + } else { + g_signal_connect(G_OBJECT(dw->drawingarea), "expose-event", + G_CALLBACK(displaywindow_expose), dw); + } gtk_window_set_resizable(GTK_WINDOW(dw->window), TRUE); gtk_widget_show_all(dw->window); @@ -3162,11 +3172,11 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, gtk_widget_grab_focus(dw->drawingarea); - g_signal_connect(GTK_OBJECT(dw->drawingarea), "button-press-event", + g_signal_connect(G_OBJECT(dw->drawingarea), "button-press-event", G_CALLBACK(displaywindow_press), dw); - g_signal_connect(GTK_OBJECT(dw->drawingarea), "button-release-event", + g_signal_connect(G_OBJECT(dw->drawingarea), "button-release-event", G_CALLBACK(displaywindow_release), dw); - g_signal_connect(GTK_OBJECT(dw->drawingarea), "key-press-event", + g_signal_connect(G_OBJECT(dw->drawingarea), "key-press-event", G_CALLBACK(displaywindow_keypress), dw); if ( dw->simple ) { diff --git a/src/geoptimiser.c b/src/geoptimiser.c index 9b03d1b9..f66ac709 100644 --- a/src/geoptimiser.c +++ b/src/geoptimiser.c @@ -674,7 +674,8 @@ static int add_distance_to_list(struct gpanel *gp, } -static int count_pixels_with_min_peaks(struct gpanel *gp, int min_num_peaks) +static int count_pixels_with_min_peaks(struct gpanel *gp, int min_num_peaks, + int max_num_peaks) { int pixel_count = 0; int ifs, iss; @@ -683,7 +684,8 @@ static int count_pixels_with_min_peaks(struct gpanel *gp, int min_num_peaks) for ( ifs=0; ifs<gp->p->w; ifs++ ) { int idx = ifs+gp->p->w*iss; - if ( gp->num_pix_displ[idx] >= min_num_peaks ) { + if ( gp->num_pix_displ[idx] >= min_num_peaks && + gp->num_pix_displ[idx] < max_num_peaks) { pixel_count += 1; } @@ -722,7 +724,7 @@ static void adjust_min_peaks_per_conn(struct rg_collection *connected, gp = &gpanels[panel_number(det, p)]; pix_count = count_pixels_with_min_peaks(gp, - min_num_peaks); + min_num_peaks, gparams->max_num_peaks_per_pix); di_count += pix_count; @@ -829,6 +831,7 @@ static int compute_avg_pix_displ(struct gpanel *gp, int idx, } else { + gp->num_pix_displ[idx] = -1; gp->avg_displ_x[idx] = -10000.0; gp->avg_displ_y[idx] = -10000.0; gp->avg_displ_abs[idx] = -10000.0; @@ -916,7 +919,8 @@ static double compute_error(struct rg_collection *connected, pix_index = ifs+p->w*iss; if ( gp->num_pix_displ[pix_index] - >= conn_data[di].num_peaks_per_pixel ) + >= conn_data[di].num_peaks_per_pixel && + gp->avg_displ_abs[pix_index] > -9999.0 ) { double cer; @@ -1304,7 +1308,8 @@ static int compute_shift_for_empty_panels(struct rg_collection *quadrants, "number of measurements.\n"); for ( di=0; di<connected->n_rigid_groups; di++ ) { - if ( conn_data[di].n_peaks_in_conn >= min_pix ) { + if ( conn_data[di].n_peaks_in_conn >= min_pix && + conn_data[di].sh_x > -9999.0 ) { aver_x[conn_data[di].num_quad] += conn_data[di].sh_x; aver_y[conn_data[di].num_quad] += conn_data[di].sh_y; n[conn_data[di].num_quad]++; @@ -2530,7 +2535,7 @@ int optimize_geometry(struct geoptimiser_params *gparams, if ( gparams->stretch_map ) { -#ifdef HAVE_SAVE_TO_PNG +#ifdef CAN_SAVE_TO_PNG STATUS("Saving stretch map - for out-of-plane rotations.\n"); @@ -2543,12 +2548,12 @@ int optimize_geometry(struct geoptimiser_params *gparams, return 1; } -#else /* HAVE_SAVE_TO_PNG */ +#else /* CAN_SAVE_TO_PNG */ STATUS("ERROR: geoptimiser was compiled without GTK and " "cairo support.\n Stretch map will not be saved.\n"); -#endif /* HAVE_SAVE_TO_PNG */ +#endif /* CAN_SAVE_TO_PNG */ } @@ -2817,7 +2822,7 @@ int main(int argc, char *argv[]) strcat(command_line, buffer); } -#ifdef HAVE_SAVE_TO_PNG +#ifdef CAN_SAVE_TO_PNG #if !GLIB_CHECK_VERSION(2,36,0) g_type_init(); #endif diff --git a/src/hdfsee.c b/src/hdfsee.c index 954eb762..754dea1c 100644 --- a/src/hdfsee.c +++ b/src/hdfsee.c @@ -109,7 +109,7 @@ void hdfsee_window_closed(DisplayWindow *dw) main_n_windows--; - if ( main_n_windows == 0 ) gtk_exit(0); + if ( main_n_windows == 0 ) gtk_main_quit(); } diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 0794750f..87fd8fe7 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -90,7 +90,7 @@ struct sandbox struct sb_shm *shared; sem_t *queue_sem; - char *tmpdir; + const char *tmpdir; /* Final output */ Stream *stream; @@ -638,7 +638,6 @@ static void start_worker_process(struct sandbox *sb, int slot) } free(sb->fhs); free(sb->fds); - free(sb->tmpdir); free(sb->running); /* Not freed because it's not worth passing them down just for * this purpose: event list file handle, @@ -858,10 +857,11 @@ static void delete_temporary_folder(const char *tmpdir, int n_proc) /* List of files which it's safe to delete */ char *files[] = {"gmon.out", "mosflm.lp", "SUMMARY", "XDS.INP", "xfel_001.img", "xfel_001.spt", "xfel.drx", - "xfel.felix", "xfel.gve", "xfel.ini", "xfel.log"}; + "xfel.felix", "xfel.gve", "xfel.ini", "xfel.log", + "IDXREF.LP", "SPOT.XDS", "xfel.newmat", "XPARM.XDS"}; /* Number of items in the above list */ - int n_files = 11; + int n_files = 15; if ( n_proc > 99999 ) return; /* Paranoia */ @@ -870,6 +870,11 @@ static void delete_temporary_folder(const char *tmpdir, int n_proc) path = calloc(len+64, 1); if ( (workerdir == NULL) || (path == NULL) ) return; + snprintf(path, 127, "%s/mosflm.lp", tmpdir); + unlink(path); + snprintf(path, 127, "%s/SUMMARY", tmpdir); + unlink(path); + for ( slot=0; slot<n_proc; slot++ ) { struct stat s; @@ -883,17 +888,62 @@ static void delete_temporary_folder(const char *tmpdir, int n_proc) unlink(path); } - rmdir(workerdir); + if ( rmdir(workerdir) ) { + ERROR("Failed to delete worker temporary folder: %s\n", + strerror(errno)); + } + + } + + if ( rmdir(tmpdir) ) { + ERROR("Failed to delete temporary folder: %s\n", strerror(errno)); + } +} + + +char *create_tempdir(const char *temp_location) +{ + char *tmpdir; + size_t ll; + struct stat s; + + if ( temp_location == NULL ) { + temp_location = ""; + } + + ll = 64+strlen(temp_location); + tmpdir = malloc(ll); + if ( tmpdir == NULL ) { + ERROR("Failed to allocate temporary directory name\n"); + return NULL; + } + snprintf(tmpdir, ll, "%s/indexamajig.%i", temp_location, getpid()); + + if ( stat(tmpdir, &s) == -1 ) { + + int r; + + if ( errno != ENOENT ) { + ERROR("Failed to stat temporary folder.\n"); + return NULL; + } + + r = mkdir(tmpdir, S_IRWXU); + if ( r ) { + ERROR("Failed to create temporary folder: %s\n", + strerror(errno)); + return NULL; + } } - rmdir(tmpdir); + return tmpdir; } void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, - Stream *stream, const char *tempdir, int serial_start) + Stream *stream, const char *tmpdir, int serial_start) { int i; struct sandbox *sb; @@ -923,6 +973,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->n_proc = n_proc; sb->iargs = iargs; sb->serial = serial_start; + sb->tmpdir = tmpdir; sb->fds = NULL; sb->fhs = NULL; @@ -957,36 +1008,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, return; } - if ( tempdir == NULL ) { - tempdir = ""; - } - - ll = 64+strlen(tempdir); - sb->tmpdir = malloc(ll); - if ( sb->tmpdir == NULL ) { - ERROR("Failed to allocate temporary directory name\n"); - return; - } - snprintf(sb->tmpdir, ll, "%s/indexamajig.%i", tempdir, getpid()); - - if ( stat(sb->tmpdir, &s) == -1 ) { - - int r; - - if ( errno != ENOENT ) { - ERROR("Failed to stat temporary folder.\n"); - return; - } - - r = mkdir(sb->tmpdir, S_IRWXU); - if ( r ) { - ERROR("Failed to create temporary folder: %s\n", - strerror(errno)); - return; - } - - } - /* Fill the queue */ pthread_mutex_lock(&sb->shared->queue_lock); sb->shared->n_events = 0; @@ -1062,8 +1083,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, } while ( !allDone ); - if ( iargs->profile ) time_accounts_print(taccs); - fclose(fh); /* Indicate to the workers that we are finished, and wake them up one @@ -1116,7 +1135,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->shared->n_crystals); delete_temporary_folder(sb->tmpdir, n_proc); - free(sb->tmpdir); munmap(sb->shared, sizeof(struct sb_shm)); free(sb); diff --git a/src/im-sandbox.h b/src/im-sandbox.h index 06c9a33c..1412d199 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -71,6 +71,8 @@ struct sb_shm int n_crystals; }; +extern char *create_tempdir(const char *temp_location); + extern void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, Stream *stream, const char *tempdir, int serial_start); diff --git a/src/indexamajig.c b/src/indexamajig.c index 4a790a17..c0334d29 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -64,6 +64,7 @@ #include "integration.h" #include "taketwo.h" #include "im-sandbox.h" +#include "image.h" static void show_help(const char *s) @@ -110,13 +111,16 @@ static void show_help(const char *s) " (peakfinder8 only) Default: 0\n" " --max-res=<n> Maximum resolution for peak search (in pixels)\n" " (peakfinder8 only) Default: 1200\n" -" --min-snr-biggest-pix=<n> (peakFinder9 only) min snr of the biggest pixel in " -" the peak\n" -" --min-snr-peak-pix=<n> (peakFinder9 only) min snr of a peak pixel\n" -" --min-sig=<n> (peakFinder9 only) minimum standard deviation of " -" the background\n" -" --min-peak-over-neighbour=<n> (peakFinder9 only) just for speed. Biggest pixel" -" in peak must be n higher than this.\n" +" --min-snr-biggest-pix=<n>\n" +" Minimum snr of the biggest pixel in the peak\n" +" (peakfinder9 only)\n" +" --min-snr-peak-pix=<n>\n" +" Minimum snr of a peak pixel (peakfinder9 only)\n" +" --min-sig=<n> Minimum standard deviation of the background\n" +" (peakfinder9 only)\n" +" --min-peak-over-neighbour=<n>\n" +" Just for speed. Biggest pixel in peak must be n\n" +" higher than this (peakfinder9 only).\n" " --no-use-saturated Reject saturated peaks\n" " --no-revalidate Don't re-integrate and check HDF5 peaks\n" " --no-half-pixel-shift\n" @@ -171,6 +175,26 @@ static void show_help(const char *s) " Default: 30\n" " --felix-tthrange-min Minimum 2theta to consider for indexing (degrees)\n" " Default: 0\n" +"\n" +" --xgandalf-sampling-pitch\n" +" Sampling pitch: 0 (loosest) to 4 (most dense)\n" +" or with secondary Miller indices: 5 (loosest) to\n" +" 7 (most dense). Default: 6\n" +" --xgandalf-grad-desc-iterations\n" +" Gradient descent iterations: 0 (few) to 5 (many)\n" +" Default: 4\n" +" --xgandalf-tolerance Relative tolerance of the lattice vectors.\n" +" Default is 0.02\n" +" --xgandalf-no-deviation-from-provided-cell\n" +" Force the fitted cell to have the same lattice\n" +" parameters as the provided one\n" +" --xgandalf-min-lattice-vector-length\n" +" Minimum possible lattice vector length in A.\n" +" Default: 30 A\n" +" --xgandalf-max-lattice-vector-length\n" +" Maximum possible lattice vector length in A.\n" +" Default: 250 A\n" +"\n" "\nIntegration options:\n\n" " --integration=<meth> Integration method (rings,prof2d)-(cen,nocen)\n" " Default: rings-nocen\n" @@ -217,6 +241,66 @@ 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; i<s->n; 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; i<s->n; i++ ) { + s->weights[i] /= w_sum; + } + + return s; +} + + int main(int argc, char *argv[]) { int c; @@ -237,7 +321,10 @@ int main(int argc, char *argv[]) char *intrad = NULL; char *pkrad = NULL; char *int_str = NULL; - char *tempdir = NULL; + char *temp_location = NULL; /* e.g. /tmp */ + char *tmpdir; /* e.g. /tmp/indexamajig.12345 */ + char *rn; /* e.g. /home/taw/indexing */ + int r; char *int_diag = NULL; char *geom_filename = NULL; struct beam_params beam; @@ -250,6 +337,7 @@ int main(int argc, char *argv[]) int if_multi = 0; int if_retry = 1; int serial_start = 1; + char *spectrum_fn = NULL; /* Defaults */ iargs.cell = NULL; @@ -310,6 +398,12 @@ int main(int argc, char *argv[]) iargs.taketwo_opts.len_tol = -1.0; iargs.taketwo_opts.angle_tol = -1.0; iargs.taketwo_opts.trace_tol = -1.0; + iargs.xgandalf_opts.sampling_pitch = 6; + iargs.xgandalf_opts.grad_desc_iterations = 4; + iargs.xgandalf_opts.tolerance = 0.02; + iargs.xgandalf_opts.no_deviation_from_provided_cell = 0; + iargs.xgandalf_opts.minLatticeVectorLength_A = 30; + iargs.xgandalf_opts.maxLatticeVectorLength_A = 250; iargs.felix_opts.ttmin = -1.0; iargs.felix_opts.ttmax = -1.0; iargs.felix_opts.min_visits = 0; @@ -417,10 +511,23 @@ int main(int argc, char *argv[]) {"serial-start", 1, NULL, 344}, {"felix-domega", 1, NULL, 345}, {"felix-max-internal-angle", 1, NULL, 346}, - {"min-snr-biggest-pix" ,1, NULL,347}, - {"min-snr-peak-pix" ,1, NULL,348}, - {"min-sig" ,1, NULL,349}, - {"min-peak-over-neighbour" ,1, NULL,350}, + {"min-snr-biggest-pix", 1, NULL, 347}, + {"min-snr-peak-pix", 1, NULL, 348}, + {"min-sig", 1, NULL, 349}, + {"min-peak-over-neighbour", 1, NULL, 350}, + {"xgandalf-sampling-pitch", 1, NULL, 351}, + {"xgandalf-sps", 1, NULL, 351}, + {"xgandalf-grad-desc-iterations", 1, NULL, 352}, + {"xgandalf-gdis", 1, NULL, 352}, + {"xgandalf-tolerance", 1, NULL, 353}, + {"xgandalf-tol", 1, NULL, 353}, + {"xgandalf-no-deviation-from-provided-cell", 0, NULL, 354}, + {"xgandalf-ndfpc", 0, NULL, 354}, + {"xgandalf-min-lattice-vector-length", 1, NULL, 355}, + {"xgandalf-min-lvl", 1, NULL, 355}, + {"xgandalf-max-lattice-vector-length", 1, NULL, 356}, + {"xgandalf-max-lvl", 1, NULL, 356}, + {"spectrum-file", 1, NULL, 357}, {0, 0, NULL, 0} }; @@ -539,7 +646,7 @@ int main(int argc, char *argv[]) break; case 317 : - tempdir = strdup(optarg); + temp_location = strdup(optarg); break; case 318 : @@ -763,6 +870,58 @@ int main(int argc, char *argv[]) iargs.min_peak_over_neighbour = strtof(optarg, NULL); break; + case 351: + if (sscanf(optarg, "%u", &iargs.xgandalf_opts.sampling_pitch) != 1) + { + ERROR("Invalid value for --xgandalf-sampling-pitch\n"); + return 1; + } + break; + + case 352: + if (sscanf(optarg, "%u", &iargs.xgandalf_opts.grad_desc_iterations) != 1) + { + ERROR("Invalid value for --xgandalf-grad-desc-iterations\n"); + return 1; + } + break; + + case 353: + if (sscanf(optarg, "%f", &iargs.xgandalf_opts.tolerance) != 1) + { + ERROR("Invalid value for --xgandalf-tolerance\n"); + return 1; + } + break; + + case 354: + iargs.xgandalf_opts.no_deviation_from_provided_cell = 1; + break; + + case 355: + if (sscanf(optarg, "%f", + &iargs.xgandalf_opts.minLatticeVectorLength_A) != 1) + { + ERROR("Invalid value for " + "--xgandalf-min-lattice-vector-length\n"); + return 1; + } + break; + + case 356: + if (sscanf(optarg, "%f", + &iargs.xgandalf_opts.maxLatticeVectorLength_A) != 1) + { + ERROR("Invalid value for " + "--xgandalf-max-lattice-vector-length\n"); + return 1; + } + break; + + case 357: + spectrum_fn = strdup(optarg); + break; + case 0 : break; @@ -791,8 +950,8 @@ int main(int argc, char *argv[]) return 1; } - if ( tempdir == NULL ) { - tempdir = strdup("."); + if ( temp_location == NULL ) { + temp_location = strdup("."); } /* Open input */ @@ -950,6 +1109,19 @@ int main(int argc, char *argv[]) iargs.cell = NULL; } + /* Load spectrum from file if given */ + if ( spectrum_fn != NULL ) { + iargs.spectrum = read_spectrum_fromfile(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; + } + /* Parse integration diagnostic */ if ( int_diag != NULL ) { @@ -994,6 +1166,19 @@ int main(int argc, char *argv[]) } + tmpdir = create_tempdir(temp_location); + if ( tmpdir == NULL ) return 1; + + /* Change into temporary folder, temporarily, to control the crap + * dropped by indexing programs during setup */ + rn = getcwd(NULL, 0); + r = chdir(tmpdir); + if ( r ) { + ERROR("Failed to chdir to temporary folder: %s\n", + strerror(errno)); + return 1; + } + if ( indm_str == NULL ) { STATUS("No indexing methods specified. I will try to "); @@ -1001,6 +1186,7 @@ int main(int argc, char *argv[]) STATUS("To disable auto-detection of indexing methods, specify "); STATUS("which methods to use with --indexing=<methods>.\n"); STATUS("Use --indexing=none to disable indexing and integration.\n"); + indm_str = detect_indexing_methods(iargs.cell); } @@ -1060,6 +1246,7 @@ int main(int argc, char *argv[]) iargs.ipriv = setup_indexing(indm_str, iargs.cell, iargs.det, iargs.tols, flags, &iargs.taketwo_opts, + &iargs.xgandalf_opts, &iargs.felix_opts); if ( iargs.ipriv == NULL ) { ERROR("Failed to set up indexing system\n"); @@ -1068,6 +1255,16 @@ int main(int argc, char *argv[]) } + /* Change back to where we were before. Sandbox code will create + * worker subdirs inside the temporary folder, and process_image will + * change into them. */ + r = chdir(rn); + if ( r ) { + ERROR("Failed to chdir: %s\n", strerror(errno)); + return 1; + } + free(rn); + /* Open output stream */ st = open_stream_for_write_4(outfile, geom_filename, iargs.cell, argc, argv, indm_str); @@ -1081,13 +1278,14 @@ int main(int argc, char *argv[]) gsl_set_error_handler_off(); create_sandbox(&iargs, n_proc, prefix, config_basename, fh, - st, tempdir, serial_start); + st, tmpdir, serial_start); free_imagefile_field_list(iargs.copyme); cell_free(iargs.cell); free(iargs.beam->photon_energy_from); free(prefix); - free(tempdir); + free(temp_location); + free(tmpdir); free_detector_geometry(iargs.det); free(iargs.hdf5_peak_path); close_stream(st); diff --git a/src/merge.c b/src/merge.c index da333f6a..1dec13fd 100644 --- a/src/merge.c +++ b/src/merge.c @@ -3,11 +3,11 @@ * * Parallel weighted merging of intensities * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -46,6 +46,7 @@ #include "cell.h" #include "utils.h" #include "reflist.h" +#include "reflist-utils.h" #include "cell-utils.h" @@ -89,6 +90,17 @@ static void *create_merge_job(void *vqargs) } +static int alloc_contribs(struct reflection_contributions *c) +{ + c->contribs = realloc(c->contribs, c->max_contrib*sizeof(Reflection *)); + c->contrib_crystals = realloc(c->contrib_crystals, + c->max_contrib*sizeof(Crystal *)); + if ( c->contribs == NULL ) return 1; + if ( c->contrib_crystals == NULL ) return 1; + return 0; +} + + /* Find reflection hkl in 'list', creating it if it's not there, under * protection of 'lock' and returning a locked reflection */ static Reflection *get_locked_reflection(RefList *list, pthread_rwlock_t *lock, @@ -109,12 +121,31 @@ static Reflection *get_locked_reflection(RefList *list, pthread_rwlock_t *lock, * So, we must check again */ f = find_refl(list, h, k, l); if ( f == NULL ) { + + struct reflection_contributions *c; + f = add_refl(list, h, k, l); lock_reflection(f); pthread_rwlock_unlock(lock); set_intensity(f, 0.0); set_temp1(f, 0.0); set_temp2(f, 0.0); + + c = malloc(sizeof(struct reflection_contributions)); + if ( c != NULL ) { + c->n_contrib = 0; + c->max_contrib = 32; + c->contribs = NULL; + c->contrib_crystals = NULL; + if ( alloc_contribs(c) ) { + set_contributions(f, NULL); + } else { + set_contributions(f, c); + } + } else { + set_contributions(f, NULL); + } + } else { /* Someone else created it */ lock_reflection(f); @@ -160,6 +191,7 @@ static void run_merge_job(void *vwargs, int cookie) signed int h, k, l; double mean, sumweight, M2, temp, delta, R; double corr, res, w; + struct reflection_contributions *c; if ( get_partiality(refl) < MIN_PART_MERGE ) continue; @@ -212,6 +244,18 @@ static void run_merge_job(void *vwargs, int cookie) set_temp2(f, M2 + sumweight * delta * R); set_temp1(f, temp); set_redundancy(f, get_redundancy(f)+1); + + /* Record this contribution */ + c = get_contributions(f); + if ( c != NULL ) { + c->contribs[c->n_contrib] = refl; + c->contrib_crystals[c->n_contrib++] = cr; + if ( c->n_contrib == c->max_contrib ) { + c->max_contrib += 64; + alloc_contribs(c); + } + } /* else, too bad! */ + unlock_reflection(f); wargs->n_reflections++; @@ -291,9 +335,18 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads, Reflection *r2; get_indices(refl, &h, &k, &l); - r2 = add_refl(full2, h, k, l); + r2 = add_refl(full2, h, k, l); copy_data(r2, refl); + } else { + + /* We do not need the contribution list any more */ + struct reflection_contributions *c; + c = get_contributions(refl); + free(c->contribs); + free(c->contrib_crystals); + free(c); + } } diff --git a/src/multihistogram.c b/src/multihistogram.c index cf2b8c2b..ba1f7e76 100644 --- a/src/multihistogram.c +++ b/src/multihistogram.c @@ -93,7 +93,7 @@ void multihistogram_free(MultiHistogram *hi) } -void multihistogram_add_value(MultiHistogram *hi, double val, int cat) +void multihistogram_add_value(MultiHistogram *hi, double val, unsigned int cat) { int i, j; @@ -104,7 +104,7 @@ void multihistogram_add_value(MultiHistogram *hi, double val, int cat) if ( j >= hi->n_bins ) j = hi->n_bins - 1; for ( i=0; i<32; i++ ) { - if ( cat & 1<<i ) hi->bins[i][j]++; + if ( cat & (unsigned)1<<i ) hi->bins[i][j]++; } } diff --git a/src/multihistogram.h b/src/multihistogram.h index e4c2ae89..0949d259 100644 --- a/src/multihistogram.h +++ b/src/multihistogram.h @@ -40,7 +40,8 @@ extern MultiHistogram *multihistogram_new(); extern void multihistogram_free(MultiHistogram *hi); extern void multihistogram_delete_all_values(MultiHistogram *hi); -extern void multihistogram_add_value(MultiHistogram *hi, double val, int cat); +extern void multihistogram_add_value(MultiHistogram *hi, double val, + unsigned int cat); extern void multihistogram_set_min(MultiHistogram *hi, double min); extern void multihistogram_set_max(MultiHistogram *hi, double max); diff --git a/src/partial_sim.c b/src/partial_sim.c index 218780fd..012b0a59 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -844,7 +844,7 @@ int main(int argc, char *argv[]) image.num_peaks = 0; image.num_saturated_peaks = 0; image.spectrum_size = 0; - image.spectrum = NULL; + image.spectrum0 = NULL; image.serial = 0; image.event = NULL; diff --git a/src/partialator.c b/src/partialator.c index 6d10bc37..86431e4c 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -199,12 +199,14 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile, STATUS("Writing two-way split to %s ", tmp); write_reflist_2(tmp, split, sym); + free_contribs(split); reflist_free(split); snprintf(tmp, 1024, "%s2", outfile); split = merge_intensities(crystals2, n_crystals2, nthreads, min_measurements, push_res, 1, 0); STATUS("and %s\n", tmp); write_reflist_2(tmp, split, sym); + free_contribs(split); reflist_free(split); } @@ -291,6 +293,7 @@ static void write_custom_split(struct custom_split *csplit, int dsn, split = merge_intensities(crystalsn, n_crystalsn, nthreads, min_measurements, push_res, 1, 0); write_reflist_2(tmp, split, sym); + free_contribs(split); reflist_free(split); write_split(crystalsn, n_crystalsn, tmp, nthreads, pmodel, @@ -318,10 +321,12 @@ static void show_help(const char *s) " --no-scale Disable scale factor (G, B) refinement.\n" " --no-Bscale Disable B factor scaling.\n" " --no-pr Disable orientation/physics refinement.\n" +" --no-deltacchalf Disable rejection based on deltaCChalf.\n" " -m, --model=<model> Specify partiality model.\n" " --min-measurements=<n> Minimum number of measurements to require.\n" " --no-polarisation Disable polarisation correction.\n" " --max-adu=<n> Saturation value of detector.\n" +" --min-res=<n> Merge only crystals which diffract above <n> A.\n" " --push-res=<n> Merge higher than apparent resolution cutoff.\n" " -j <n> Run <n> analyses in parallel.\n" " --no-free Disable cross-validation (testing only).\n" @@ -702,9 +707,15 @@ static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals, static void all_residuals(Crystal **crystals, int n_crystals, RefList *full, double *presidual, double *pfree_residual, - double *plog_residual, double *pfree_log_residual) + double *plog_residual, double *pfree_log_residual, + int *pn_used) { int i; + int n_used = 0; + int n_nan_linear = 0; + int n_nan_linear_free = 0; + int n_nan_log = 0; + int n_nan_log_free = 0; *presidual = 0.0; *pfree_residual = 0.0; @@ -723,6 +734,11 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full, log_r = log_residual(crystals[i], full, 0, NULL, NULL); free_log_r = log_residual(crystals[i], full, 1, NULL, NULL); + if ( isnan(r) ) n_nan_linear++; + if ( isnan(free_r) ) n_nan_linear_free++; + if ( isnan(log_r) ) n_nan_log++; + if ( isnan(free_log_r) ) n_nan_log_free++; + if ( isnan(r) || isnan(free_r) || isnan(log_r) || isnan(free_log_r) ) continue; @@ -730,7 +746,28 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full, *pfree_residual += free_r; *plog_residual += log_r; *pfree_log_residual += free_log_r; + + n_used++; + } + + if ( n_nan_linear ) { + ERROR("WARNING: %i crystals had NaN linear residuals\n", + n_nan_linear); + } + if ( n_nan_linear_free ) { + ERROR("WARNING: %i crystals had NaN linear free residuals\n", + n_nan_linear_free); + } + if ( n_nan_log ) { + ERROR("WARNING: %i crystals had NaN log residuals\n", + n_nan_log); } + if ( n_nan_log_free ) { + ERROR("WARNING: %i crystals had NaN log free residuals\n", + n_nan_log_free); + } + + *pn_used = n_used; } @@ -738,13 +775,15 @@ static void show_all_residuals(Crystal **crystals, int n_crystals, RefList *full) { double dev, free_dev, log_dev, free_log_dev; + int n; all_residuals(crystals, n_crystals, full, - &dev, &free_dev, &log_dev, &free_log_dev); + &dev, &free_dev, &log_dev, &free_log_dev, &n); STATUS("Residuals:" " linear linear free log log free\n"); STATUS(" "); - STATUS("%15e %15e %15e %15e\n", dev, free_dev, log_dev, free_log_dev); + STATUS("%15e %15e %15e %15e (%i crystals)\n", + dev, free_dev, log_dev, free_log_dev, n); } @@ -875,6 +914,9 @@ int main(int argc, char *argv[]) double force_radius = -1.0; char *audit_info; int scaleflags = 0; + double min_res = 0.0; + int do_write_logs = 0; + int no_deltacchalf = 0; /* Long options */ const struct option longopts[] = { @@ -901,6 +943,7 @@ int main(int argc, char *argv[]) {"operator", 1, NULL, 9}, {"force-bandwidth", 1, NULL, 10}, {"force-radius", 1, NULL, 11}, + {"min-res", 1, NULL, 12}, {"no-scale", 0, &no_scale, 1}, {"no-Bscale", 0, &no_Bscale, 1}, @@ -912,6 +955,7 @@ int main(int argc, char *argv[]) {"no-free", 0, &no_free, 1}, {"output-every-cycle", 0, &output_everycycle, 1}, {"no-logs", 0, &no_logs, 1}, + {"no-deltacchalf", 0, &no_deltacchalf, 1}, {0, 0, NULL, 0} }; @@ -1059,6 +1103,16 @@ int main(int argc, char *argv[]) force_radius *= 1e9; break; + case 12 : + errno = 0; + min_res = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value for --min-res.\n"); + return 1; + } + min_res = 1e10/min_res; + break; + case 0 : break; @@ -1159,7 +1213,14 @@ int main(int argc, char *argv[]) scaleflags |= SCALE_NO_B; } - if ( !no_logs ) { + /* Decide whether or not to create stuff in the pr-logs folder */ + if ( !(no_logs || (no_pr && pmodel == PMODEL_UNITY)) ) { + do_write_logs = 1; + } else { + do_write_logs = 0; + } + + if ( do_write_logs ) { int r = mkdir("pr-logs", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); if ( r ) { if ( errno == EEXIST ) { @@ -1170,6 +1231,12 @@ int main(int argc, char *argv[]) return 1; } } + } else { + struct stat s; + if ( stat("pr-logs", &s) != -1 ) { + ERROR("WARNING: pr-logs folder exists, but I will not " + "write anything in it with these settings.\n"); + } } /* Read the custom split list (if applicable) */ @@ -1251,6 +1318,8 @@ int main(int argc, char *argv[]) n_crystals_seen++; if ( n_crystals_seen <= start_after ) continue; + if ( crystal_get_resolution_limit(cur.crystals[i]) < min_res ) continue; + crystals_new = realloc(crystals, (n_crystals+1)*sizeof(Crystal *)); if ( crystals_new == NULL ) { @@ -1349,9 +1418,10 @@ int main(int argc, char *argv[]) } /* Check rejection and write figures of merit */ - check_rejection(crystals, n_crystals, full, max_B); + check_rejection(crystals, n_crystals, full, max_B, no_deltacchalf); show_all_residuals(crystals, n_crystals, full); - if ( !no_pr && !no_logs ) { + + if ( do_write_logs ) { write_pgraph(full, crystals, n_crystals, 0, ""); write_logs_parallel(crystals, n_crystals, full, 0, nthreads, scaleflags); @@ -1369,6 +1439,7 @@ int main(int argc, char *argv[]) /* Create new reference if needed */ if ( reference == NULL ) { + free_contribs(full); reflist_free(full); if ( !no_scale ) { scale_all(crystals, n_crystals, nthreads, @@ -1379,10 +1450,10 @@ int main(int argc, char *argv[]) push_res, 1, 0); } /* else full still equals reference */ - check_rejection(crystals, n_crystals, full, max_B); + check_rejection(crystals, n_crystals, full, max_B, no_deltacchalf); show_all_residuals(crystals, n_crystals, full); - if ( !no_logs ) { + if ( do_write_logs ) { write_pgraph(full, crystals, n_crystals, i+1, ""); } @@ -1417,6 +1488,7 @@ int main(int argc, char *argv[]) /* Final merge */ STATUS("Final merge...\n"); if ( reference == NULL ) { + free_contribs(full); reflist_free(full); if ( !no_scale ) { scale_all(crystals, n_crystals, nthreads, scaleflags); @@ -1431,7 +1503,7 @@ int main(int argc, char *argv[]) /* Write final figures of merit (no rejection any more) */ show_all_residuals(crystals, n_crystals, full); - if ( !no_pr && !no_logs ) { + if ( do_write_logs ) { write_pgraph(full, crystals, n_crystals, -1, ""); write_logs_parallel(crystals, n_crystals, full, -1, nthreads, scaleflags); @@ -1463,6 +1535,7 @@ int main(int argc, char *argv[]) reflist_free(crystal_get_reflections(crystals[i])); crystal_free(crystals[i]); } + free_contribs(full); reflist_free(full); free_symoplist(sym); free(outfile); diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 6c9eb880..5b25100e 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -368,6 +368,7 @@ int main(int argc, char *argv[]) char *outfile = NULL; char *geometry = NULL; char *spectrum_str = NULL; + char *spectrum_fn = NULL; GradientMethod grad; SpectrumType spectrum_type; int ndone = 0; /* Number of simulations done (images or not) */ @@ -429,6 +430,7 @@ int main(int argc, char *argv[]) {"photon-energy", 1, NULL, 9}, {"nphotons", 1, NULL, 10}, {"beam-radius", 1, NULL, 11}, + {"spectrum-file", 1, NULL, 12}, {0, 0, NULL, 0} }; @@ -590,6 +592,9 @@ int main(int argc, char *argv[]) } break; + case 12 : + spectrum_fn = strdup(optarg); + break; case 0 : break; @@ -695,6 +700,8 @@ int main(int argc, char *argv[]) strcasecmp(spectrum_str, "twocolours") == 0 || strcasecmp(spectrum_str, "twocolors") == 0) { spectrum_type = SPECTRUM_TWOCOLOUR; + } else if ( strcasecmp(spectrum_str, "fromfile") == 0) { + spectrum_type = SPECTRUM_FROMFILE; } else { ERROR("Unrecognised spectrum type '%s'\n", spectrum_str); return 1; @@ -812,7 +819,7 @@ int main(int argc, char *argv[]) powder.det = image.det; powder.beam = NULL; powder.lambda = 0.0; - powder.spectrum = NULL; + powder.spectrum0 = NULL; powder.dp = malloc(image.det->n_panels*sizeof(float *)); if ( powder.dp == NULL ) { ERROR("Failed to allocate powder data\n"); @@ -855,6 +862,11 @@ int main(int argc, char *argv[]) STATUS(" X-ray spectrum: two colour, " "separation %.5f %%\n", image.bw*100.0); break; + + case SPECTRUM_FROMFILE: + STATUS(" X-ray spectrum: from %s\n", + spectrum_fn); + break; } if ( random_size ) { STATUS(" Crystal size: random, between " @@ -966,15 +978,20 @@ int main(int argc, char *argv[]) switch ( spectrum_type ) { case SPECTRUM_TOPHAT : - image.spectrum = generate_tophat(&image); + image.spectrum0 = generate_tophat(&image); break; case SPECTRUM_SASE : - image.spectrum = generate_SASE(&image, rng); + image.spectrum0 = generate_SASE(&image, rng); break; case SPECTRUM_TWOCOLOUR : - image.spectrum = generate_twocolour(&image); + image.spectrum0 = generate_twocolour(&image); + break; + + case SPECTRUM_FROMFILE : + image.spectrum0 = generate_spectrum_fromfile(&image, + spectrum_fn); break; } diff --git a/src/post-refinement.c b/src/post-refinement.c index bde9c02c..0a799b9e 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -69,8 +69,8 @@ const char *str_prflag(enum prflag flag) case PRFLAG_EARLY : return "early rejection"; - case PRFLAG_CC : - return "low CC"; + case PRFLAG_DELTACCHALF : + return "negative delta CC½"; case PRFLAG_BIGB : return "B too big"; diff --git a/src/post-refinement.h b/src/post-refinement.h index 71a6d7f3..e2873b76 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -50,7 +50,7 @@ enum prflag PRFLAG_FEWREFL = 16, PRFLAG_SOLVEFAIL = 17, PRFLAG_EARLY = 18, - PRFLAG_CC = 19, + PRFLAG_DELTACCHALF = 19, PRFLAG_BIGB = 20, PRFLAG_SCALEBAD = 21, }; diff --git a/src/process_image.c b/src/process_image.c index b0a03e46..10863fd3 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -272,6 +272,11 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, image.bw = 0.00000001; } + /* Set beam spectrum for pink beam data */ + if ( iargs->spectrum != NULL ) { + image.spectrum = iargs->spectrum; + } + if ( image_feature_count(image.features) < iargs->min_peaks ) { r = chdir(rn); if ( r ) { @@ -280,6 +285,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, return; } free(rn); + image.hit = 0; if ( iargs->stream_nonhits ) { goto streamwrite; @@ -287,6 +293,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, goto out; } } + image.hit = 1; /* Index the pattern */ time_accounts_set(taccs, TACC_INDEXING); @@ -317,7 +324,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, if ( iargs->fix_profile_r < 0.0 ) { for ( i=0; i<image.n_crystals; i++ ) { - refine_radius(image.crystals[i], &image); + if ( refine_radius(image.crystals[i], &image) ) { + ERROR("WARNING: Radius determination failed\n"); + } } } diff --git a/src/process_image.h b/src/process_image.h index 52f063fd..d4dcbf04 100644 --- a/src/process_image.h +++ b/src/process_image.h @@ -41,6 +41,7 @@ struct index_args; #include "im-sandbox.h" #include "time-accounts.h" #include "taketwo.h" +#include "xgandalf.h" #include "felix.h" @@ -109,7 +110,9 @@ struct index_args int overpredict; int profile; /* Whether or not to do wall clock profiling */ struct taketwo_options taketwo_opts; + struct xgandalf_options xgandalf_opts; struct felix_options felix_opts; + struct spectrum *spectrum; }; diff --git a/src/rejection.c b/src/rejection.c index 9d41f9b4..55546bbc 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -3,11 +3,11 @@ * * Crystal rejection for scaling * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -94,74 +94,231 @@ void early_rejection(Crystal **crystals, int n) } -static void check_cc(Crystal *cr, RefList *full) +static void calculate_refl_mean_var(RefList *full) { - RefList *list = crystal_get_reflections(cr); Reflection *refl; RefListIterator *iter; - double G = crystal_get_osf(cr); - double B = crystal_get_Bfac(cr); - UnitCell *cell = crystal_get_cell(cr); - double *vec1, *vec2; - int n, i; - double cc; - - n = num_reflections(list); - vec1 = malloc(n*sizeof(double)); - vec2 = malloc(n*sizeof(double)); - if ( (vec1 == NULL) || (vec2 == NULL) ) { - ERROR("Not enough memory to check CCs\n"); - return; - } + signed int oh = 0; + signed int ok = 0; + signed int ol = 0; - i = 0; - for ( refl = first_refl(list, &iter); + /* Iterate over all reflections */ + for ( refl = first_refl(full, &iter); refl != NULL; refl = next_refl(refl, iter) ) { + struct reflection_contributions *c; + int j; signed int h, k, l; - double pobs, pcalc; - double res, corr, Ipart; - Reflection *match; - - if ( !get_flag(refl) ) continue; /* Not free-flagged */ + double res; + double K; + double Ex = 0.0; + double Ex2 = 0.0; get_indices(refl, &h, &k, &l); - res = resolution(cell, h, k, l); - if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue; - match = find_refl(full, h, k, l); - if ( match == NULL ) continue; + /* next_refl() will iterate over multiple copies of the same + * unique reflection, but we are only interested in seeing each + * unique reflection once. */ + if ( (h==oh) && (k==ok) && (l==ol) ) continue; + oh = h; ok = k; ol = l; + + /* We use the mean (merged) intensity as the reference point + * for shifting the data in the variance calculation */ + K = get_intensity(refl); + + c = get_contributions(refl); + assert(c != NULL); + + /* Calculate the resolution just once, using the cell from the + * first crystal to contribute, otherwise it takes too long */ + res = resolution(crystal_get_cell(c->contrib_crystals[0]), + h, k, l); + + /* Mean of contributions */ + for ( j=0; j<c->n_contrib; j++ ) { + + double Ii, G, B; + + G = crystal_get_osf(c->contrib_crystals[j]); + B = crystal_get_Bfac(c->contrib_crystals[j]); + Ii = get_intensity(c->contribs[j]); + + /* Total (multiplicative) correction factor */ + Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j]) + / get_partiality(c->contribs[j]); + + Ex += Ii - K; + Ex2 += (Ii - K) * (Ii - K); + + } + + if ( c->n_contrib < 2 ) continue; + + set_temp1(refl, Ex); + set_temp2(refl, Ex2); + } + +} + + +static double calculate_cchalf(RefList *template, RefList *full, + Crystal *exclude, int *pnref) +{ + Reflection *trefl; + RefListIterator *iter; + double sig2E, sig2Y; + int n = 0; + double wSum = 0.0; + double wSum2 = 0.0; + double mean = 0.0; + double S = 0.0; + double all_sum_var = 0.0; + signed int oh = 0; + signed int ok = 0; + signed int ol = 0; + + /* "template" is the list of reflections to be included in CChalf */ + for ( trefl = first_refl(template, &iter); + trefl != NULL; + trefl = next_refl(trefl, iter) ) + { + signed int h, k, l; + double refl_mean, refl_var; + double Ex, Ex2, K; + int n_removed = 0; + double w = 1.0; + double meanOld; + double res; + struct reflection_contributions *c; + Reflection *refl; + Reflection *exrefl; + + /* The values we need are stored in the "full" list, not the + * template list */ + get_indices(trefl, &h, &k, &l); + + if ( (h==oh) && (k==ok) && (l==ol) ) continue; + oh = h; ok = k; ol = l; + + refl = find_refl(full, h, k, l); + if ( refl == NULL ) continue; /* However, there might not have + * been enough measurements for + * it to appear in "full" */ + + /* We use the mean (merged) intensity as the reference point + * for shifting the data in the variance calculation */ + K = get_intensity(refl); + Ex = get_temp1(refl); + Ex2 = get_temp2(refl); + c = get_contributions(refl); + assert(c != NULL); + + /* Resolution from first contributing crystal, like above */ + res = resolution(crystal_get_cell(c->contrib_crystals[0]), + h, k, l); + + /* Remove contribution(s) from the excluded crystal. + * If the crystal is marked as bad, we should not remove it + * because it did not contribute in the first place. */ + if ( exclude != NULL && !crystal_get_user_flag(exclude) ) { + exrefl = find_refl(crystal_get_reflections(exclude), h, k, l); + } else { + exrefl = NULL; + } + + while ( exrefl != NULL ) { + + double G, B; + double Ii = get_intensity(exrefl); + + G = crystal_get_osf(exclude); + B = crystal_get_Bfac(exclude); - /* Calculated partiality */ - pcalc = get_partiality(refl); + /* Total (multiplicative) correction factor */ + Ii *= 1.0/G * exp(B*res*res) * get_lorentz(exrefl) + / get_partiality(exrefl); - /* Observed partiality */ - corr = G * exp(B*res*res) * get_lorentz(refl); - Ipart = get_intensity(refl) * corr; - pobs = Ipart / get_intensity(match); + /* Remove contribution of this reflection */ + Ex -= Ii - K; + Ex2 -= (Ii - K)*(Ii - K); + + n_removed++; + + exrefl = next_found_refl(exrefl); + } - vec1[i] = pobs; - vec2[i] = pcalc; - i++; + if ( c->n_contrib - n_removed < 2 ) continue; + + refl_mean = K + (Ex / (c->n_contrib - n_removed)); + refl_var = (Ex2 - (Ex*Ex)/(c->n_contrib - n_removed)) / (c->n_contrib - n_removed - 1); + refl_var /= (c->n_contrib - n_removed) / 2.0; + + all_sum_var += refl_var; + n++; + + /* Running variance calculation to get sig2Y */ + wSum += w; + wSum2 += w*w; + meanOld = mean; + mean = meanOld + (w/wSum) * (refl_mean - meanOld); + S += w * (refl_mean - meanOld) * (refl_mean - mean); } - cc = gsl_stats_correlation(vec1, 1, vec2, 1, i); - //printf("%f\n", cc); - if ( cc < 0.5 ) crystal_set_user_flag(cr, PRFLAG_CC); + sig2E = all_sum_var / n; + sig2Y = S / (wSum - 1.0); - free(vec1); - free(vec2); + if ( pnref != NULL ) { + *pnref = n; + } + return (sig2Y - 0.5*sig2E) / (sig2Y + 0.5*sig2E); } -static UNUSED void check_ccs(Crystal **crystals, int n_crystals, RefList *full) +static void check_deltacchalf(Crystal **crystals, int n, RefList *full) { + double cchalf; int i; + double *vals; + double mean, sd; + int nref = 0; + + calculate_refl_mean_var(full); - for ( i=0; i<n_crystals; i++ ) { - check_cc(crystals[i], full); + cchalf = calculate_cchalf(full, full, NULL, &nref); + STATUS("Overall CChalf = %f %% (%i reflections)\n", cchalf*100.0, nref); + + vals = malloc(n*sizeof(double)); + if ( vals == NULL ) { + ERROR("Not enough memory for deltaCChalf check\n"); + return; } + + for ( i=0; i<n; i++ ) { + double cchalf, cchalfi; + RefList *template = crystal_get_reflections(crystals[i]); + cchalf = calculate_cchalf(template, full, NULL, &nref); + cchalfi = calculate_cchalf(template, full, crystals[i], &nref); + //STATUS("Frame %i:", i); + //STATUS(" With = %f ", cchalf*100.0); + //STATUS("Without = %f", cchalfi*100.0); + //STATUS(" Delta = %f ", (cchalf - cchalfi)*100.0); + //STATUS("(nref = %i)\n", nref); + vals[i] = cchalf - cchalfi; + progress_bar(i, n-1, "Calculating deltaCChalf"); + } + + mean = gsl_stats_mean(vals, 1, n); + sd = gsl_stats_sd_m(vals, 1, n, mean); + STATUS("deltaCChalf = %f ± %f %%\n", mean*100.0, sd*100.0); + + for ( i=0; i<n; i++ ) { + if ( vals[i] < mean-2.0*sd ) { + crystal_set_user_flag(crystals[i], PRFLAG_DELTACCHALF); + } + } + + free(vals); } @@ -192,14 +349,16 @@ static void show_duds(Crystal **crystals, int n_crystals) } -void check_rejection(Crystal **crystals, int n, RefList *full, double max_B) +void check_rejection(Crystal **crystals, int n, RefList *full, double max_B, + int no_deltacchalf) { int i; int n_acc = 0; - /* Check according to CCs FIXME: Disabled */ - //if ( full != NULL ) check_ccs(crystals, n, full); - + /* Check according to delta CC½ */ + if ( !no_deltacchalf && (full != NULL) ) { + check_deltacchalf(crystals, n, full); + } for ( i=0; i<n; i++ ) { if ( fabs(crystal_get_Bfac(crystals[i])) > max_B ) { diff --git a/src/rejection.h b/src/rejection.h index 8979313b..5a0c3e9f 100644 --- a/src/rejection.h +++ b/src/rejection.h @@ -3,11 +3,11 @@ * * Crystal rejection for scaling * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -39,6 +39,6 @@ extern void early_rejection(Crystal **crystals, int n); extern void check_rejection(Crystal **crystals, int n, RefList *full, - double max_B); + double max_B, int no_deltacchalf); #endif /* REJECTION_H */ diff --git a/src/scaling.c b/src/scaling.c index dcfa7a59..35481982 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -46,6 +46,7 @@ #include "cell.h" #include "cell-utils.h" #include "scaling.h" +#include "reflist-utils.h" struct scale_args @@ -166,6 +167,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags) meanB /= n_crystals; STATUS("Mean B = %e\n", meanB); + free_contribs(full); reflist_free(full); niter++; diff --git a/src/time-accounts.c b/src/time-accounts.c index d3e989fa..15a984cd 100644 --- a/src/time-accounts.c +++ b/src/time-accounts.c @@ -173,8 +173,8 @@ static const char *taccname(enum timeaccount acc) case TACC_ENDCHECK : return "Checking end"; case TACC_WAKEUP : return "Waking up workers"; case TACC_WAITPID : return "Waiting on workers"; - case TACC_HDF5OPEN : return "Opening HDF5 file"; - case TACC_HDF5READ : return "Reading HDF5 file"; + case TACC_HDF5OPEN : return "Opening image file"; + case TACC_HDF5READ : return "Reading image file"; case TACC_FILTER : return "Image filters"; case TACC_RESRANGE : return "Resolution range"; case TACC_PEAKSEARCH : return "Peak search"; diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c index 8c666775..aee7a068 100644 --- a/tests/gpu_sim_check.c +++ b/tests/gpu_sim_check.c @@ -157,8 +157,8 @@ int main(int argc, char *argv[]) gpu_image.det = det; cpu_image.beam = NULL; gpu_image.beam = NULL; - cpu_image.spectrum = NULL; - gpu_image.spectrum = NULL; + cpu_image.spectrum0 = NULL; + gpu_image.spectrum0 = NULL; cpu_image.lambda = ph_en_to_lambda(eV_to_J(6000)); gpu_image.lambda = ph_en_to_lambda(eV_to_J(6000)); @@ -167,8 +167,8 @@ int main(int argc, char *argv[]) cpu_image.nsamples = 10; gpu_image.nsamples = 10; - cpu_image.spectrum = generate_tophat(&cpu_image); - gpu_image.spectrum = generate_tophat(&gpu_image); + cpu_image.spectrum0 = generate_tophat(&cpu_image); + gpu_image.spectrum0 = generate_tophat(&gpu_image); start = get_hires_seconds(); if ( get_diffraction_gpu(gctx, &gpu_image, 8, 8, 8, cell, 0, 0) ) { |