From 895bb2a9cf514bf81affa215bc1ea86bd42757fa Mon Sep 17 00:00:00 2001 From: Cornelius Gati Date: Tue, 5 Feb 2013 12:11:19 +0100 Subject: Fixed negative intensities in SPOT.XDS --- libcrystfel/src/xds.c | 402 ++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 355 insertions(+), 47 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c index 96c4aca3..f40ade79 100644 --- a/libcrystfel/src/xds.c +++ b/libcrystfel/src/xds.c @@ -69,18 +69,220 @@ struct xds_data { int rbufpos; int rbuflen; + /* High-level stuff */ + char newmatfile[128]; + char spotfile[128]; + int step; + int finished_ok; + UnitCell *target_cell; }; +static void xds_parseline(const char *line, struct image *image, + struct xds_data *xds) +{ + + char *copy; + int i; + + copy = strdup(line); + for ( i=0; ipty, xds->rbuffer+xds->rbufpos, + xds->rbuflen-xds->rbufpos); + if ( (rval == -1) || (rval == 0) ) return 1; + + xds->rbufpos += rval; + assert(xds->rbufpos <= xds->rbuflen); + + while ( (!no_string) && (xds->rbufpos > 0) ) { + + int i; + int block_ready = 0; + + /* See if there's a full line in the buffer yet */ + for ( i=0; irbufpos-1; i++ ) { + /* Means the last value looked at is rbufpos-2 */ + + if ( (xds->rbuffer[i] == '\r') + && (xds->rbuffer[i+1] == '\n') ) { + block_ready = 1; + break; + } + + } + + if ( block_ready ) { + + unsigned int new_rbuflen; + unsigned int endbit_length; + char *block_buffer = NULL; + + block_buffer = malloc(i+1); + memcpy(block_buffer, xds->rbuffer, i); + block_buffer[i] = '\0'; + + if ( block_buffer[0] == '\r' ) { + memmove(block_buffer, block_buffer+1, i); + } + + xds_parseline(block_buffer, image, xds); + free(block_buffer); + endbit_length = i+2; + + /* Now the block's been parsed, it should be + * forgotten about */ + memmove(xds->rbuffer, + xds->rbuffer + endbit_length, + xds->rbuflen - endbit_length); + + /* Subtract the number of bytes removed */ + xds->rbufpos = xds->rbufpos + - endbit_length; + new_rbuflen = xds->rbuflen - endbit_length; + if ( new_rbuflen == 0 ) new_rbuflen = 256; + xds->rbuffer = realloc(xds->rbuffer, + new_rbuflen); + xds->rbuflen = new_rbuflen; + + } else { + + if ( xds->rbufpos==xds->rbuflen ) { + + /* More buffer space is needed */ + xds->rbuffer = realloc( + xds->rbuffer, + xds->rbuflen + 256); + xds->rbuflen = xds->rbuflen + 256; + /* The new space gets used at the next + * read, shortly... */ + + } + no_string = 1; + + } + + } + + return 0; +} + + + + + + +static int read_newmat(const char *filename, struct image *image) +{ + FILE * fh; + float axstar, aystar, azstar, bxstar, bystar, bzstar, cxstar, cystar, czstar; + char asx[11], asy[11], asz[11]; + char bsx[11], bsy[11], bsz[11]; + char csx[11], csy[11], csz[11]; + char *rval, line[1024]; + int r; + + fh = fopen("IDXREF.LP", "r"); + if ( fh == NULL ) { + ERROR("Couldn't open 'IDXREF.LP'\n"); + return 1; + } + + do { + rval = fgets(line, 1023, fh); + + } while (strcmp(line, " # COORDINATES OF REC. BASIS VECTOR LENGTH 1/LENGTH\n") != 0 ); + + /*free line after chunk*/ + rval = fgets(line, 1023, fh); + rval = fgets(line, 1023, fh); + + memcpy(asx, line+7, 10); asx[10] = '\0'; + memcpy(asy, line+17, 10); asy[10] = '\0'; + memcpy(asz, line+27, 10); asz[10] = '\0'; + + rval = fgets(line, 1023, fh); + + memcpy(bsx, line+7, 10); bsx[10] = '\0'; + memcpy(bsy, line+17, 10); bsy[10] = '\0'; + memcpy(bsz, line+27, 10); bsz[10] = '\0'; + + rval = fgets(line, 1023, fh); + + 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"); + return 1;} + + fclose(fh); + + //printf("asx='%s'\n", asx); + //printf("asy='%s'\n", asy); + //printf("asz='%s'\n", asz); + //printf("cxstar='%f'\n", cxstar); + //printf("cystar='%f'\n", cystar); + //printf("czstar='%f'\n", czstar); + //printf("bsx='%s'\n", bsx); + //printf("bsy='%s'\n", bsy); + //printf("bsz='%s'\n", bsz); + //printf("axstar='%f'\n", axstar); + //printf("aystar='%f'\n", aystar); + //printf("azstar='%f'\n", azstar); + //printf("csx='%s'\n", csx); + //printf("csy='%s'\n", csy); + //printf("csz='%s'\n", csz); + //printf("bxstar='%f'\n", bxstar); + //printf("bystar='%f'\n", bystar); + //printf("bzstar='%f'\n", bzstar); + + image->candidate_cells[0] = cell_new(); + + cell_set_reciprocal(image->candidate_cells[0], + axstar*10e9, aystar*10e9, azstar*10e9, + bxstar*10e9, bystar*10e9, bzstar*10e9, + -cxstar*10e9, -cystar*10e9, -czstar*10e9); + image->ncells = 1; + + //cell_print(image->candidate_cells[0]); + + return 0; +} + + static void write_spot(struct image *image, const char *filename) { FILE *fh; int i; - double fclen = 67.8e-3; /* fake camera length in m */ + //double fclen = 170.0e-3; /* fake camera length in m */ int n; - fh = fopen(filename, "w"); + fh = fopen("SPOT.XDS", "w"); if ( !fh ) { - ERROR("Couldn't open temporary file '%s'\n", filename); + ERROR("Couldn't open temporary file '%s'\n", "SPOT.XDS"); return; } @@ -100,30 +302,38 @@ static void write_spot(struct image *image, const char *filename) xs = (f->fs-p->min_fs)*p->fsx + (f->ss-p->min_ss)*p->ssx; ys = (f->fs-p->min_fs)*p->fsy + (f->ss-p->min_ss)*p->ssy; - rx = (xs + p->cnx) / p->res; - ry = (ys + p->cny) / p->res; + rx = ((xs + p->cnx) / p->res); + ry = ((ys + p->cny) / p->res); + + //printf("xs=%f ys=%f ----> rx=%f ry=%f\n", xs, ys, rx, ry); + + //x = rx*fclen/p->clen; + //y = ry*fclen/p->clen; /* Peak positions in m */ - x = rx*fclen/p->clen; - y = ry*fclen/p->clen; /* Peak positions in m */ + //printf("x=%f y=%f\n", x, y); - fprintf(fh, "%10.2f %10.2f %10.2f %10.0f", - x*1e3, y*1e3, 0.0, 0.0); + x = (rx / 70e-6) + 1500; + y = (ry / 70e-6) + 1500; + + if (f->intensity <= 0) continue; + + fprintf(fh, "%10.2f %10.2f %10.2f %10.0f.\n", + x, y, 0.5, f->intensity); } fclose(fh); } -static char *write_INP(struct image *image) +static char *write_inp(struct image *image) { FILE *fh; char *filename; - //double tt; - + filename = malloc(1024); if ( filename == NULL ) return NULL; - snprintf(filename, 1023, filename, image->id); + snprintf(filename, 1023, "XDS.INP"); fh = fopen(filename, "w"); if ( !fh ) { @@ -132,59 +342,157 @@ static char *write_INP(struct image *image) return NULL; } - //get_q_for_panel(image->det->furthest_out_panel, - // image->det->furthest_out_fs, - // image->det->furthest_out_ss, - // &tt, 1.0/image->lambda); - fprintf(fh, "JOB= IDXREF\n"); - fprintf(fh, "ORGX= 1540.64\n"); - fprintf(fh, "ORGY= 1543.78\n"); - fprintf(fh, "DETECTOR_DISTANCE= 50.000\n"); - fprintf(fh, "OSCILLATION_RANGE= 0.100\n"); - fprintf(fh, "X-RAY_WAVELENGTH= 1.77000\n"); + fprintf(fh, "ORGX= 1500\n"); + fprintf(fh, "ORGY= 1500\n"); + fprintf(fh, "DETECTOR_DISTANCE= 65.00\n"); //IMPORTANT + fprintf(fh, "OSCILLATION_RANGE= 0.300\n"); + fprintf(fh, "X-RAY_WAVELENGTH= 6.2\n"); // IMPORTANT fprintf(fh, "NAME_TEMPLATE_OF_DATA_FRAMES=/home/dikay/insu/data_exp1/ins_ssad_1_???.img \n"); fprintf(fh, "DATA_RANGE=1 1\n"); fprintf(fh, "SPOT_RANGE=1 1\n"); - fprintf(fh, "SPACE_GROUP_NUMBER=0\n"); - fprintf(fh, "UNIT_CELL_CONSTANTS= 70 80 90 90 90 90\n"); - fprintf(fh, "NX= 3072\n"); - fprintf(fh, "NY= 3072\n"); - fprintf(fh, "QX= .073242\n"); - fprintf(fh, "QY= .073242\n"); + //fprintf(fh, "SPACE_GROUP_NUMBER=94\n"); //CatB 94 + //fprintf(fh, "UNIT_CELL_CONSTANTS= 125.4 125.4 54.56 90 90 90\n"); //CatB 125.4 125.4 54.56 90 90 90 + //fprintf(fh, "SPACE_GROUP_NUMBER=194\n"); //PS1 194 + //fprintf(fh, "UNIT_CELL_CONSTANTS=281 281 165.2 90 90 120\n"); //PS1 281 281 165.2 90 90 120 + fprintf(fh, "SPACE_GROUP_NUMBER= 0\n"); //LYS 96 + fprintf(fh, "UNIT_CELL_CONSTANTS= 0 0 0 0 0 0\n"); //LYS 77.32 77.32 38.16 90 90 90 + fprintf(fh, "NX= 3000\n"); + fprintf(fh, "NY= 3000\n"); + fprintf(fh, "QX= 0.07\n"); + fprintf(fh, "QY= 0.07\n"); + fprintf(fh, "INDEX_ORIGIN=0 0 0\n"); fprintf(fh, "DIRECTION_OF_DETECTOR_X-AXIS=1 0 0\n"); fprintf(fh, "DIRECTION_OF_DETECTOR_Y-AXIS=0 1 0\n"); fprintf(fh, "INCIDENT_BEAM_DIRECTION=0 0 1\n"); - fprintf(fh, "ROTATION_AXIS=1 0 0\n"); - fprintf(fh, "DETECTOR= CCDCHESS\n"); + fprintf(fh, "ROTATION_AXIS=0 1 0\n"); + fprintf(fh, "DETECTOR= CSPAD\n"); fprintf(fh, "MINIMUM_VALID_PIXEL_VALUE= 1\n"); - fprintf(fh, "OVERLOAD= 65500\n"); + fprintf(fh, "OVERLOAD= 200000000\n"); + fprintf(fh, "INDEX_ERROR= 0.2\n"); + //fprintf(fh, "INDEX_QUALITY= 0.5\n"); + //fprintf(fh, "REFINE(IDXREF)= ALL\n"); + //fprintf(fh, "MINIMUM_NUMBER_OF_PIXELS_IN_A_SPOT= 1\n"); + //fprintf(fh, "MAXIMUM_ERROR_OF_SPOT_POSITION= 20.0\n"); fclose(fh); return filename; +} void run_xds(struct image *image, UnitCell *cell) { - //unsigned int opts; - //int status; - //int rval; - struct xds_data *xds; - char *INP_filename; - //char gff_filename[1024]; + + unsigned int opts; + int status; + int rval; + struct xds_data *xds; + char *inp_filename; + + xds = malloc(sizeof(struct xds_data)); + if ( xds == NULL ) { + ERROR("Couldn't allocate memory for xds data.\n"); + return; + } + + xds->target_cell = cell; - write_INP(image); - INP_filename = write_INP(image); + write_inp(image); + inp_filename = write_inp(image); - if ( INP_filename == NULL ) { + if ( inp_filename == NULL ) { ERROR("Failed to write XDS.INP file for XDS.\n"); return; } - xds = malloc(sizeof(struct xds_data)); - if ( xds == NULL ) { - ERROR("Couldn't allocate memory for xds data.\n"); + snprintf(xds->spotfile, 127, "SPOT.XDS"); + write_spot(image, xds->spotfile); + + xds->pid = forkpty(&xds->pty, NULL, NULL, NULL); + + if ( xds->pid == -1 ) { + ERROR("Failed to fork for XDS\n"); + free(xds); return; - }} - //free(xds); -//} + } + if ( xds->pid == 0 ) { + + /* Child process: invoke XDS */ + struct termios t; + + /* Turn echo off */ + tcgetattr(STDIN_FILENO, &t); + t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); + tcsetattr(STDIN_FILENO, TCSANOW, &t); + + execlp("xds", "", (char *)NULL); + ERROR("Failed to invoke XDS.\n"); + _exit(0); + + } + + xds->rbuffer = malloc(256); + xds->rbuflen = 256; + xds->rbufpos = 0; + + /* Set non-blocking */ + opts = fcntl(xds->pty, F_GETFL); + fcntl(xds->pty, F_SETFL, opts | O_NONBLOCK); + + //xds->step = 1; /* This starts the "initialisation" procedure */ + xds->finished_ok = 0; + + do { + + fd_set fds; + struct timeval tv; + int sval; + + FD_ZERO(&fds); + FD_SET(xds->pty, &fds); + + tv.tv_sec = 30; + tv.tv_usec = 0; + + sval = select(xds->pty+1, &fds, NULL, NULL, &tv); + + if ( sval == -1 ) { + + const int err = errno; + + switch ( err ) { + + case EINTR: + STATUS("Restarting select()\n"); + break; + + default: + ERROR("select() failed: %s\n", strerror(err)); + rval = 1; + + } + + } else if ( sval != 0 ) { + rval = xds_readable(image, xds); + } else { + ERROR("No response from XDS..\n"); + rval = 1; + } + + } while ( !rval ); + + close(xds->pty); + free(xds->rbuffer); + waitpid(xds->pid, &status, 0); + + //FIXME// + //if ( xds->finished_ok == 1 ) { + // ERROR("XDS doesn't seem to be working properly.\n"); + //} else { + /* Read the XDS NEWMAT file and get cell if found */ + read_newmat(xds->newmatfile, image); + + //} + + free(xds); +} -- cgit v1.2.3