#!/proj/axaf/bin/perl -w

##################################################################################
#
# align_evt: Align event files prior to merging, based on detected X-ray sources
# Copyright (C) 2001  Tom Aldcroft
# 
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
##################################################################################

use Getopt::Long;
use POSIX;

# Version control

$VERSION = '$Id: align_evt,v 1.1.1.1 2002/05/21 18:47:34 aldcroft Exp $';  # ';

# Set up some constants

$d2r = atan(1.0) / 45.;
$d2a = 3600.;

# Now do the main processing

get_options();

run_celldetect();

make_cat();

foreach $evt2 (@evt2) {
    if ($evt2 ne $baseline_evt2) {
	($d_ra, $d_dec) = xcorr_sources($evt2);
	print_dmhedit($evt2, $d_ra, $d_dec);
    }
}

##****************************************************************************
sub get_options {
##****************************************************************************
    $par{quiet}     = 0;
    $par{snr}      = 4;
    $par{det_rad}  = 4.0;	# Arcmin
    $par{match_rad}= 3.0;	# Arcsec
    $par{force}    = 0;
    $par{longest}  = 0;

    GetOptions( \%par, 
	       'loud!',
	       'snr=f',
	       'det_rad=f',	# Radius for celldetect
	       'match_rad=f',
	       'force!',
	       'longest!',
	       'help!',
	       ) ||
		   exit( 1 );

    usage( 0 ) if $par{help};

    $par{loud} = ! $par{quiet};
    if ($par{loud}) {
	print "COMMAND LINE PARAMETERS\n";
	foreach (sort keys %par) {
	    printf "  %-16s = %s\n", $_, $par{$_};
	}
	print "\n";
    }

    @evt2 = @ARGV;		# Evt2 files are the rest of the arguments
}

##****************************************************************************
sub run_celldetect {
##****************************************************************************
    for $evt2 (@evt2) {
	my ($hdr) = fits_read_header($evt2,'EVENTS');
	$hdr{$evt2} = $hdr;	# Save for other routines

	my $src_evt2 = "source_evt2.fits";
	($src2 = $evt2) =~ s/evt/src/;
	$src2 =~ s/\.gz$//;
	$src2{$evt2} = $src2;
	next if (-e $src2 && !$par{force});

# Extract a small radius around each celestial location source in field

	print "Extracting photons at $hdr->{RA_NOM}, $hdr->{DEC_NOM}\n" if ($par{loud});
	($ra_hms, $dec_hms) = dec2hms($hdr->{RA_NOM}, $hdr->{DEC_NOM});
	if ($par{det_rad} < 10) {
	    $dmcopy = "dmcopy \"${evt2}[events][(x,y)=circle($ra_hms,$dec_hms,$par{det_rad}\')]\" $src_evt2 clobber=yes";
	} else {
	    $dmcopy = "cp $evt2 $src_evt2";
	}
	print "$dmcopy\n" if ($par{loud});
	system($dmcopy);

# Find sources in the small field

	$celldet = "celldetect $src_evt2 $src2 thresh=$par{snr} clobber=yes";
	print "$celldet\n" if ($par{loud});
	system($celldet);
	unlink ($src_evt2) if (-e $src_evt2);
    }

}

##***************************************************************************
sub make_cat {
##***************************************************************************
    my $exposure = 0;
    if ($par{longest}) {
	foreach $evt2 (@evt2) {
	    if ($hdr{$evt2}->{EXPOSURE} > $exposure) {
		$exposure = $hdr{$evt2}->{EXPOSURE};
		$baseline_src2 = $src2{$evt2};
		$baseline_evt2 = $evt2;
	    }
	}
    } else {
	$baseline_src2 = $src2{$evt2[0]};
	$baseline_evt2 = $evt2[0];
    }
    
    print "\n";
    print "Event file $baseline_evt2 used for baseline alignment\n\n";
    foreach (get_cell_detect($baseline_src2)) {
	my ($ra, $dec, undef, $snr) = split;
	next unless ($snr >= $par{snr});
	push @cat, {ra  => $ra,
		    dec => $dec,
		    counts => $counts,
		    snr    => $snr};
    }
}


##***************************************************************************
sub xcorr_sources {
##***************************************************************************
    $evt2 = shift;
    $src2 = $src2{$evt2};
    @celldet = get_cell_detect($src2);

    foreach (@celldet) {
	($ra, $dec, $counts, $snr) = split;
	next unless ($snr >= $par{snr});

	foreach $cat (@cat) {
	    $delta_rad = 3600 * radec_dist($ra, $dec, $cat->{ra}, $cat->{dec});
	    next unless ($delta_rad < $par{match_rad});
	    push @{$match{$src2}}, { ra => $ra,
				     dec => $dec,
				     snr => $snr,
				     d_dec => ($cat->{dec} - $dec) * $d2a,
				     d_ra  => ($cat->{ra} - $ra) * cos($cat->{dec}*$d2r) * $d2a,
				 };
	    last;
	}
    }

    @d_ra = ();
    @d_dec= ();
    print "\nFollowing sources in $evt2 match baseline catalog:\n" if ($par{loud});
    printf("%10s %10s %6s %6s %6s\n", "RA", "Dec", "D_RA", "D_Dec", "SNR") if ($par{loud});
    foreach $m (@{$match{$src2}}) {
	printf("%10.5f %10.5f %6.2f %6.2f %6.2f\n", $m->{ra}, $m->{dec},
	       $m->{d_ra}, $m->{d_dec}, $m->{snr}) if ($par{loud});
	push @d_ra, $m->{d_ra};
	push @d_dec, $m->{d_dec};
    }
    print "\n" if ($par{loud});
    
    return (@d_ra) ? (median(@d_ra), median(@d_dec)) : (0.0, 0.0)
}

