Generate XDS.INP: Difference between revisions

#!/bin/bash                                                                  

# purpose: generate XDS.INP                                                  

#                                                                            

# tested with some datasets from ALS, SSRL, SLS, ESRF, BESSY, SPring-8 and PF; only MAR, ADSC/SMV, PILATUS, Eiger, RAXIS (in-house), Bruker (PHOTON II) detectors;  

# for other detectors, values marked with XXX must be manually filled in.                                

#                                                                                                        

# revision 0.03 . Kay Diederichs 2/2010                                                                  

# revision 0.04 . Kay Diederichs 4/2010 - include alternative ORGX, ORGY calculations for ADSC            

# revision 0.05 . Kay Diederichs 5/2010 - grep for "Corrected" in addition to "marccd"; needed for BESSY  

# revision 0.06 . KD 6/2010 - add UNTRUSTED_RECTANGLE and UNTRUSTED_ELLIPSE; use `whereis catmar` and so on  

# revision 0.07 . KD 6/2010 - decide about ORGX/Y info in MAR header being pixels or mm; other fixes      

# revision 0.08 . KD 6/2010 - fixes for Pilatus 6M                                                        

# revision 0.09 . KD 6/2010 - get rid of requirement for mccd_xdsparams.pl and/or catmar; rather use "od"  

# revision 0.10 . Tim Gruene 7/2010 - set link 'images' to image directory if path exceeds 72 characters  

# revision 0.11 . KD 7/2010 - for MarCCD: look for distance info at different byte position                

# revision 0.12 . KD 7/2010 - fix for negative PHISTART                                                    

# revision 0.13 . KD 8/2010 - store correct NX NY QX QY in XDS.INP                                        

# revision 0.14 . KD 1/2011 - SENSOR_THICKNESS for Pilatus; MINIMUM_NUMBER_OF_PIXELS_IN_A_SPOT=3          

# revision 0.16 . KD 3/2011 - SENSOR_THICKNESS=0.01 for ADSC and MarCCD. Add comment about SILICON=        

# revision 0.17 . KD 3/2011 - make it work for .bz2 frames; improve screen output                          

# revision 0.18 . KD 4/2011 - faster by doing "strings" only once; revert "images/${1##/*/}" "correction"  

# revision 0.19 . KD 6/2011 - bugfix for 0.18                                                              

# revision 0.20 . KD 7/2011 - redirect stderr of /bin/ls to /dev/null                                      

# revision 0.48 . Kip Guja 11/2014 add detector serial number for ALS 5.0.2 to beam center convention 1  

REVISION="0.77 (19-Dec-2017)"

#                                                                                                          

# usage: e.g. generate_XDS.INP "/file/system/frms/mydata_1_???.img"                                                      

# make sure to have the two quotation marks !                                                              

# the ? are wildcards for the frame numbers.                                                                

#                                                                                                          

# known problems:                                                                                          

# - the same might be a problem for MAR headers, too (not sure about this)  

# - on Mac OS X, Xcode (from http://developer.apple.com/tools/xcode) might be needed - not sure about this                             

#                                                                                                              

# notes for debugging of the script:                                                                            

# - add the -v option to the first line, to see where an error occurs                                          

# - comment out the removal of tmp1 and tmp2 in the last line                                                  

#                                                                                                              

# ====== Start of script ======                                                                                

echo generate_XDS.INP version $REVISION . Obtain the latest version from                                        

echo http://strucbio.biologie.uni-konstanz.de/xdswiki/index.php/generate_XDS.INP                                

if [ "$1" == "help" ] || [ "$1" == "-help" ] || [ "$1" == "-h" ]; then                                          

   echo usage: generate_XDS.INP \"/file/system/frms/mydata_1_???.img\"  \(_with_ the quotation marks!\)                      

   echo if the frames are compressed with bzip2, leave out the .bz2 extension!                                  

   exit                                                                                                          

fi                                                                                                              

#                                                                                                              

# defaults:                                                                                                    

#                                                                                                              

DETECTOR="XXX MINIMUM_VALID_PIXEL_VALUE=XXX OVERLOAD=XXX"                                                      

ORGX=XXX                                                                                                        

ORGY=XXX                                                                                                        

DETECTOR_DISTANCE=XXX                                                                                          

OSCILLATION_RANGE=XXX                                                                                          

X_RAY_WAVELENGTH=XXX                                                                                            

QX=XXX                                                                                                          

QY=XXX                                                                                                          

NX=XXX                                                                                                          

NY=XXX                                                                                                          

SENSOR_THICKNESS=0  

TRUSTED_REGION="0.0 1.2 ! partially use corners of detector (0 1.4143: use all pixels)"                                                                                          

