#!/usr/bin/perl
# Maxim Markevitch <maxim % head-cfa.harvard.edu>, last change 11 Jun 2010
#
# This simple script creates a model background event file for approximate
# correction of the ACIS readout artifact (a.k.a. out-of-time events or
# trailed images, see POG Sec. 6.12.4), under the assumption of no
# pileup. Its intended use is for spectral analysis of the outskirts of
# galaxy clusters with sharp central brightness peaks. (For piled-up point
# sources, see the CIAO tool acisreadcorr.) 
# 
# The script takes an input level-1 event file with your source, randomizes
# CHIPY values and reruns acis_process_events to recalculate the sky coords
# and PI values, applying the CTI correction and the gain map appropriate
# for the new CHIPY values. Note that this is not physically correct,
# because the true gain and CTI values for an artifact event correspond to
# the source position rather than the recorded, shifted event position. But
# the PI values of the artifact events in the data are similarly calculated
# incorrectly from their recorded coordinates (because acis_process_events
# has no way of distinguishing artifact events from regular events), so we
# are right to do it the same way for our model background.
# 
# Spectra and images extracted from this event file should be normalized by
# factor 0.0128 or 0.0132 for the 6-CCD or 5-CCD ACIS mode, respectively
# (=41ms / 3.2s or 3.1s, where 41ms is the readout time, 3.2s or 3.1s is the
# NO-SUBARRAY frame time for 6 or 5 CCD modes). To this end, the script
# optionally appends a fake GTI extension containing the observation's clean
# exposure (supplied as an argument) multiuplied by 78.05 (=3.2s/41ms) so
# that most standard tools could automatically calculate the right
# normalization when extracting spectra from this file. This assumes a 6-CCD
# mode; MULTIPLY THE EXPOSURE ARGUMENT BY 0.9688 (3.1s/3.2s) FOR 5-CCD
# OBSERVATIONS. Note that to use the fake GTI option, the input event file
# should already be cleaned of any flare intervals, etc., so that no
# additional filtering will be required when extracting spectra or images.
# The existing GTI extensions in the input file are not copied to the output
# file.
# 
# Optionally, the script also reapplies the bad pixel filter.
# 
# Check the structure of the resulting file (GTI extensions, etc.) as I
# may not have kept up with changes of the acis_process_events behavior.
#
# Command-line arguments (in this order):
#	1. input event filename
#	2. output filename
#	3. aspect filename (*aoff1* or *asol1*)
#	4. bad pixel ASCII table ("none" to skip bad pixel re-filtering)
#	5. clean exposure (see above; "0" to skip the fake GTI step)
#	6. gain file (full filename or "default" = current CALDB)
#	7. CTI file (full filename, "default" or "none")
#
# For CTI correction of randomized events, the input event file should
# contain the necessary columns: PHAS, CCD_ID, CHIPX, CHIPY. To skip it (not
# recommended), set the CTI file argument to "none".
# 
# The script calls various ftools, acis_process_events and (optionally)
# A.Vikhlinin's badpixfilter (/home/maxim/BIN/[SOL,LIN]/badpixfilter at CfA).
#
##############################################################################
# CHANGES:
# 6/10/10: replaced default gain and cti files with CALDB, added instruction
#          for 5-ccd mode normalization
# 2/11/5: changed default gainfile
# 8/5/03: added CTI correction; added "; *" to the 'fcopy' column editing 
#         syntax to keep up with the FTOOLS improvement (as of FITSIO 2.4)
#         -- uncomment the old 'fcopy' line below if needed.
# 7/22/01: now do not copy existing GTI exts. from input to output; also,
#          append only one GTI ext. (used to copy 2, GTI and STDGTI).
# 3/23/01: for acis_process_events, give aoff1 file as the alignmentfile
#          argument.
# 12/7/00: 1. in ciao-2, acis_process_events interface changed, so now call all
# 	   parameters explicitly (par=value) to be compatible with old ciao
#          2. added gainfile as optional 6th argument.
# 4/27/00: added PI recalculation for the randomized coords.
# 3/29/00: original version
##############################################################################


####### EDIT THESE FILENAMES IF NOT AT CFA: #######
#
$defaultgain="CALDB";
$defaultcti="CALDB";

$badpixfilter="badpixfilter"; 
# (make sure it's on the path)

