#!/usr/local/bin/perl 
select STDERR; $| = 1; select STDOUT; $| = 1;

#############################################################
#                    Large_detect.pl                        #
#                                                           #
#  This Perl script runs CIAO celldetect and wavdetect      #
#  programs repeatedly over an ACIS event list covering     #
#  a large region (e.g. the entire ACIS-I array).           #
#  The resulting source lists are consolidated into         #
#  master lists for scientific study.                       #
#                                                           #
#    ver 0    Y.Maeda                Sep. 3rd, 1999         #
#    ver 1    Y.Maeda & E.Feigelson  March 5 2000           #
#############################################################
 
###  The command line for running this script is:
###      Large_detect.pl  input.evt  out_root binfactor
###  input.evt is an ACIS event list FITS file
###  out_root is name for output files
###  binfactor is an input for celldetect & wavdetect
###      binfactor=1 is recommended

### The user may reset the following parameters which define
### the size of image segments to be analyzed in each run of celldetect
### and wavdetect.  Note wavdetect wants 2**n size segments.

$imgsize = 1024;
$margin = 128;

###  Output files are:
###      outroot_celldetect_srclist.fits
###      outroot_wavdetect_srclist.fits
###  and many output files from the detection programs
###      for individual image segments.  

###  Read command line 

system("rm srclist_celldetect srclist_wavdetect");

$pardir = "./";

$evtfile = $ARGV[0];
##print "$evtfile\n";

$outroot = $ARGV[1];
$binfactor = $ARGV[2];

### Calculate (x,y) center of the image under study 

$fstatistic = "fstatistic $evtfile x -";
print "$fstatistic\n";
system("$fstatistic");

$rangexmin = `pget fstatistic min`; chop($rangexmin);
$rangexmax = `pget fstatistic max`; chop($rangexmax);

$fstatistic = "fstatistic $evtfile y -";
print "$fstatistic\n";
system("$fstatistic");

$rangeymin = `pget fstatistic min`; chop($rangeymin);
$rangeymax = `pget fstatistic max`; chop($rangeymax);

$rangex = $rangexmax - $rangexmin;
$rangey = $rangeymax - $rangeymin;

$centerx = ($rangexmax + $rangexmin) / 2; $centerx = int($centerx);
$centery = ($rangeymax + $rangeymin) / 2; $centery = int($centery);

$nx = $rangex / ($imgsize - $margin) / $binfactor ; $nx = $nx + 1; $nx = int($nx);
$ny = $rangey / ($imgsize - $margin) / $binfactor ; $ny = $ny + 1; $ny = int($ny);

###  Obtain parameters for celldetect and wavdetect from parameter files
###  First search for *.par files in local directory, then in ~/PFILES

$HOME = `printenv HOME`; chop($HOME);

$_ = `ls $pardir/celldetect.par`;
if($_=~ /celldetect.par/)
{
}
else
{
    system("cp $HOME/PFILES/celldetect.par $pardir");
    system("pedit $pardir/celldetect.par");
}
    system("cp $pardir/celldetect.par $HOME/PFILES ");

$_ = `ls $pardir/wavdetect.par`;
if($_=~ /wavdetect.par/)
{
}
else
{
    system("cp $HOME/PFILES/wavdetect.par $pardir");
    system("pedit $pardir/wavdetect.par");
}
    system("cp $pardir/wavdetect.par $HOME/PFILES ");


	$sigthresh = `pget $pardir/wavdetect sigthresh`; chop($sigthresh);
	$bkgsigthresh = `pget $pardir/wavdetect bkgsigthresh`; chop($bkgsigthresh);
	$scales = `pget $pardir/wavdetect scales`; chop($scales);

###  Loop over image segments.  For example:
#	$xmin = ($d*1024 - $margin*$d + 1)*$binfactor;
#	$xmax = (($d+1)*1024 - $margin*$d)*$binfactor;
#	$ymin = ($e*1024 - $margin*$e + 1)*$binfactor;
#	$ymax = (($e+1)*1024 - $margin*$e)*$binfactor;