MNOPIAS=3  

dname=`echo "$1" | xargs dirname`

bname=`echo "$1" | xargs basename`

# see how we are called:                                                                                        

# don't accept the "direct beam" shot at SLS/Pilatus PX-I and PX-II  

find -H $dname -name "$bname" -or -name "${bname}.bz2" -or -name "${bname}.gz" -or -name "${bname}.xz" | egrep -v "_00000.cbf|_000.img" | sort > tmp1

   

  data_half=`echo "scale=0; $data_num/2" | bc -l`                      

DET=XXX            

IFS=$'\n'  

  strings $FIRSTFRAME > tmp2  

  # the reason is that FIRSTFRAME may contain a blank, which makes some commands fail                                              

  egrep -qi 'marccd|Corrected' tmp2 && DET=mccd                                

  grep -q PILATUS tmp2            && DET=pilatus                            

  grep -q BEAM_CENTER_X tmp2      && DET=adsc                                

  grep -q mar345 tmp2 && DET=MAR345    

  grep -q BRUKER tmp2 && grep -q CBF tmp2  && DET=Bruker-cbf  

  grep -q CMOS1 tmp2 && DET=adsc-CMOS1                  

# identify other detector types in the same way  

   echo "this is not a MAR, ADSC/SMV or PILATUS detector - fill in XXX values manually!"

   DETECTOR="XXX MINIMUM_VALID_PIXEL_VALUE=XXX OVERLOAD=XXX"                          

elif [ "$DET" == "mccd" ]; then  

   SENSOR_THICKNESS=0.01                                          

   # use first frame of dataset to obtain parameters              

   let SKIP=1024+80                                                      

   let SKIP=$SKIP+4                                                                        

   DETECTOR_DISTANCE=`echo "scale=3; $DETECTOR_DISTANCE/1000" | bc -l`                                    

   let SKIP=1024+256+128+256+4                                                                            

   ORGX=$(od -t dI -j $SKIP -N 4 tmp2 | head -1 | awk '{print $2}')            

   ORGX=`echo "scale=2; $ORGX/1000" | bc -l `                                                            

   let SKIP=$SKIP+4                                                                                      

   ORGY=$(od -t dI -j $SKIP -N 4 tmp2 | head -1 | awk '{print $2}')            

   ORGY=`echo "scale=2; $ORGY/1000" | bc -l `                                                            

# fixed Aug 30, 2017 after IUCr2017 @ Hyderabad  

   OSCILLATION_RANGE=$(od -t dI -j $SKIP -N 4 tmp2 | head -1 | awk '{print $2}')  

   OSCILLATION_RANGE=`echo "scale=3; $OSCILLATION_RANGE/1000" | bc -l`  

   let SKIP=1024+256+128+256+128+4                                                              

   QX=$(od -t dI -j $SKIP -N 4 tmp2 | head -1 | awk '{print $2}')    

   QX=`echo "scale=10; $QX/1000000" |bc -l `                                                    

   let SKIP=$SKIP+4                                                                              

   QY=$(od -t dI -j $SKIP -N 4 tmp2 | head -1 | awk '{print $2}')    

   QY=`echo "scale=10; $QY/1000000" |bc -l `                                                    

   X_RAY_WAVELENGTH=`echo "scale=5; $X_RAY_WAVELENGTH/100000" | bc -l`                                  

   NXBYFOUR=`echo "scale=0; $NX/4" | bc -l `                              

   ORGXINT=`echo "scale=0; $ORGX/1" | bc -l `                            

   if [ $ORGXINT -lt $NXBYFOUR ]; then                                    

     ORGX=`echo "scale=1; $ORGX/$QX" | bc -l`                            

     ORGY=`echo "scale=1; $ORGY/$QY" | bc -l`                            

   else                                                                              

     echo MARCCD detector: header ORGX, ORGY seem to be in pixel units              

   fi                                                                                

   SENSOR_THICKNESS=0.01                                                      

       NY=`grep SIZE2 tmp2 | tail -1 | sed s/SIZE2=//`                                            

       QY=$QX                                              

# at ESRF, PF, ALS 5.0.2 and ... (pls fill in!) the following should be used:        

       echo - at ESRF, PF, APS Argonne BLs use: ORGX=$ORGX1 ORGY=$ORGY1                  

! ORGX= $ORGX1 ORGY= $ORGY1 ! For ESRF, PF, APS Argonne ..."

       echo - at e.g. ALS 5.0.3 use: ORGX=$ORGX2 ORGY=$ORGY2                                                        

