#!/usr/bin/perl -w

use SAS;
require "esas_parameters_init.pl";
require "esas_GVariables.pm";

&init();

&run();

#========================================================================
sub init() {
#.. Create object to store parameters
#.. Define Global Variables
    use vars qw($gv);
    
    $gv = esas_GVariables->new();   
    
	#.. Init Program Parameters
	&cheesebands_parameters_init(@ARGV);
		

	# $task    = "cheese-bands";
	# $version = "4.0";
	# $date    = "2010-03-02";
	# $author  = "S. L. Snowden";
	# print "\n";
	# print "Task: $task \n";
	# print "Version: $version \n";
	# print "Date: $date \n";
	# print "Author: $author \n";

	# Included in SAS: Date
	# SAS Version: 
	# A. Ibarra Ibaibarriaga
}

## @method void run()
# Run Cheese
# @return Void

#========================================================================
sub run() {
	# Set the variables		

	$nexp=0;
	foreach ( $gv->Cheesebands_prefixm ) {
			$nexp = $nexp + 1;
			$prefix[$nexp] = $_;
			print "mos prefix $nexp: $prefix[$nexp]\n";
	}
	
	$nexpp=0;
	foreach ($gv->Cheesebands_prefixp ) {			
			$nexpp = $nexpp + 1;
			$prefixp[$nexpp] = $_;
			print "pn prefix $nexpp: $prefixp[$nexpp]\n";
	}
 	
	$verb  = $gv->Cheesebands_verb;
	print "verb: $verb\n";
	$scale = $gv->Cheesebands_scale;
	print "scale: $scale\n";
	$ratet = $gv->Cheesebands_ratet;
	$ratet  = 1.0e-14 * $ratet;
	print "ratet: $ratet\n";
	$rates = $gv->Cheesebands_rates;
	$rates  = 1.0e-14 * $rates;
	print "rates: $rates\n";
	$rateh = $gv->Cheesebands_rateh;
	$rateh  = 1.0e-14 * $rateh;
	print "rateh: $rateh\n";
	$dist = $gv->Cheesebands_dist;
	print "dist: $dist\n";
	
	$clobber = $gv->Cheesebands_clobber;
	print "clobber: $clobber\n";
	
	$nb=0;
	foreach ( $gv->Cheesebands_elowlist ) {
			$nb = $nb + 1;
			$elow[$nb] = $_;
	}
	$nb=0;
	foreach ( $gv->Cheesebands_ehighlist ) {
			$nb = $nb + 1;
			$ehigh[$nb] = $_;
			print "Band $nb: $elow[$nb] $ehigh[$nb]\n";
	}

	open( JUNK_HNDL, ">>command.csh" );

	#  First create an attitude file if necessary

	$affile = "atthk.fits";
	if ( -r $affile ) {
		print "attitude file already exists\n";
		print "use the old file\n";
		print "   \n";
	}
	else {
		$comm = "atthkgen atthkset=atthk.fits timestep=1 ";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
	}

	#  Next extract the MOS images in soft and hard bands

	for ( $nn = 1 ; $nn <= $nexp ; $nn++ ) {
          print "nn: $nn\n";
          print "prefix: $prefix[$nn]\n";
          $infile = "mos".$prefix[$nn]."-obj-im-s.fits";
          if(-r $infile) {
             print "Soft event image already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "evselect table=mos".$prefix[$nn]."-clean.fits:EVENTS "
                     . "withfilteredset=yes "
		     . "expression='(PATTERN<=12)&&(FLAG == 0)&&"
                     . "(PI in [".$elow[1].":".$ehigh[1]."])' "
		     . "filtertype=expression imagebinning='imageSize' imagedatatype='Int32' "
		     . "imageset=".$infile." squarepixels=yes "
                     . "ignorelegallimits=yes withxranges=yes withyranges=yes "
                     . "xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 "
		     . "ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 "
		     . "updateexposure=yes filterexposure=yes verbosity=" . $verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	    $comm = "evselect table=mos".$prefix[$nn]."-clean.fits:EVENTS "
              . "withfilteredset=yes "
	      . "expression='(PATTERN<=12)&&(FLAG == 0)&&"
              . "(PI in [".$elow[1].":".$ehigh[1]."])' "
	      . "filtertype=expression imagebinning='imageSize' imagedatatype='Int32' "
	      . "imageset=".$infile." squarepixels=yes "
              . "ignorelegallimits=yes withxranges=yes withyranges=yes "
              . "xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 "
	      . "ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 "
	      . "updateexposure=yes filterexposure=yes verbosity=" . $verb;
	    print JUNK_HNDL "$comm\n";
	    print JUNK_HNDL "   \n";
	    print "$comm\n";
	    print "   \n";
	    system "$comm";
            print "   \n";
          }
		
          $infile = "mos".$prefix[$nn]."-obj-im-h.fits";
          if(-r $infile) {
             print "Hard event image already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "evselect table=mos".$prefix[$nn]."-clean.fits:EVENTS "
                     . "withfilteredset=yes "
		     . "expression='(PATTERN<=12)&&(FLAG == 0)&&"
                     . "(PI in [".$elow[2].":".$ehigh[2]."])' "
		     . "filtertype=expression imagebinning='imageSize' imagedatatype='Int32' "
		     . "imageset=".$infile." squarepixels=yes "
                     . "ignorelegallimits=yes withxranges=yes withyranges=yes "
                     . "xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 "
		     . "ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 "
		     . "updateexposure=yes filterexposure=yes verbosity=" . $verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "evselect table=mos".$prefix[$nn]."-clean.fits:EVENTS "
                 . "withfilteredset=yes "
	         . "expression='(PATTERN<=12)&&(FLAG == 0)&&"
                 . "(PI in [".$elow[2].":".$ehigh[2]."])' "
		     . "filtertype=expression imagebinning='imageSize' imagedatatype='Int32' "
		     . "imageset=".$infile." squarepixels=yes "
                     . "ignorelegallimits=yes withxranges=yes withyranges=yes "
                     . "xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 "
		     . "ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 "
		     . "updateexposure=yes filterexposure=yes verbosity=" . $verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
          }

	#  Now make mos exposure maps for the FOV images

          $infile = "mos".$prefix[$nn]."-exp-im-s.fits";
          if(-r $infile) {
             print "Soft band exposure map already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "eexpmap attitudeset=atthk.fits "
		     . "eventset=mos".$prefix[$nn]."-clean.fits:EVENTS "
		     . "expimageset=".$infile." imageset=mos".$prefix[$nn]."-obj-im-s.fits "
		     . "pimax=1250 pimin=400 withdetcoords=no verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "eexpmap attitudeset=atthk.fits "
	          . "eventset=mos".$prefix[$nn]."-clean.fits:EVENTS "
	          . "expimageset=".$infile." imageset=mos".$prefix[$nn]."-obj-im-s.fits "
	          . "pimax=1250 pimin=400 withdetcoords=no verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
          }

          $infile = "mos".$prefix[$nn]."-exp-im-h.fits";
          if(-r $infile) {
             print "Hard band exposure map already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "eexpmap attitudeset=atthk.fits "
		     . "eventset=mos".$prefix[$nn]."-clean.fits:EVENTS "
		     . "expimageset=".$infile." "
		     . "imageset=mos".$prefix[$nn]."-obj-im-h.fits "
		     . "pimax=1250 pimin=400 withdetcoords=no verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "eexpmap attitudeset=atthk.fits "
	          . "eventset=mos".$prefix[$nn]."-clean.fits:EVENTS "
	          . "expimageset=".$infile." "
	          . "imageset=mos".$prefix[$nn]."-obj-im-h.fits "
	          . "pimax=1250 pimin=400 withdetcoords=no verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
          }

	#  Make masks

          $infile = "mos".$prefix[$nn]."-mask-im-s.fits";
          if(-r $infile) {
             print "Soft band mask already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "emask detmaskset=".$infile." "
		     . "expimageset=mos".$prefix[$nn]."-exp-im-s.fits "
		     . "threshold1=0.1 threshold2=0.5 verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "emask detmaskset=".$infile." "
	          . "expimageset=mos".$prefix[$nn]."-exp-im-s.fits "
	          . "threshold1=0.1 threshold2=0.5 verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
          }

          $infile = "mos".$prefix[$nn]."-mask-im-h.fits";
          if(-r $infile) {
             print "Hard band mask already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "emask detmaskset=".$infile." "
		     . "expimageset=mos".$prefix[$nn]."-exp-im-h.fits "
		     . "threshold1=0.1 threshold2=0.5 verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "emask detmaskset=".$infile." "
	          . "expimageset=mos".$prefix[$nn]."-exp-im-h.fits "
	          . "threshold1=0.1 threshold2=0.5 verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
          }
	}

	#  Next extract the pn images in the selected band eV band

	for ( $nnp = 1 ; $nnp <= $nexpp ; $nnp++ ) {
          print "nnp: $nnp\n";
          print "prefixp: $prefixp[$nnp]\n";
	
          $infile = "pn".$prefixp[$nnp]."-obj-im-s.fits";
          if(-r $infile) {
             print "PN soft image already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "evselect table=pn".$prefixp[$nnp]."-clean.fits:EVENTS "
		     . "withfilteredset=yes "
		     . "expression='(PATTERN<=12)&&(FLAG == 0)&&"
                     . "(PI in [".$elow[1].":".$ehigh[1]."])&&(DETY in [-16510:14345])' "
		     . "filtertype=expression imagebinning='imageSize' imagedatatype='Int32' "
		     . "imageset=".$infile." squarepixels=yes "
		     . "ignorelegallimits=yes withxranges=yes withyranges=yes "
		     . "xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 "
		     . "ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 "
		     . "updateexposure=yes filterexposure=yes verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "evselect table=pn".$prefixp[$nnp]."-clean.fits:EVENTS "
	          . "withfilteredset=yes "
	          . "expression='(PATTERN<=12)&&(FLAG == 0)&&"
                  . "(PI in [".$elow[1].":".$ehigh[1]."])&&(DETY in [-16510:14345])' "
	          . "filtertype=expression imagebinning='imageSize' imagedatatype='Int32' "
	          . "imageset=".$infile." squarepixels=yes "
	          . "ignorelegallimits=yes withxranges=yes withyranges=yes "
	          . "xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 "
	          . "ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 "
	          . "updateexposure=yes filterexposure=yes verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
	  }
	
          $infile = "pn".$prefixp[$nnp]."-obj-im-h.fits";
          if(-r $infile) {
             print "PN hard image already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "evselect table=pn".$prefixp[$nnp]."-clean.fits:EVENTS "
		     . "withfilteredset=yes "
		     . "expression='(PATTERN<=12)&&(FLAG == 0)&&"
                     . "(PI in [".$elow[2].":".$ehigh[2]."])&&(DETY in [-16510:14345])' "
		     . "filtertype=expression imagebinning='imageSize' imagedatatype='Int32' "
		     . "imageset=".$infile." squarepixels=yes "
		     . "ignorelegallimits=yes withxranges=yes withyranges=yes "
		     . "xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 "
		     . "ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 "
		     . "updateexposure=yes filterexposure=yes verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "evselect table=pn".$prefixp[$nnp]."-clean.fits:EVENTS "
	          . "withfilteredset=yes "
	          . "expression='(PATTERN<=12)&&(FLAG == 0)&&"
                  . "(PI in [".$elow[2].":".$ehigh[2]."])&&(DETY in [-16510:14345])' "
	          . "filtertype=expression imagebinning='imageSize' imagedatatype='Int32' "
	          . "imageset=".$infile." squarepixels=yes "
	          . "ignorelegallimits=yes withxranges=yes withyranges=yes "
	          . "xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 "
	          . "ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 "
	          . "updateexposure=yes filterexposure=yes verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
	  }

	#  Now make pn exposure maps for the FOV image

          $infile = "pn".$prefixp[$nnp]."-exp-im-s.fits";
          if(-r $infile) {
             print "PN soft exposure map already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "eexpmap attitudeset=atthk.fits "
		     . "eventset=pn".$prefixp[$nnp]."-clean.fits:EVENTS "
		     . "expimageset=".$infile." "
		     . "imageset=pn".$prefixp[$nnp]."-obj-im-s.fits "
		     . "pimax=1250 pimin=400 withdetcoords=no verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "eexpmap attitudeset=atthk.fits "
	          . "eventset=pn".$prefixp[$nnp]."-clean.fits:EVENTS "
	          . "expimageset=".$infile." "
	          . "imageset=pn".$prefixp[$nnp]."-obj-im-s.fits "
	          . "pimax=1250 pimin=400 withdetcoords=no verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
	  }

          $infile = "pn".$prefixp[$nnp]."-exp-im-h.fits";
          if(-r $infile) {
             print "PN hard exposure map already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "eexpmap attitudeset=atthk.fits "
		     . "eventset=pn".$prefixp[$nnp]."-clean.fits:EVENTS "
		     . "expimageset=".$infile." "
		     . "imageset=pn".$prefixp[$nnp]."-obj-im-h.fits "
		     . "pimax=1250 pimin=400 withdetcoords=no verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "eexpmap attitudeset=atthk.fits "
	          . "eventset=pn".$prefixp[$nnp]."-clean.fits:EVENTS "
	          . "expimageset=".$infile." "
	          . "imageset=pn".$prefixp[$nnp]."-obj-im-h.fits "
	          . "pimax=1250 pimin=400 withdetcoords=no verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
	  }

	#  Make a mask

          $infile = "pn".$prefixp[$nnp]."-mask-im-s.fits";
          if(-r $infile) {
             print "PN soft mask already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "emask detmaskset=".$infile." "
		     . "expimageset=pn".$prefixp[$nnp]."-exp-im-s.fits "
		     . "threshold1=0.1 threshold2=0.5 verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "emask detmaskset=".$infile." "
	          . "expimageset=pn".$prefixp[$nnp]."-exp-im-s.fits "
	          . "threshold1=0.1 threshold2=0.5 verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
	  }

          $infile = "pn".$prefixp[$nnp]."-mask-im-h.fits";
          if(-r $infile) {
             print "PN hard mask already exists: $infile\n";
             if($clobber==1) {
                print "Clobber the existing file\n";
                print "   \n";
		$comm = "emask detmaskset=".$infile." "
		     . "expimageset=pn".$prefixp[$nnp]."-exp-im-h.fits "
		     . "threshold1=0.1 threshold2=0.5 verbosity=".$verb;
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
             } else {
                print "Use the existing file\n";
             }
          } else {
	     $comm = "emask detmaskset=".$infile." "
	          . "expimageset=pn".$prefixp[$nnp]."-exp-im-h.fits "
	          . "threshold1=0.1 threshold2=0.5 verbosity=".$verb;
	     print JUNK_HNDL "$comm\n";
	     print JUNK_HNDL "   \n";
	     print "$comm\n";
	     print "   \n";
	     system "$comm";
	     print "   \n";
	  }
	}

	#  Set the parameter names

	$tempdet = 0;
	for ( $nn = 1 ; $nn <= $nexp ; $nn++ ) {
		if ( $nn == 1 ) {
			$imlist    = "mos" . $prefix[$nn] . "-obj-im-s.fits " 
			      . "mos" . $prefix[$nn] . "-obj-im-h.fits";
			$explist   = "mos" . $prefix[$nn] . "-exp-im-s.fits " 
			      . "mos" . $prefix[$nn] . "-exp-im-h.fits";
			$bkglist   = "mos" . $prefix[$nn] . "-bkgimage-s.fits " 
			      . "mos" . $prefix[$nn] . "-bkgimage-h.fits";
			$piminlist = $elow[1]." ".$elow[2];
			$pimaxlist = $ehigh[1]." ".$ehigh[2];
#			$piminlist = "400 2000";
#			$pimaxlist = "1250 7200";
			$ecflist   = "1.2 1.2";
			$det       = substr( $prefix[$nn], 0, 1 );

			if ( $tempdet != $det ) {
				$masklist = "mos" . $prefix[$nn] . "-mask-im-h.fits";
				$tempdet  = $det;
			}
		}
		else {
			$imlist    = $imlist . " mos" . $prefix[$nn] . "-obj-im-s.fits " 
			      . "mos" . $prefix[$nn] . "-obj-im-h.fits";
			$explist   = $explist . " mos" . $prefix[$nn] . "-exp-im-s.fits " 
			      . "mos" . $prefix[$nn] . "-exp-im-h.fits";
			$bkglist   = $bkglist . " mos" . $prefix[$nn] . "-bkgimage-s.fits " 
			      . "mos" . $prefix[$nn] . "-bkgimage-h.fits";
			$piminlist = $piminlist." ".$elow[1]." ".$elow[2];
			$pimaxlist = $pimaxlist." ".$ehigh[1]." ".$ehigh[2];
#			$piminlist = "400 2000";
#			$pimaxlist = "1250 7200";
			$ecflist   = $ecflist . " 1.2 1.2";
			$det       = substr( $prefix[$nn], 0, 1 );

			if ( $tempdet != $det ) {
				$masklist = $masklist . " mos" . $prefix[$nn] . "-mask-im-h.fits";
				$tempdet  = $det;
			}
		}
	}

	if ( $nexpp >= 1 ) {
		for ( $nnp = 1 ; $nnp <= $nexpp ; $nnp++ ) {
			if ( $nnp == 1 ) {
				if ( $nexp >= 1 ) {
					$imlist    = $imlist . " pn" . $prefixp[$nnp] . "-obj-im-s.fits " 
					     . "pn" . $prefixp[$nnp] . "-obj-im-h.fits";
					$explist   = $explist . " pn" . $prefixp[$nnp] . "-exp-im-s.fits " 
					     . "pn" . $prefixp[$nnp] . "-exp-im-h.fits";
					$bkglist   = $bkglist . " pn" . $prefixp[$nnp] . "-bkgimage-s.fits " 
					     . "pn" . $prefixp[$nnp] . "-bkgimage-h.fits";
					$piminlist = $piminlist." ".$elow[1]." ".$elow[2];
					$pimaxlist = $pimaxlist." ".$ehigh[1]." ".$ehigh[2];
#					$piminlist = "400 2000";
#					$pimaxlist = "1250 7200";
					$masklist  = $masklist . " pn" . $prefixp[$nnp] . "-mask-im-h.fits";
					$ecflist   = $ecflist . " 3.2 3.2";
				}
				else {
					$imlist    = "pn" . $prefixp[$nnp] . "-obj-im-s.fits " 
					     . "pn" . $prefixp[$nnp] . "-obj-im-h.fits";
					$explist   = "pn" . $prefixp[$nnp] . "-exp-im-s.fits " 
					     . "pn" . $prefixp[$nnp] . "-exp-im-h.fits";
					$bkglist   = "pn" . $prefixp[$nnp] . "-bkgimage-s.fits " 
					     . "pn" . $prefixp[$nnp] . "-bkgimage-h.fits";
					$piminlist = $elow[1]." ".$elow[2];
					$pimaxlist = $ehigh[1]." ".$ehigh[2];
#					$piminlist = "400 2000";
#					$pimaxlist = "1250 7200";
					$masklist  = "pn" . $prefixp[$nnp] . "-mask-im-h.fits";
					$ecflist   = "3.2 3.2";
				}
			}
			else {
				$imlist  = $imlist . " pn" . $prefixp[$nnp] . "-obj-im-s.fits " 
					     . "pn" . $prefixp[$nnp] . "-obj-im-h.fits";
				$explist = $explist . " pn" . $prefixp[$nnp] . "-exp-im-s.fits " 
					     . "pn" . $prefixp[$nnp] . "-exp-im-h.fits";
				$bkglist = $bkglist . " pn" . $prefixp[$nnp] . "-bkgimage.fits " 
					     . "pn" . $prefixp[$nnp] . "-bkgimage-h.fits";
				$piminlist = $piminlist." ".$elow[1]." ".$elow[2];
				$pimaxlist = $pimaxlist." ".$ehigh[1]." ".$ehigh[2];
				$masklist  =
				  $masklist . " pn" . $prefixp[$nnp] . "-mask-im.fits-h";
				$ecflist = $ecflist . " 3.2 3.2";
			}
		}
	}

	#  Do the first box detection

	print "imlist: $imlist\n";
	print JUNK_HNDL "imlist: $imlist\n";
	print "explist: $explist\n";
	print JUNK_HNDL "explist: $explist\n";
	print "masklist: $masklist\n";
	print JUNK_HNDL "masklist: $masklist\n";
	print "bkglist: $bkglist\n";
	print JUNK_HNDL "bkglist: $bkglist\n";
	print "piminlist: $piminlist\n";
	print JUNK_HNDL "piminlist: $piminlist\n";
	print "pimaxlist: $pimaxlist\n";
	print JUNK_HNDL "pimaxlist: $pimaxlist\n";
	print "ecflist: $ecflist\n";
	print JUNK_HNDL "ecflist: $ecflist\n";

	$comm =
	    "eboxdetect detmasksets='".$masklist."' withdetmask=yes imagebuffersize=900 ".
	    "expimagesets='".$explist."' withexpimage=yes nruns=3 likemin=8 ". 
            "boxsize=5 ecf='".$ecflist."' imagesets='".$imlist."' verbosity=".$verb." ". 
            "boxlistset=boxlist.fits pimin='".$piminlist."' pimax='".$pimaxlist."'";

	print "$comm\n";
	print "   \n";
	system "$comm";
	print "   \n";
	print JUNK_HNDL "$comm\n";
	print JUNK_HNDL "   \n";

	# Create the spline background maps for the second call to eboxdetect

	for ( $nn = 1 ; $nn <= $nexp ; $nn++ ) {
		$comm = "esplinemap boxlistset=boxlist.fits "
		  . "imageset=mos".$prefix[$nn]."-obj-im-s.fits "
		  . "expimageset=mos".$prefix[$nn]."-exp-im-s.fits "
		  . "withexpimage=yes "
		  . "detmaskset=mos".$prefix[$nn]."-mask-im-s.fits withdetmask=yes "
		  . "bkgimageset=mos".$prefix[$nn]. "-bkgimage-s.fits "
		  . "pimin=400 pimax=1250 idband=1 verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "esplinemap boxlistset=boxlist.fits "
		  . "imageset=mos".$prefix[$nn]."-obj-im-h.fits "
		  . "expimageset=mos".$prefix[$nn]."-exp-im-h.fits "
		  . "withexpimage=yes "
		  . "detmaskset=mos".$prefix[$nn]."-mask-im-h.fits "
		  . "withdetmask=yes "
		  . "bkgimageset=mos".$prefix[$nn]. "-bkgimage-h.fits "
		  . "pimin=2000 pimax=7200 idband=2 verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
	}

	for ( $nnp = 1 ; $nnp <= $nexpp ; $nnp++ ) {
		$comm = "esplinemap boxlistset=boxlist.fits "
		  . "imageset=pn".$prefixp[$nnp]."-obj-im-s.fits "
		  . "expimageset=pn".$prefixp[$nnp]."-exp-im-s.fits "
		  . "withexpimage=yes "
		  . "detmaskset=pn".$prefixp[$nnp]."-mask-im-s.fits "
		  . "withdetmask=yes "
		  . "bkgimageset=pn".$prefixp[$nnp]."-bkgimage-s.fits "
		  . "pimin=400 pimax=1250 idband=1 verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
		$comm =
		    "esplinemap boxlistset=boxlist.fits "
		  . "imageset=pn".$prefixp[$nnp]."-obj-im-h.fits "
		  . "expimageset=pn".$prefixp[$nnp]."-exp-im-h.fits "
		  . "withexpimage=yes "
		  . "detmaskset=pn".$prefixp[$nnp]."-mask-im-h.fits "
		  . "withdetmask=yes "
		  . "bkgimageset=pn".$prefixp[$nnp]."-bkgimage-h.fits "
		  . "pimin=400 pimax=1250 idband=2 verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";
	}

#  Do the second box detection

	print "imlist: $imlist\n";
	print JUNK_HNDL "imlist: $imlist\n";
	print "explist: $explist\n";
	print JUNK_HNDL "explist: $explist\n";
	print "masklist: $masklist\n";
	print JUNK_HNDL "masklist: $masklist\n";
	print "bkglist: $bkglist\n";
	print JUNK_HNDL "bkglist: $bkglist\n";
	print "piminlist: $piminlist\n";
	print JUNK_HNDL "piminlist: $piminlist\n";
	print "pimaxlist: $pimaxlist\n";
	print JUNK_HNDL "pimaxlist: $pimaxlist\n";
	print "ecflist: $ecflist\n";
	print JUNK_HNDL "ecflist: $ecflist\n";

	$comm = "eboxdetect detmasksets='".$masklist."' withdetmask=yes ".
            "imagebuffersize=900 expimagesets='".$explist."' withexpimage=yes ".
            "nruns=3 likemin=8 boxsize=5 ecf='".$ecflist."' imagesets='".$imlist."' ".
            "boxlistset=boxlist-f.fits pimin='".$piminlist."' pimax='". $pimaxlist."' ". 
            "boxsize=5 ecf='".$ecflist. "' imagesets='". $imlist."' usemap=yes ".
            "bkgimagesets='".$bkglist."' verbosity=".$verb;

	print "$comm\n";
	print "   \n";
	system "$comm";
	print "   \n";
	print JUNK_HNDL "$comm\n";
	print JUNK_HNDL "   \n";

#  Do the maximum likelihood detection

	$comm = "emldetect imagesets='".$imlist."' bkgimagesets='".$bkglist."' ".
	    "boxlistset=boxlist-f.fits mllistset=emllist.fits mlmin=8 ".
	    "dmlextmin=10 maxextent=20 expimagesets='".$explist."' ".
	    "withexpimage=yes detmasksets='".$masklist."' withdetmask=yes ".
            "imagebuffersize=900 scut=0.9 ecut=0.68 ecf='".$ecflist."' ".
            "mergedlistset=mergedlist.fits pimin='".$piminlist."' ".
	    "pimax='".$pimaxlist."' verbosity=".$verb;

	print "$comm\n";
	print "   \n";
	system "$comm";
	print "   \n";
	print JUNK_HNDL "$comm\n";
	print JUNK_HNDL "   \n";
#@	}

#	#  Make a backup copy of the EML source list
#
#        $comm = "cp emllist.fits emllist-ori.fits";
#	print "$comm\n";
#	print "   \n";
#	system "$comm";
#	print "   \n";
#
#	#  Fix the source list for arbitrarily low background rates
#
#       $comm = "/Home/lhea1/snowden/esas/eml_fix emllist.fits";
#	print "$comm\n";
#	print "   \n";
#	system "$comm";
#	print "   \n";

	#  Create the background regions selection expression

	for ( $nn = 1 ; $nn <= $nexp ; $nn++ ) {

		$det = substr( $prefix[$nn], 0, 1 ) + 1;

		$comm = "region eventset=mos".$prefix[$nn]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(FLUX >= ".$rates.")&&(ID_INST == ".$det.")&&"
		  . "(ID_BAND == 1)&&(DIST_NN >= ".$dist. ")' "
		  . "bkgregionset=mos".$prefix[$nn]."-bkg_region-det-s.fits  "
		  . "energyfraction=".$scale." radiusstyle=enfrac "
		  . "outunit=detxy verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=mos".$prefix[$nn]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(FLUX >= ".$rateh.")&&(ID_INST == ".$det.")&&"
		  . "(ID_BAND == 2)&&(DIST_NN >= ".$dist. ")' "
		  . "bkgregionset=mos".$prefix[$nn]."-bkg_region-det-h.fits  "
		  . "energyfraction=".$scale." radiusstyle=enfrac "
		  . "outunit=detxy verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=mos".$prefix[$nn]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(FLUX >= ".$ratet.")&&(ID_INST == ".$det.")&&"
		  . "(ID_BAND == 0)&&(DIST_NN >= ".$dist. ")' "
		  . "bkgregionset=mos".$prefix[$nn]."-bkg_region-det.fits  "
		  . "energyfraction=".$scale." radiusstyle=enfrac "
		  . "outunit=detxy verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=mos".$prefix[$nn]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(DIST_NN >= ".$dist.")&&(FLUX >= ".$rates.")"
		  . "&&(ID_INST == ".$det.")&&(ID_BAND == 1)' "
		  . "bkgregionset=mos".$prefix[$nn]."-bkg_region-sky-s.fits  "
		  . "radiusstyle=contour bkgratestyle=col nosrcellipse=yes "
		  . "bkgfraction=".$scale." outunit=xy verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=mos".$prefix[$nn]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(DIST_NN >= ".$dist.")&&(FLUX >= ".$rateh.")"
		  . "&&(ID_INST == ".$det.")&&(ID_BAND == 2)' "
		  . "bkgregionset=mos".$prefix[$nn]."-bkg_region-sky-h.fits  "
		  . "radiusstyle=contour bkgratestyle=col nosrcellipse=yes "
		  . "bkgfraction=".$scale." outunit=xy verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=mos".$prefix[$nn]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(DIST_NN >= ".$dist.")&&(FLUX >= ".$ratet.")"
		  . "&&(ID_INST == ".$det.")&&(ID_BAND == 0)' "
		  . "bkgregionset=mos".$prefix[$nn]."-bkg_region-sky.fits  "
		  . "radiusstyle=contour bkgratestyle=col nosrcellipse=yes "
		  . "bkgfraction=".$scale." outunit=xy verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		#  Create the masks

		$comm = "make_mask inimage=mos".$prefix[$nn]."-obj-im-s.fits "
		  . "inmask=mos".$prefix[$nn]."-mask-im-s.fits "
		  . "outmask=mos".$prefix[$nn]."-cheese-s.fits "
		  . "reglist=mos".$prefix[$nn]."-bkg_region-sky-s.fits";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "make_mask inimage=mos".$prefix[$nn]."-obj-im-h.fits "
		  . "inmask=mos".$prefix[$nn]."-mask-im-h.fits "
		  . "outmask=mos".$prefix[$nn]."-cheese-h.fits "
		  . "reglist=mos".$prefix[$nn]."-bkg_region-sky-h.fits";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "make_mask inimage=mos".$prefix[$nn]."-obj-im-h.fits "
		  . "inmask=mos".$prefix[$nn]."-mask-im-h.fits "
		  . "outmask=mos".$prefix[$nn]."-cheese.fits "
		  . "reglist=mos".$prefix[$nn]."-bkg_region-sky.fits";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

	}

	#  Do it again for the pn

	for ( $nnp = 1 ; $nnp <= $nexpp ; $nnp++ ) {
		$det = 1;
		$comm = "region eventset=pn".$prefixp[$nnp]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(DIST_NN >= ".$dist.")&&(FLUX >= ".$rates.")"
		  . "&&(ID_INST == ".$det.")&&(ID_BAND == 1)' "
		  . "bkgregionset=pn".$prefixp[$nnp]."-bkg_region-det-s.fits  "
		  . "energyfraction=".$scale." radiusstyle=enfrac "
		  . "outunit=detxy verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=pn".$prefixp[$nnp]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(DIST_NN >= ".$dist.")&&(FLUX >= ".$rateh.")"
		  . "&&(ID_INST == ".$det.")&&(ID_BAND == 2)' "
		  . "bkgregionset=pn".$prefixp[$nnp]."-bkg_region-det-h.fits  "
		  . "energyfraction=".$scale." radiusstyle=enfrac "
		  . "outunit=detxy verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=pn".$prefixp[$nnp]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(DIST_NN >= ".$dist.")&&(FLUX >= ".$ratet.")"
		  . "&&(ID_INST == ".$det.")&&(ID_BAND == 0)' "
		  . "bkgregionset=pn".$prefixp[$nnp]."-bkg_region-det.fits  "
		  . "energyfraction=".$scale." radiusstyle=enfrac "
		  . "outunit=detxy verbosity=".$verb;
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=pn".$prefixp[$nnp]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(DIST_NN >= ".$dist.")&&(FLUX >= ".$rates.")"
		  . "&&(ID_INST == ".$det.")&&(ID_BAND == 1)' "
		  . "bkgregionset=pn".$prefixp[$nnp]."-bkg_region-sky-s.fits  "
		  . "radiusstyle=contour bkgratestyle=col nosrcellipse=yes "
		  . "bkgfraction=".$scale." outunit=xy verbosity=".$verb;

		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=pn".$prefixp[$nnp]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(DIST_NN >= ".$dist.")&&(FLUX >= ".$rateh.")"
		  . "&&(ID_INST == ".$det.")&&(ID_BAND == 2)' "
		  . "bkgregionset=pn".$prefixp[$nnp]."-bkg_region-sky-h.fits "
		  . "radiusstyle=contour bkgratestyle=col nosrcellipse=yes "
		  . "bkgfraction=".$scale." outunit=xy verbosity=".$verb;

		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "region eventset=pn".$prefixp[$nnp]."-clean.fits "
		  . "operationstyle=global srclisttab=emllist.fits:SRCLIST "
		  . "expression='(DIST_NN >= ".$dist.")&&(FLUX >= ".$ratet.")"
		  . "&&(ID_INST == ".$det.")&&(ID_BAND == 0)' "
		  . "bkgregionset=pn".$prefixp[$nnp]."-bkg_region-sky.fits "
		  . "radiusstyle=contour bkgratestyle=col nosrcellipse=yes "
		  . "bkgfraction=".$scale." outunit=xy verbosity=".$verb;

		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		#  Create the masks

		$comm = "make_mask inimage=pn".$prefixp[$nnp]."-obj-im-s.fits "
		  . "inmask=pn".$prefixp[$nnp]."-mask-im-s.fits "
		  . "outmask=pn".$prefixp[$nnp]."-cheese-s.fits "
		  . "reglist=pn".$prefixp[$nnp]."-bkg_region-sky-s.fits";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "make_mask inimage=pn".$prefixp[$nnp]."-obj-im-h.fits "
		  . "inmask=pn".$prefixp[$nnp]."-mask-im-h.fits "
		  . "outmask=pn".$prefixp[$nnp]."-cheese-h.fits "
		  . "reglist=pn".$prefixp[$nnp]."-bkg_region-sky-h.fits";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		$comm = "make_mask inimage=pn".$prefixp[$nnp]."-obj-im-h.fits "
		  . "inmask=pn".$prefixp[$nnp]."-mask-im-h.fits "
		  . "outmask=pn".$prefixp[$nnp]."-cheese.fits "
		  . "reglist=pn".$prefixp[$nnp]."-bkg_region-sky.fits";
		print "$comm\n";
		print "   \n";
		system "$comm";
		print "   \n";
		print JUNK_HNDL "$comm\n";
		print JUNK_HNDL "   \n";

		#  Get rid of a few more extra files

		system "rm -f filtered.fits";

	}
}