##***************************************************************************
sub print_dmhedit {
##***************************************************************************
    my ($evt2, $d_ra, $d_dec) = @_;
    my $h = $hdr{$evt2};
    my %val = ();
    my %ind = ();

    foreach (keys %{$h}) {
	if ($h->{$_} =~ /(RA|DEC)---?TAN/) {
	    $coor = $1;
	    die "Unexpected keyword name $_ with value matching '/(RA|DEC)---?TAN/'\n"
		unless (/TCTYP(\d+)/);
	    $ind{$coor} = $1;
	    $val{$coor} = $h->{"TCRVL$1"};
	}
    }

    if ($par{loud}) {
	print "Current values are:  \n" ;
	print "RA:  TCRVL$ind{RA}  = $val{RA}\n";
	print "Dec: TCRVL$ind{DEC} = $val{DEC}\n\n";
    }
    $new_ra = $val{RA}+$d_ra/3600./cos($val{DEC}*$d2r);
    $new_dec = $val{DEC}+$d_dec/3600.;


    print "Cut and paste the following to align $evt2 to $baseline_evt2:\n\n";
    print " punlearn dmhedit\n";
    print " dmhedit infile=$evt2 filelist=none operation=add unit=degrees key=TCRVL$ind{RA} value=$new_ra\n";
    print " dmhedit infile=$evt2 filelist=none operation=add unit=degrees key=TCRVL$ind{DEC} value=$new_dec\n";
    print " dmhedit infile=$evt2 filelist=none operation=add unit=degrees key=RA_NOM value=$new_ra\n";
    print " dmhedit infile=$evt2 filelist=none operation=add unit=degrees key=DEC_NOM value=$new_dec\n";
    print "\n";
}

##***************************************************************************
sub get_cell_detect {
##***************************************************************************
# Read the cell detect source file, and then eliminate any sources
# which are within $par{excl_rad} of another source
    my $src2 = shift;

    my @celldet_good = ();
    $cmd = "dmlist \"${src2}\[col ra,dec,net_counts,snr\]\" data,clean";
    print "$cmd\n" if ($par{loud});
    my @celldet = `$cmd`;
    @celldet = grep !/^\#/, @celldet;
    for $i (0 .. $#celldet) {
	$bad = 0;
	for $j (0 .. $#celldet) {
	    next if ($i == $j);
	    ($ra1, $dec1) = split ' ', $celldet[$i];
	    ($ra2, $dec2) = split ' ', $celldet[$j];
	    $bad = 1 if (radec_dist($ra1, $dec1, $ra2, $dec2)*3600 < $par{match_rad});
	}
	push @celldet_good, $celldet[$i] unless ($bad);
    }
    
    return @celldet_good;
}

##***************************************************************************
sub median {
##***************************************************************************
    my @a = sort {$a <=> $b} @_;
    my $mid = int(@a / 2);
    return (@a % 2) ? $a[$mid] : ($a[$mid-1]+$a[$mid])/2.0;
}

###################################################################################
sub fits_read_header {
###################################################################################
    my @keys;
    my $file = shift;
    my $ext = shift;
    $file = "$file\[$ext\]" if ($ext);
    @keys = `dmlist $file header,raw,clean`;

    my $noquote = "[^']";
    my %keys = ();
    foreach (@keys) {
	if (/(\S+)\s*= '(${noquote}+)'/ or /(\S+)\s*= ([^\/]+)/) {   
	    my $key = $1;
	    my $keyval = $2;
	    $key =~ s/\s+$//;
	    $keyval =~ s/\s+$//;
	    $key =~ s/^\s+//;
	    $keyval =~ s/^\s+//;
	    $keys{$key} = $keyval;
	}
    }
    
    return \%keys;
}

##**********************************************************************
sub radec_dist {
##**********************************************************************
    my ($a1, $d1, $a2, $d2) = @_;
    my $d2r = 3.14159265358979/180.;

    return(0.0) if ($a1==$a2 && $d1==$d2);

    return acos( cos($d1*$d2r)*cos($d2*$d2r) * cos(($a1-$a2)*$d2r) +
              sin($d1*$d2r)*sin($d2*$d2r)) / $d2r;
}