# this 3rd alternative convention should be used at the following beamlines (pls complete the list): ALS 8.2.2, ...  

       echo - at e.g. SPring-8 use: ORGX=$ORGX4 ORGY=$ORGY4  

       # For convention 1; Known PF detectors = 449: NW12A Q210, 472: NE3A Q270, 474: BL17A Q270, 912: BL5A Q315, 923: ALS BL5.0.2 Q315, 933: AichiSR BL2S1 Q315, 916: APS 24 IDE

       # find DETECTOR_DISTANCE and OSCILLATION_RANGE:                                                                              

       DETECTOR_DISTANCE=`grep ^DISTANCE tmp2 | sed s/DISTANCE=//`                                                                    

       OSCILLATION_RANGE=`grep OSC_RANGE tmp2 | sed s/OSC_RANGE=//`      

   SENSOR_THICKNESS=0.01                                                      

   NY=`grep SIZE2 tmp2 | tail -1 | sed s/SIZE2=//`                                            

   QY=$QX                                              

   DETECTOR_DISTANCE=`grep ^DISTANCE tmp2 | sed s/DISTANCE=/-/`                                                                    

   OSCILLATION_RANGE=`grep OMEGA_DELTA tmp2 | sed s/OMEGA_DELTA=//`   

   DIRECTION_OF_DETECTOR_X_AXIS="-1 0 0"                                                                

   DETECTOR="PILATUS MINIMUM_VALID_PIXEL_VALUE=0 OVERLOAD= 1048576 !PILATUS"

   QX=0.172 QY=0.172                                                        

   grep Pixel_size tmp2 && QX=`awk '/Pixel_size/{print 1000*$2}' tmp2`

   grep -q Pixel_size tmp2 && QY=`awk '/Pixel_size/{print 1000*$5}' tmp2`                                                    

       SENSOR_THICKNESS=`grep thickness tmp2 | sed -e s/'Silicon sensor, thickness'// | awk '{print $1*1000}'`

       # find X_RAY_WAVELENGTH:                                                                              

       X_RAY_WAVELENGTH=`grep Wavelength tmp2 | sed -e s/Wavelength// -e s/A// | awk '{print $1}'`          

       NY=`grep X-Binary-Size-Second-Dimension tmp2 | awk '{print $2}'`  

       # find DETECTOR_DISTANCE and OSCILLATION_RANGE:

       DETECTOR_DISTANCE=`awk '/distance/{print $2}' tmp2`

       OSCILLATION_RANGE=`awk '/Angle/{print $2}' tmp2`

        

           rotation_axis="0 -1 0"

