#!/bin/csh -f

set verify
set echo

# (c)2011 Rigaku Americas Corporation
# This script was written by 
#
# JW Pflugrath
# Rigaku Americas Corporation
#
# as part of the ACA 2011 session "Data Processing with the Pros"
#
# This script shows how to process the 2QVO diffraction images from
# the APS/ID-?? which has PDB entry 2QVO
#
# Images are found in the directory $IMAGE_DIR
#

setenv IMAGE_DIR /ac1/jwp/2QVO/2qvo

# There are a few scans there, only some of which are treated in this script

#a Process scan a.
#a Images are ACA10_AF1382_1.0001 to ACA10_AF1382_1.360
#a (a second scan has images are ACA10_AF1382_2.0001 to ACA10_AF1382_2.360)
#a From the image header, they appear to be 1 degrees per image.
#a The beam center in the image header is close enough.
#a

#a One can use dtdisplay/dtprocess (the GUI set), but one can also script 
#a all this, which is what this script is about.  So start with the
#a sample_script.scom found in ${DTREK_ROOT}/bin and modify it ...
#a ... or use the scripts created by dtprocess
#a
#a OTHER THINGS to CHECK
#a Is the beam center? -> Yes
#a Is the distance correct? -> Yes
#a Is the wavelength correct? -> Yes
#a Is the direction of the rotation axis correct-> NO!
#a   This was determined by predicting the spots on images 2 through 10
#a   AFTER indexing from spots on image 1.  Image 1 matches, but 2-10 do not.
#a   This means the direction of the rotation axis is not correctly specified
#a   to d*TREK.  
#a Steps to take:
#a   1. View image header with dtprocess. See that SCAN_ROTATION_AXIS_NAME is
#a      Phi.
#a   2. Change the direction of the Phi axis by
#a   3. See that the Phi axis is the THIRD axis specified in 
#a       CRYSTAL_GONIO_NAMES
#a   4. Change the vector specifying the Phi axis by multiplying by -1.
#a      Keyword is CRYSTAL_GONIO_VECTORS, last 3 numbers must be multiplied
#a      by -1.  We can use the d*TREK override mechanism by specifying an
#a      environment variable with the same name as the keyword, thus:
setenv CRYSTAL_GONIO_VECTORS "-1 0 0 0 0 1 1 0 0"
#a The above environment variable will be used INSTEAD of the header
#a keyword in any d*TREK module where the env variable exists,
#a such as THIS SCRIPT.  Be careful.  If this environment variable remains
#a defined globally, it will override other d*TREK sessions.

#a 
###############
###############
###############

# This sample script demonstrates one of many possible ways to use d*TREK 
# without the graphical user interface. Note that comment lines begin with 
# a hash or number sign #
# Examine each line that begins with setenv and possibly edit it (I recommend
# using the emacs or xemacs editor), save it, then enter 'sample_script.com' 
# on the command line. This script performs for one scan of protein images:
#    Find spots
#    Index spots
#    Refine spots and reflections
#    Predict reflections
#    Integrate reflections
#    Scale and average reflection

# Not all the options below need to be used.  For example, the -spacegroup
# option for the dtindex step is usually not needed.  However, one should 
# understand what the options of the dtintegrate step are doing.  Please read
# the overview*.pdf for more help on parameters.

# All output file names will begin with the following prefix.  These output
# files are the log files, the .head files, and the .ref files.
# Use a different prefix to process a different way or different images
# in the same directory.
#
setenv DTREK_PREFIX     a_

# Each scan (0-90, 180-270, different wavelengths) should have a unique batch
# prefix.
#
setenv BATCH_PREFIX     1

# The image name is used to derive a templace for all the images.
# Use the first image in a scan, if possible
#
setenv IMAGE_NAME       ${IMAGE_DIR}/ACA10_AF1382_1.0001

# The starting header can come from the first image or another .head file.
# It should have the correct direct beam position, detector distance, detector
# swing angle, wavelength, and other data collection parameters.  If these
# are already correct in the first image of the scan, then just use the first
# image. Otherwise you may wish to use dtdisplay and/or dtprocess to create 
# a .head file with the correct values and specify that file here.
#
setenv START_HEAD       ${IMAGE_DIR}/ACA10_AF1382_1.0001

# Specify the image sequence to use for finding spots for indexing.  Usually
# the first image is OK, the first two are often better.
#
setenv IMAGE_SEQ_FIND   "1 1"

# Specify the image sequence to use for integration.  Usually all images in the
# scan is what you want.
#
setenv IMAGE_SEQ_INT    "1 360"


# Create a beam.mask file to mark beamstop shadows and other bad pixels 
# with dtdisplay/Edit/Erase (edit mask) mode, # then File/Save as... beam.mask

#a
#a With this scan, the beam stop shadow is visible, so either manually
#a create an image with the bad pixels masked or use dtbeammask with 
#a known areas
#a For this MAR CCD, you may wish to mask the seams as well.
#a The seams create a tic-tac-toe game
#a SPECIAL NOTE: This detector has 2 edges with WIDE inactive areas

dtbeammask ${IMAGE_NAME} -circle 1980 2040 70 \
           -rect 0 0 4096 106 -rect 0 0 150 4096 \
           a_beam.mask
setenv MARCCD_NONUNF_INFO a_beam.mask
setenv MARCCD_NONUNF_TYPE  Simple_mask

#setenv CCD_NONUNF_INFO  beam.mask
#setenv CCD_NONUNF_TYPE  Simple_mask
##setenv RX_NONUNF_INFO  None            # Use these for no mask at all.
##setenv RX_NONUNF_TYPE  None            # (R-AXIS only)
#setenv RX_NONUNF_INFO  beam.mask
#setenv RX_NONUNF_TYPE  Simple_mask