##***************************************************************************
sub dec2hms {
##***************************************************************************
    # Converts between sexigesimal (HMS) and decimal RA and Dec.  The
    # direction of conversion is given by the number and form of inputs
    # Returns two-element array of (RA, Dec) in either case.

    $_ = join ' ', @_;
    s/[,:dhms]/ /g;
    @arg = split;

    if (@arg == 2) {
	my $ra = shift;
	my $dec = shift;
	my ($rah, $ram, $ras);
	my ($dec_sign, $decd, $decm, $decs);
	my ($ra_hms, $dec_hms);

	$ra += 360.0 if ($ra < 0);
	$ra /= 15.;
	$rah = floor($ra);
	$ram = floor(($ra - $rah) * 60.);
	$ras = ($ra - $rah - $ram / 60.) * 60. * 60.;
	
	$dec_sign = ($dec < 0);
	$dec = abs($dec);
	$decd = floor($dec);
	$decm = floor(($dec - $decd ) * 60.);
	$decs = ($dec - $decd - $decm / 60) * 60. * 60.;
	
	$ra_hms = sprintf "%d:%02d:%06.3f", $rah, $ram, $ras;
	$dec_hms = sprintf "%s%d:%02d:%05.2f", $dec_sign ? '-' : '+', $decd, $decm, $decs;

	return ($ra_hms,$dec_hms);
    } elsif (@arg == 6) {
	my ($rah, $ram, $ras, $decd, $decm, $decs) = @arg;
	$ra = 15.0*($rah + $ram/60. + $ras/3600.);
	$dec = abs($decd) + $decm/60. + $decs/3600.;
	$dec = -$dec if ($decd < 0.0);
	return (sprintf("%12.7f",$ra), sprintf("%12.6f", $dec));
    } else {
	warn "hms2dec: Error -- enter either 6 or 2 arguments\n";
    }
}


##***************************************************************************
sub usage
##***************************************************************************
{
  my ( $exit ) = @_;

  local $^W = 0;
  require Pod::Text;
  Pod::Text::pod2text( '-75', $0 );
  exit $exit;
}

=pod

=head1 SYNOPSIS

align_evt - Generate dmhedit commands to align Chandra event files
for multiple observations of the same field, based on detected sources.

=head1 USAGE

B<align_evt>  [I<options>] <event file 1> <event file 2> ...

=head1 DESCRIPTION


This program generates dmhedit commands to align Chandra event files
for multiple observations of the same field.  This is done in the
following way:

=over 4

=item 1.

Generate a list of detected X-ray sources within a specified radius of
the center for each event file using B<celldetect>.  This step is skipped if a
corresponding source file is already available (i.e. for an event file named
<root>_evt2.fits, the corresponding source file must be named
<root>_src2.fits).

=item 2.

A reference catalog of sources if taken from either the longest observation
or the first observation on the command (<event file 1>).  This is controlled
with the -longest switch.  The reference catalog is the baseline to which
the remaining event files are aligned.

=item 3.

The sources from each event file are cross-correlated with the reference 
catalog.  The RA and Dec offset for sources which match to within a specified 
match radius are recorded.  For each observation, the median of all the offsets
is used to create B<dmhedit> commands to align the event file to the
reference event file. 

=item 4.

The user can then cut and paste the B<dmhedit> commands to align
the event files.

=back

=head1 OPTIONS

=over 4

=item B<-help>

Print this help information.

=item B<-quiet>

Run in quiet mode to suppress the printing of processing messages

=item B<-snr> <Signal_to_Noise>

Specify the minimum signal to noise ratio for sources to be used in
cross correlation

=item B<-det_rad> <Detection_Radius>

Distance from center (arcmin) for sources to be used in cross-correlation.  
Default value is 4 arcmin. 

=item B<-match_rad> <Match_Radius>

Maximum radial offset (arcsec) for cross-correlation between a source and 
a candidate match in the reference catalog.  Default value is 3 arcsec.

=item B<-force>

Force program to run celldetect for each event file even if corresponding
source files are found.

=item B<-longest>

Use the longest observation for the reference catalog.  Default is to use
the first observation in the command line.

=back

=head1 EXAMPLES

    To run B<align_evt> you first need to be in the CIAO
    environment so tools like B<dmlist> and
    B<celldetect> are in your path.

Align the files acisf02234N001_evt2.fits and acisf02423N001_evt2.fits, using
the first one in the list for the baseline alignment:

  align_evt acisf02234N001_evt2.fits acisf02423N001_evt2.fits

Align all the evt2 files in the current directory, using the longest 
observation for the baseline alignment.  Run in 'quiet' mode to suppress 
printing processing information:

  align_evt -longest -quiet *evt2.fits

Align all the evt2 files in the current directory, using only bright sources
with a celldetect S/N ratio of at least 10:

  align_evt -snr 10.0  *evt2.fits

=head1 AUTHOR

Tom Aldcroft ( taldcroft@cfa.harvard.edu )

=cut