# insert similar code for Petra P14 here

   OVERLOAD=`h5dump -d "/entry/instrument/detector/detectorSpecific/countrate_correction_count_cutoff" $FIRSTFRAME | awk '/\(0\):/{print $2}' `

   DETECTOR="EIGER MINIMUM_VALID_PIXEL_VALUE=0 OVERLOAD= $OVERLOAD !EIGER"

   echo Data from a Eiger hdf5

      SENSOR_THICKNESS=`h5dump -d "/entry/instrument/detector/sensor_thickness" $FIRSTFRAME | awk '/\(0\): [0-9]/{print $2*1000}'`

      X_RAY_WAVELENGTH=`h5dump -d "/entry/instrument/beam/incident_wavelength" $FIRSTFRAME | awk '/\(0\): [0-9]/{print $2}'`

      NX=`h5dump -d "/entry/instrument/detector/detectorSpecific/x_pixels_in_detector" $FIRSTFRAME | awk '/\(0\): [0-9]/{print $2}'`

      NY=`h5dump -d "/entry/instrument/detector/detectorSpecific/y_pixels_in_detector" $FIRSTFRAME | awk '/\(0\): [0-9]/{print $2}'`

      # find ORGX and ORGY:

      ORGX=`h5dump -d "/entry/instrument/detector/beam_center_x" $FIRSTFRAME | awk '/\(0\): [0-9]/{print $2}'`

      ORGY=`h5dump -d "/entry/instrument/detector/beam_center_y" $FIRSTFRAME | awk '/\(0\): [0-9]/{print $2}'`

      DETECTOR_DISTANCE=`h5dump -d "/entry/instrument/detector/detector_distance" $FIRSTFRAME | awk '/\(0\): [0-9]/{print $2*1000}'`

      OSCILLATION_RANGE=`h5dump -d "/entry/sample/goniometer/omega_range_average" $FIRSTFRAME | awk '/\(0\): [0-9]/{print $2}'`

   echo Data from a RAXIS detector

   NX=$(python -c 'import struct; f=open("'$FIRSTFRAME'","rb");f.seek(768);print "%.4d"%struct.unpack(">i",f.read(4))')

   NY=$(python -c 'import struct; f=open("'$FIRSTFRAME'","rb");f.seek(772);print "%.4d"%struct.unpack(">i",f.read(4))')

 

   DETECTOR_DISTANCE=$(python -c 'import struct; f=open("'$FIRSTFRAME'","rb");f.seek(344);print "-%.4f"%struct.unpack(">f",f.read(4))')

 

   ORGX=$(python -c 'import struct; f=open("'$FIRSTFRAME'","rb");f.seek(540);print "%.4f"%struct.unpack(">f",f.read(4))')

   ORGY=$(python -c 'import struct; f=open("'$FIRSTFRAME'","rb");f.seek(544);print "%.4f"%struct.unpack(">f",f.read(4))')

 

   OSCILLATION_RANGE=$(python -c 'import struct; f=open("'$FIRSTFRAME'","rb");f.seek(524);phis,phie=struct.unpack(">ff",f.read(8));print "%.4f"%(phie-phis)')

 

   QX=$(python -c 'import struct; f=open("'$FIRSTFRAME'","rb");f.seek(776);print "%.6f"%struct.unpack(">f",f.read(4))')

   QY=$(python -c 'import struct; f=open("'$FIRSTFRAME'","rb");f.seek(780);print "%.6f"%struct.unpack(">f",f.read(4))')

 

   X_RAY_WAVELENGTH=$(python -c 'import struct; f=open("'$FIRSTFRAME'","rb");f.seek(292);print "%.6f"%struct.unpack(">f",f.read(4))')

 

     # Find rotation axis. Warning: currently not support flipping (det(tmpmat)<0)

     # I'm not sure this method is really valid - but at least mosflm seems to read this to determine rotation axis.

    tmpmat=(`grep ${dname}SPATIAL_DISTORTION_VECTORS= tmp2 | tail -1 | sed -e "s/.*=//"`)

    rotx=`echo "scale=6; -1.0*${tmpmat[0]}" | bc -l`

  if [ `echo "(${tmpmat[0]}*${tmpmat[3]}-(${tmpmat[1]}*${tmpmat[2]}))/1"|bc` -lt 0 ]; then

     echo "WARNING!! not-supported SPATIAL_DISTORTION_VECTORS header detected."

   else # should be RX_

    DETECTOR="RAXIS MINIMUM_VALID_PIXEL_VALUE=0"

     dtrek_det="raxis"

 

   # find DETECTOR_DISTANCE  

    DETECTOR_DISTANCE="-${GONIO_VALUES[$idx]}"

    # TODO: Check unit!! - ${GONIO_UNITS[$idx]}

   # find OSCILLATION_RANGE

  NX=`awk '/FORMAT/{print $2}' tmp2`

  NY=$NX

  QX=`awk '/PIXEL/{print $3/1000.}' tmp2`

  QY=`awk '/PIXEL/{print $5/1000.}' tmp2`

  DETECTOR_DISTANCE=`grep 'DISTANCE' tmp2 | awk '{print $2}'`              

  X_RAY_WAVELENGTH=`grep 'WAVELENGTH' tmp2 | awk '{print $2}'`  

  OSCILLATION_RANGE=`grep 'PHI' tmp2 | awk '{print $5-$3}'`  

  TRUSTED_REGION="0 0.99"      

  MNOPIAS=6  

  TRUSTED_REGION="0 1.42"      

  pol_frac=`awk '/polarizn_source_ratio/{print $2}' tmp2`

  OVERLOAD=`awk '/_array_intensities.overload/{print $2}' tmp2`

  NX=`awk '/X-Binary-Size-Fastest-Dimension/{print $2}' tmp2`

  NY=`awk '/X-Binary-Size-Second-Dimension/{print $2}' tmp2`

  QX=`awk '/ELEMENT_X ELEMENT_X/{print $4}' tmp2`

  QY=`awk '/ELEMENT_Y ELEMENT_Y/{print $4}' tmp2`

  DETECTOR_DISTANCE=`grep "0 0 ? ? ?" tmp2 | awk '/DX /{print $2}'`  

     rotation_axis="$R1 $R2 $R3"

# here we could check if DELTAPHI is <0, and if so, negate it and rotation_axis

     rotation_axis="0 -1 0"

# here we could check if DELTAOMEGA is <0, and if so, negate it and rotation_axis

   TWOTHETA=`awk '/TWOTHETA \? \? \?/{print $5}' tmp2`

   exit 1                    

fi                          

echo ORGX= $ORGX ORGY= $ORGY - check these values with adxv !

echo DETECTOR_DISTANCE= $DETECTOR_DISTANCE  ! only read by XYCORR, IDXREF                    

