#! /usr1/local/bin/perl5 
$version ="1.12";
$date  ="2001-07-10";
$author = "kaa";

# v1.10   switched to CALDB tree for fefs and added option that if an arf is
#         given then rmfs will be calculated even if the chip is s3.
# v1.11   switched to using the new unified fefs (one per temperature).
# v1.12   added support for dates after 12/30/2000.

# This script creates ACIS rmf files for an arbitrary input spectral file.
# It works out which FEFs should be used. If the RMFs already exist in our
# database it uses addrmf to add them together, if not it runs mkrmf to
# create them.

# this is where we store the RMFs and the FEFs

$rmfdir = "/chandra/calib_data/acis/s3/rmf/00may/";
$fefdir = "/FTP/caldb/data/chandra/acis/cpf/fefs/";

# Get and parse the command line option 

    if ($ARGV[0]=~/-(\S+)/){
	$option=$1;
	if($option=~/h/){$help =1;}
	shift(@ARGV);
    }

if($help){
    print "\n acismakermf version $version $date $author\n\n";
    print "usage: acismakermf [-h] input_spectral_file chip_id output_rmf_file arf_file fef_file\n\n";
    print "This script is an intelligent builder of rmfs. It will run\n";
    print "mkrmf if necessary otherwise will add together rmfs in the\n";
    print "local database. If mkrmf is run then the arf_file is used to\n";
    print "define the energy bins for the response (axis1). If the arf is\n";
    print "specified on the command line then mkrmf will be run even if\n";
    print "pre-calculated rmfs are available. Instead of an arf file an\n";
    print "energy specification string, eg 0.1:11.0:0.01, can be given. The\n";
    print "fef file used can be set by the fef_file argument - if the \n";
    print "argument is not given then the default is used.\n";
    exit(0);
}

# we need HEASOFT...

if($ENV{'LHEASOFT'} !~/\S/) 
{
    print "\n You need to set up HEASOFT to use this script.\n\n";
    exit(0);
}

# we also need to make sure acisrmf is in the path...

$rmfexec = `which acisrmf`;
if(substr($rmfexec,0,2) eq 'no')
{
    print "\n You need to add acisrmf to your path.\n\n";
    exit(0);
} 

# Check that the user gave the correct number of arguments.

if(@ARGV != 3 && @ARGV != 4 && @ARGV != 5)
{
    print "\n usage : acismakermf [-h] input_spectral_file chip_id output_rmf_file arf_file fef_file\n";
    print " type acismakermf -h to get more information\n\n";
    exit(0);
}

$infile  = $ARGV[0];
($chip    = $ARGV[1]) =~ tr/a-z/A-Z/;
$outfile = $ARGV[2];
$arffile = $ARGV[3];
$feffile = $ARGV[4];

# convert the chip_id to a number if the user specifies I# or S#

if ( substr($chip,0,1) eq "I" ) {
    $chip = substr($chip,1,1);
} elsif ( substr($chip,0,1) eq "S" ) {
    $chip = substr($chip,1,1) + 4
}

# First read the DATE-OBS keyword from the input file to find the observation
# date

$command = "fkeypar $infile DATE-OBS\n";
system($command);
$obsdate = `pget fkeypar value`;
$year  = substr($obsdate,1,4);
$month = substr($obsdate,6,2);
$day   = substr($obsdate,9,2);

print "\n","Observation performed on ",$month,"/",$day,"/",$year,"\n";

# The following somewhat messy code decides the temperature of the focal plane.

if($year > 2000) {

  print "Detector temperature is -120 C\n";
  $ccdtemp = -120;

} elsif($year == 2000) {

  if( ($month == 1 && $day > 29) || $month > 1 ) {
     print "Detector temperature is -120 C\n";
     $ccdtemp = -120;
  } else {
     print "Detector temperature is -110 C\n";
     $ccdtemp = -110;
  }

} elsif($year == 1999) {

  if ( ($month == 9 && $day > 16) || $month > 10 ) {
     print "Detector temperature is -110 C\n";
     $ccdtemp = -110;
  } elsif ( ($month == 8 && $day > 23) || ($month == 9 && $day <= 16) ) {
     print "Detector temperature is -100 C and CTI is changing rapidly\n";
     $ccdtemp = -100;
  } elsif ( $month == 8 && $day > 11 && $day <= 23 ) {
     print "Detector temperature is -100 C\n";
     $ccdtemp = -100;
  } else {
     print "Detector temperature is -90 C\n";
     $ccdtemp = -90;
  }

}     

# Now set up to run acisrmf. If the chip is S3 and no arffile was specfied
# then we just have to add together the rmfs in the database.

if ( $chip == 7 && @ARGV < 4 ) {

# create the list of rmfs to add

    $command = "acisrmf infile=$infile rmfdir=$rmfdir ccdtemp=$ccdtemp outfile=temp.lst fefoption=no";
    print "\n",$command,"\n";
    system($command);

# now add them

    $command = "addrmf list=\"\@temp.lst\" rmffile=$outfile";
    print "\n",$command,"\n";
    system($command);

# and tidy up

    $command = "rm temp.lst";
    system($command);
    exit(0);

}