# Also, depending on version of FTOOLS, may need to edit the line with 
# fcopy below.
#
###################################################


if (@ARGV eq 7) {
  ($in, $out, $aspect, $badpix, $expos, $gain, $cti) = @ARGV;

  if (($gain eq "default") || ($gain eq "DEFAULT")) {
    $gain=$defaultgain;
    print "Using default gain file:\n";
    print "$gain\n";
  }

  if (($cti eq "none") || ($cti eq "NONE")) {
    $applycti="no";
    $regrade="no";
    print "Skipping CTI correction\n";
  }
  else {
    $applycti="yes";
    $regrade="yes";
  }

  if (($cti eq "default") || ($cti eq "DEFAULT")) {
    $cti=$defaultcti;
    print "Using default CTI file:\n";
    print "$cti\n";
  }
}
else {
  die "Usage: \nmake_readout_bg in.fits out.fits aspectfile badpixtable expos gainfile ctifile\nMore instructions inside script.\n";
}


## Randomize chipy:

print "Randomizing chipy...\n";

unlink "TMPJNK", $out;
`fextract $in"[events]" TMPJNK`; ## (to get rid of all existing gti exts.)

#######
## old (ver<2.4) fitsio, uncomment if the next line doesn't work:
#`fcopy TMPJNK"[events][col chipy=(int)(2+1022*random())]" $out`;

## new (ver>2.4) fitsio:
`fcopy TMPJNK"[events][col chipy=(int)(2+1022*random()); *]" $out`;
#######

unlink "TMPJNK";


## Now recalculate sky coords and PI. Have to set doevtgrade=yes for CTI
## correction even though the effect on grades is negligible for our current
## purposes (and would need grade-unfiltered evt file to have any use for
## the new grades):

print "Running acis_process_events to recalculate sky coords and PI values...\n";

`acis_process_events infile=$out outfile=TMPJNK acaofffile=$aspect gainfile=$gain ctifile=$cti apply_cti=$applycti doevtgrade=$regrade calculate_pi=yes eventdef="{s:ccd_id,s:node_id,s:chip,l:pha,f:energy,s:pi,s:fltgrade,s:grade,f:det,f:sky,d:time}"`;

print "\n";
rename "TMPJNK", $out;


## Reapply bad pixel map, if given:

if (($badpix ne "NONE") && ($badpix ne "none")) {
  print "Reapplying bad pixel filter...\n";

  `$badpixfilter $out TMPJNK $badpix`;
  rename "TMPJNK", $out;
}


## Append a fake GTI extension containing an interval of length equal to
## exposure*3.2/0.041 and starting at TSTART:

if ($expos > 0) {
  print "Appending a fake GTI extension -- see instructions inside this script\n  on 6-CCD and 5-CCD ACIS modes.\n";

  chomp($line = `fkeyprint $in"[0]" TSTART | grep "TSTART  ="`);
  ($j,$j, $tstart, $j) = split(' ',$line,4);

  # this is for 6-ccd modes; 5-ccd mode ontime should be expos*3.1/0.041
  $ontime=$expos*3.2/0.041;
  $tstop=$tstart+$ontime;

  ## create template files for creating gti fits file:

  unlink 'TMPJNKtempl', 'TMPJNKdata';
  `echo "start 1d\nstop 1d" > TMPJNKtempl`;
  `echo $tstart $tstop > TMPJNKdata`;
  
  ## create gti fits file:
  #`fcreate TMPJNKtempl TMPJNKdata TMPJNKfits1 extname="STDGTI"`;
  `fcreate TMPJNKtempl TMPJNKdata TMPJNKfits2 extname="GTI"`;
  
  ## append them to the file:
  #`fappend TMPJNKfits1"[1]" $out`;
  `fappend TMPJNKfits2"[1]" $out`;

  ## also remove TSTOP and add ONTIME keywords, just in case:
  #`fparkey dummy $out"[0]" -TSTOP`;
  `fparkey dummy $out"[1]" -TSTOP`;
  `fparkey $ontime $out"[0]" ONTIME add=yes`;
  `fparkey $ontime $out"[1]" ONTIME add=yes`;

  unlink "TMPJNKtempl", "TMPJNKdata", "TMPJNKfits1", "TMPJNKfits2";
}

print "Done.\n";