# If desired, insert -reso commands at the integrate and scaling steps.
# The default is to process all reflections that impinge on the detector.
#
setenv DT_INTEGRESO    ""                 # "-reso 100 2.6"
setenv DT_SCALERESO    ""                 # "-reso 100 2.5"

# If you wish to specify a spacegroup, then do so in the next two lines.
setenv DT_SPACEGROUP   ""                  # if you DO NOT know spacegroup
##setenv DT_SPACEGROUP   "-spacegroup 75"       # if you know spacegroup

# Sometimes after a run of this script you wish to rescale with a different
# resolution or spacegroup.  Change resolution and spacegroup above, then
# run the script with 'sample_script.com SCALE' in order to jump ahead
# to the dtscaleaverage step.
#
if (SCALE == $1) then
  goto SCALE_ONLY
endif

##### Generally, you do not need to change anything below this line. ####
# Merge info from the image header into our starting .head file
#
dtheadermerge ${START_HEAD} ${IMAGE_NAME} +SCAN\* \
              +ROTATION\* +\*_GONIO_VALUES +SOURCE_WAVELENGTH \
      ${DTREK_PREFIX}newmerge.head >& ${DTREK_PREFIX}dtheadermerge.log

# Find spots in some of the images
#
dtfind ${DTREK_PREFIX}newmerge.head -seq ${IMAGE_SEQ_FIND} -sigma 3 \
  -min 50 -filter 6 -window 0 0 -display -out \
  ${DTREK_PREFIX}dtfind.head >& ${DTREK_PREFIX}dtfind.log

# Autoindex the spots that were found.  We know the spacegroup.
#
dtindex ${DTREK_PREFIX}dtfind.head ${DTREK_PREFIX}dtfind.ref -beamcheck 20 18 \
       -maxresid 3.0 -sigma 5 ${DT_SPACEGROUP} >& ${DTREK_PREFIX}dtindex.log

# Refine the autoindexing solution with the dtfind spot list
#
dtrefine ${DTREK_PREFIX}dtindex.head ${DTREK_PREFIX}dtfind.ref \
    +CrysAll +DetAll -sigma 5.0 -rej 0.5 0.5 0.5 -cycles 500 -display \
	 -go -go -go >& ${DTREK_PREFIX}dtrefine1.log

# Refine the autoindexing solution from images
#
dtrefine ${DTREK_PREFIX}dtrefine.head -seq ${IMAGE_SEQ_FIND} \
    +CrysAll +DetAll -sigma 5.0 -rej 0.5 0.5 0.5 -cycles 500 -display \
	 -go -go >& ${DTREK_PREFIX}dtrefine2.log

# Predict to set mosaicity where we want it to start at
#
dtpredict ${DTREK_PREFIX}dtrefine.head -seq ${IMAGE_SEQ_FIND} \
    -mosaicity 1.5000 -display >& ${DTREK_PREFIX}dtpredict.log

# Integrate the images.  Note: One should edit the options
#
dtintegrate ${DTREK_PREFIX}dtpredict.head -window 0 0 -pad 1 \
   ${DT_INTEGRESO} -mosaicitymodel 1.000 0.000 -profit 50 7 -batch 1 4 \
   -batchprefix ${BATCH_PREFIX} -display -prerefine 2 -seq ${IMAGE_SEQ_INT} \
    >& ${DTREK_PREFIX}dtintegrate.log

# Scale and average the integrated and profile-fitted reflns
#
SCALE_ONLY:

dtscaleaverage ${DTREK_PREFIX}dtintegrate.head ${DTREK_PREFIX}dtprofit.ref \
  ${DT_SCALERESO} -sigma 5.0 -anom -errormodel  -reject .0075 -batchscale \
  -reqab spherical 4 3 \
 ${DTREK_PREFIX}dtscale.ref >& ${DTREK_PREFIX}dtscaleaverage.log

rm ${DTREK_PREFIX}*tmp*
rm ${DTREK_PREFIX}*_?.head
rm ${DTREK_PREFIX}*_??.head

#a BOTTOM LINE:
#a   The "gotchas" in this processing are
#a   1. Phi axis rotates the other way
#a   2. Detector has inactive edges
#a   3. Determine the spacegroup:
#a         dtcell a_dtprofit.ref 
#a      shows the spacegroup to be P4(2)
#a   4. Re-run dtscaleaverage in this spacegroup
#a   5. Plot the results with jdtplot and observe that scale factors for
#a      image 105, 175, and 355 do not make sense, so reject these images
#a   6. We expect an anomalous signal, so make Reduced ChiSq A close to 1
#a      by using the -chisq 1.2 option.  Otherwise sigmas will be too high.
#a   7. Cut off resolution where <I/sigI> of the OBSERVATIONS is 1.5 and 
#a      completeness is above 95%.  That is: 2.4 Ang

#a Here is the final scaling script for scan 1:
dtscaleaverage a_dtintegrate.head a_dtprofit.ref -spacegroup 77 -reject sigma 5.0 -batchrestrain 0.002 -reso 999.0 2.4 -chisq  1.2000 -rejectbatch id 10105,10175,10355 -anom -ref a_dtunavg.ref -errormodel  -reject fraction .0075 -batchscale -reqab spherical 4 3 a_dtscale.ref >& a2_dtscaleaverage.log

#a Plot the log file:
# jdtplot b2_dtscaleaverage.log