# It is not S3 so we will have to make the rmfs ourselves. Check that 
# CIAO is really set up...

if($ENV{'ASCDS_BIN'} !~/\S/) 
{
    print "\n You need to set up CIAO to use this script.\n\n";
    exit(0);
}

# We will run acisrmf to find out which parts of the FEF to use.

if ( $ccdtemp == -120 ) {

   $dstring = "2000-01-29";
   $vstring = "0002";

} elsif ( $ccdtemp == -110 ) {

   $dstring = "1999-09-16";
   $vstring = "0001";

} elsif ( $ccdtemp == -100 ) {

   $dstring = "1999-08-13";
   $vstring = "0001";

} elsif ( $ccdtemp == -90 ) {

   $dstring = "1999-07-22";
   $vstring = "0001";

}

# If the user has specified an input file then set that as fefnam to
# override the default 

if ( @ARGV != 5 ) {

   $command = "acisrmf infile=$infile rmfdir=$fefdir ccdtemp=$ccdtemp outfile=temp.lst fefoption=yes fefnam=default chipid=$chip dstring=$dstring vstring=$vstring";

} else {

   $command = "acisrmf infile=$infile rmfdir=$fefdir ccdtemp=$ccdtemp outfile=temp.lst fefoption=yes fefnam=$feffile chipid=$chip dstring=$dstring vstring=$vstring";

}
print "\n",$command,"\n";
system($command);

# temp.lst is now a list of the fefs to use to create rmfs. So open and read it
# also open temp2.lst to make a list of the rmfs to be fed to addrmf.

open(LISTFILE,"temp.lst");
open(RMFLIST,">temp2.lst");

# read the spectrum file to find out how many channels we have and set the
# binning appropriately.

$command = "fkeypar $infile\[SPECTRUM\] NAXIS2\n";
system($command);
$nchans = `pget fkeypar value`;
$nchans = int($nchans);

# read the CHANTYPE so we can check that the user is not mixing PI and PHA

$command = "fkeypar $infile\[SPECTRUM\] CHANTYPE\n";
system($command);
($chantype = `pget fkeypar value`) =~ tr/A-Z/a-z/;

$irmf = 1;
while ($fefinput=<LISTFILE>) {

# grab the FEF name and the fraction from the input string

    @input = split(/]/, $fefinput);
    $fefname = $input[0].']';
    $rmffrac = $input[1];

# check whether this fef is for pha or pi.

    if ( index($fefname,'pha') >= 0 ) {
        if ( index($chantype,'pi') >= 0 ) {
	    print "FEF appears to be PHA but spectrum is PI\n";
            exit(0);
        }
        $chantype = 'pha';
        $nbin = int(4096/$nchans);
        $unbinned = $nchans * $nbin
    } else {
        if ( index($chantype,'pha') >= 0 ) {
	    print "FEF appears to be PI but spectrum is PHA\n";
            exit(0);
        }
        $chantype = 'pi';
        $nbin = int(1024/$nchans);
        $unbinned = $nchans * $nbin
    }

# run mkrmf to make an rmf for this fef. Now the way to run mkrmf has 
# completely changed between CIAO v1 and v2. To find out which is running 
# check for the ATOMDB parameter which only exists in CIAO v2. Note that
# we check whether an energy specification string was given by looking for
# a colon in $arffile.

    if ($ENV{'ATOMDB'} !~/\S/) {

       if ( index($arffile,':') >= 0 ) {

          $command = "mkrmf in_file='$fefname' out_file=file$irmf.rmf output_fmt=legacy axis1=energy,'$arffile' axis2=$chantype,1,$nchans,$nchans,li low_thres = 1e-05";

      } else {

          $command = "mkrmf in_file='$fefname' out_file=file$irmf.rmf output_fmt=legacy axis1=energy,$arffile axis2=$chantype,1,$nchans,$nchans,li low_thres = 1e-05";

      }

    } else {

       if ( index($arffile,':') >= 0 ) {

          $command = "mkrmf infile='$fefname' outfile=file$irmf.rmf outfmt=legacy axis1=\"energy=$arffile\" axis2=\"$chantype=1:$unbinned:$nbin\" thresh=1e-05 logfile=rmf.log";

       } else {

          $command = "mkrmf infile='$fefname' outfile=file$irmf.rmf outfmt=legacy axis1=\"energy=grid($arffile\[specresp\]\[cols energ_lo,energ_hi\])\" axis2=\"$chantype=1:$unbinned:$nbin\" thresh=1e-05 logfile=rmf.log";

       }

    }

    print "\n",$command,"\n";
    system($command);

    print RMFLIST "file",$irmf,".rmf"," ",$rmffrac;

    ++$irmf;

}

close(LISTFILE);
close(RMFLIST);

# now add the RMFs

$command = "addrmf list=\"\@temp2.lst\" rmffile=$outfile";
print "\n",$command,"\n";
system($command) &&
die "Oh dear, addrmf failed - I have saved the rmf files as file*.rmf\n and the input to addrmf as temp2.lst\n";

while ($irmf > 1) {

    --$irmf;
    $command = "rm -f file$irmf.rmf";
    system($command);
}

$command = "rm -f temp.lst temp2.lst rmf.log";
system($command);