echo OSCILLATION_RANGE= $OSCILLATION_RANGE  ! only read by IDXREF                

echo X-RAY_WAVELENGTH= $X_RAY_WAVELENGTH    ! only read by IDXREF                  

# now we know everything that is required to generate XDS.INP

ORGX= $ORGX ORGY= $ORGY  ! values from frame header; only read by XYCORR, IDXREF    

DETECTOR_DISTANCE= $DETECTOR_DISTANCE !read by XYCORR, IDXREF. Negative if detector normal points to crystal.              

OSCILLATION_RANGE= $OSCILLATION_RANGE                  

X-RAY_WAVELENGTH= $X_RAY_WAVELENGTH                    

! REFERENCE_DATA_SET=xxx/XDS_ASCII.HKL ! e.g. to ensure consistent indexing

DATA_RANGE=$DATA_RANGE                                                    

SPOT_RANGE=$SPOT_RANGE                                                    

! BACKGROUND_RANGE=1 10 ! rather use defaults (first 5 degree of rotation)  

! close spots/long cell axis: reduce SEPMIN and CLUSTER_RADIUS from their defaults of 6 and 3

! SEPMIN=4 CLUSTER_RADIUS=2 ! should be default for Pilatus and other detectors with low PSF

! since XDS 01-MAR-2015, POSITION supersedes DISTANCE. DISTANCE only for compatibility w/ old XDS.

! nowadays headers are usually correct so refine DISTANCE/POSITION in INTEGRATE but not IDXREF

REFINE(IDXREF)=CELL BEAM ORIENTATION AXIS POSITION ! works best with XDS >= BUILT=20170720

REFINE(INTEGRATE)= DISTANCE POSITION BEAM ORIENTATION ! AXIS CELL

! REFINE(CORRECT)=CELL BEAM ORIENTATION AXIS DISTANCE POSITION ! Default is: refine everything

 

! Better: read the article http://strucbio.biologie.uni-konstanz.de/xdswiki/index.php/SILICON .

 

if [ "$DET" == "raxis" -o "$dtrek_det" == "raxis" ]; then

DIRECTION_OF_DETECTOR_X-AXIS=1 0 0

DIRECTION_OF_DETECTOR_Y-AXIS=0 -1 0

!FRACTION_OF_POLARIZATION=0.98   ! uncomment if synchrotron; only used by CORRECT

POLARIZATION_PLANE_NORMAL=1 0 0  ! only used by CORRECT

  echo 'ROTATION_AXIS=1 0 0  ! Australian Synchrotron, SERCAT ID-22 (?), APS 19-ID (?), ESRF BM30A, SPring-8, SSRF BL17U need -1 0 0. Diamond ID24 needs 0 -1 0' >> XDS.INP

  echo 'ROTATION_AXIS=-1 0 0  ! if this is wrong, please contact author.' >> XDS.INP

DIRECTION_OF_DETECTOR_X-AXIS=-1 0 0

DIRECTION_OF_DETECTOR_Y-AXIS= 0 1 0

DIRECTION_OF_DETECTOR_X-AXIS=$DIRECTION_OF_DETECTOR_X_AXIS

DIRECTION_OF_DETECTOR_Y-AXIS=0 1 0

  if [ -e /usr/local/lib64/dectris-neggia.so ]; then

    echo LIB=/usr/local/lib64/dectris-neggia.so >> XDS.INP

    echo !LIB=/usr/local/lib64/dectris-neggia.so >> XDS.INP

    echo /usr/local/lib64/dectris-neggia.so was not found!

NUMBER_OF_PROFILE_GRID_POINTS_ALONG_ALPHA/BETA=13 ! Default is 9 - Increasing may improve data  

NUMBER_OF_PROFILE_GRID_POINTS_ALONG_GAMMA=13      ! accuracy, particularly if finely-sliced on phi,  

!                                                  and does not seem to have any downsides.  

   if [ $NX == "1475" ]; then

! flat_field correction; they cut away one additional pixel adjacent  

   elif [ $NX == "2463" ]; then

   elif [ $NX == "3110"  -a $NY == "3269" ]; then

   elif [ $NX == "4150" -a $NY == "4371" ]; then

  UNTRUSTED_RECTANGLE= 2069 2082      0 4371  

echo Full documentation, including complete detector templates, is at

echo http://www.mpimf-heidelberg.mpg.de/~kabsch/xds . More documentation in XDSwiki

echo After running xds, inspect, using XDS-Viewer, at least the beamstop mask in

echo BKGPIX.cbf, and the agreement of predicted and observed spots in FRAME.cbf!