$i = -1;
for($d=1;$d<=$nx;$d++)
{
    
    for($e=1;$e<=$ny;$e++)
    {
	$i = $i + 1;

	$xmin = ($centerx - ($imgsize - $margin) * $nx / 2 - $margin / 2 + ($d-1)*($imgsize - $margin) + 1)*$binfactor;
	$xmax = ($centerx - ($imgsize - $margin) * $nx / 2 - $margin / 2 + ($d-1)*($imgsize - $margin) + $imgsize)*$binfactor;
	$ymin = ($centery - ($imgsize - $margin) * $ny / 2 - $margin / 2 + ($e-1)*($imgsize - $margin) + 1)*$binfactor;
	$ymax = ($centery - ($imgsize - $margin) * $ny / 2 - $margin / 2 + ($e-1)*($imgsize - $margin) + $imgsize)*$binfactor;
	
	$cutxmin = $xmin + $margin/2;
	$cutxmax = $xmax - $margin/2;
	$cutymin = $ymin + $margin/2;
	$cutymax = $ymax - $margin/2;

	$region[$d][$e] = "\[bin x=$xmin:$xmax:$binfactor, y=$ymin:$ymax:$binfactor\]";

	$command = $region[$d][$e];

### Celldetect is run for two reasons: to obtain celldetect sources
### if desired, and to determine which segments of the image have no events

$bkgvalue = `pget $pardir/celldetect bkgvalue`; chop($bkgvalue);
$thresh = `pget $pardir/celldetect thresh`; chop($thresh);
    $celldetect = "celldetect \"$evtfile$command\" $outroot\_$i\_celldetect\_srclist.fits thresh=$thresh bkgvalue=$bkgvalue";
print "$celldetect\n";
system("$celldetect");
	
### The following lines avoid running wavdetect when no data are present.  

$lscom = "$outroot\_$i\_celldetect\_srclist.fits";
$ls = `ls $outroot\_$i\_celldetect\_srclist.fits`;

if($ls=~ /$lscom/)
	{
### The following lines throw away duplicate celldetect sources in overlapping margin regions.

$fselect = "fselect $outroot\_$i\_celldetect\_srclist.fits $outroot\_$i\_celldetect\_srclist\_shapeup.fits \"X.gt.$cutxmin&&X.le.$cutxmax&&Y.gt.$cutymin&&Y.le.$cutymax\"";
print "$fselect\n";
system("$fselect");
system("echo $outroot\_$i\_celldetect\_srclist\_shapeup.fits \>\> srclist_celldetect");

### Run wavdetect over one image segment.

	$wavdetect = "wavdetect infile=\"$evtfile$command\" scellfile=$outroot\_$i\_wavdetect\_scell.fits imagefile=$outroot\_$i\_wavdetect\_image.fits defnbkgfile=$outroot\_$i\_wavdetect\_defnbkg.fits outfile=$outroot\_$i\_wavdetect\_srclist.fits sigthresh=$sigthresh bkgsigthresh=$bkgsigthresh scales=\"$scales\"";
	print "$wavdetect\n";
	system("$wavdetect");


	$lscom = "$outroot\_$i\_wavdetect\_srclist.fits";
	$ls = `ls $outroot\_$i\_wavdetect\_srclist.fits`;

### The following lines throw away duplicate wavdetect sources in overlapping margin regions.

	if($ls=~ /$lscom/)
	{
	    $fselect = "fselect $outroot\_$i\_wavdetect\_srclist.fits $outroot\_$i\_wavdetect\_srclist\_shapeup.fits \"X.gt.$cutxmin&&X.le.$cutxmax&&Y.gt.$cutymin&&Y.le.$cutymax\"";
	    print "$fselect\n";
	    system("$fselect");
	    
	    system("echo $outroot\_$i\_wavdetect\_srclist\_shapeup.fits \>\> srclist_wavdetect");
	}
    }
    }
}

### Consolidate the source lists from individual image segments into master
### celldetect and wavdetect source lists.

$fmerge = "fmerge infiles=\@srclist\_celldetect outfile=$outroot\_celldetect\_srclist.fits columns=\"-\"";
print "$fmerge\n";
	system("$fmerge");

$fmerge = "fmerge infiles=\@srclist\_wavdetect outfile=$outroot\_wavdetect\_srclist.fits columns=\"-\"";
print "$fmerge\n";
	system("$fmerge");

###  Sort the mast source lists by RA.

$fsort = "fsort infile=$outroot\_celldetect\_srclist.fits columns=\"ra dec\" method=heap";
print "$fsort\n";
system("$fsort");

$fsort = "fsort infile=$outroot\_wavdetect\_srclist.fits columns=\"ra dec\" method=heap";
print "$fsort\n";
	system("$fsort");

### END