Index: trunk/pstamp/scripts/detectability_respond.pl
===================================================================
--- trunk/pstamp/scripts/detectability_respond.pl	(revision 27640)
+++ trunk/pstamp/scripts/detectability_respond.pl	(revision 27640)
@@ -0,0 +1,519 @@
+#! /usr/bin/env perl
+#
+#
+# Create a response to a MOPS_DETECTABILITY_QUERY
+#
+#
+
+use strict;
+use warnings;
+
+use Getopt::Long qw( GetOptions );
+use Pod::Usage qw( pod2usage );
+
+use PS::IPP::Config qw($PS_EXIT_SUCCESS
+                       $PS_EXIT_UNKNOWN_ERROR
+                       $PS_EXIT_SYS_ERROR
+                       $PS_EXIT_CONFIG_ERROR
+                       $PS_EXIT_PROG_ERROR
+                       $PS_EXIT_DATA_ERROR
+                       $PS_EXIT_TIMEOUT_ERROR
+                       metadataLookupStr
+                       metadataLookupBool
+                       caturi
+                       );
+
+use Sys::Hostname;
+use Carp;
+use File::Basename;
+use File::Copy;
+use PS::IPP::PStamp::RequestFile qw( :standard );
+use PS::IPP::PStamp::Job qw( :standard );
+use Astro::FITS::CFITSIO qw( :constants );
+Astro::FITS::CFITSIO::PerlyUnpacking(1);
+
+use IPC::Cmd 0.36 qw( can_run run );
+
+#
+# Set up
+#
+my $host = hostname();
+
+my $EXTVER = 1.0;
+my ($missing_tools,$request_file,$tmp_dir);
+my ($req_id,$req_name,$product);
+my ($dbname,$need_magic);
+my ($input,$output,$workdir,$verbose,$save_temps);
+GetOptions(
+    'input=s'          =>     \$request_file,
+    'output=s'        =>      \$output,
+    'dbname=s'        =>      \$dbname,
+    'workdir=s'       =>      \$workdir,
+    'verbose'         =>      \$verbose,
+    'save-temps'      =>      \$save_temps,
+    ) or pod2usage(2);
+
+pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
+pod2usage( -msg => "Required options: --input --output --dbname",
+	   -exitval => 3,
+    ) unless
+    defined $request_file and defined $output and defined $dbname;
+
+my $detect_query_read = can_run('detect_query_read') or (warn "Can't find detect_query_read" and $missing_tools = 1);
+my $psphotForced = can_run('psphotForced') or (warn "Can't find psphotForced" and $missing_tools = 1);
+if ($missing_tools) {
+    warn("Can't find required tools.");
+    exit($PS_EXIT_CONFIG_ERROR);
+}
+
+unless (defined $verbose) {
+    $verbose = 0;
+}
+
+
+my $ipprc = PS::IPP::Config->new();
+$tmp_dir = "/data/${host}.0/tmp/";
+#
+# Parse Input
+#
+
+my $dqr_command = "$detect_query_read --input $request_file";
+my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+    run(command => $dqr_command, verbose => $verbose);
+unless ($success) {
+    warn("Unable to perform $dqr_command error code: $error_code");
+}
+
+my @column_names = ();
+my %response = ();
+my $section = '';
+my $Nrows = 0;
+foreach my $entry (split /\n/, (join "", @$stdout_buf)) {
+    if ($entry =~ /^#/) {
+	@column_names = split /\s+/, $entry;
+	shift(@column_names);  # Dump the hash sign.
+	if ($section eq 'HEADER') {
+	    $section = 'CONTENT';
+	}
+	else {
+	    $section = 'HEADER';
+	}
+    }
+    else {
+	# HEADER: 
+	# QUERY_ID FPA_ID MJD_OBS FILTER OBSCODE STAGE
+	# CONTENT:
+	# ROWNUM RA1_DEG DEC1_DEG RA2_DEG DEC2_DEG MAG
+	my @columns = split /\s+/, $entry;
+	for (my $i = 0; $i <= $#columns; $i++) {
+	    $response{$section}{$column_names[$i]}[$Nrows] = $columns[$i];
+	    print "$section $column_names[$i] $Nrows $columns[$i]\n";
+	}
+	$Nrows++;
+    }
+}
+
+#
+# Identify target images.  This should properly collate targets on a single imfile.
+#
+my %image_list_hash;
+for (my $i = 1; $i < $Nrows; $i++) {
+    print "$i $Nrows $response{CONTENT}{RA1_DEG}[$i] $response{CONTENT}{DEC1_DEG}[$i]\n";
+
+    my $image_set_tmp  = find_image_set($response{HEADER}{FPA_ID}[0],$response{HEADER}{STAGE}[0],
+					$response{HEADER}{MJD_OBS}[0],$response{HEADER}{FILTER}[0],
+					$response{CONTENT}{RA1_DEG}[$i],
+					$response{CONTENT}{DEC1_DEG}[$i],
+					$response{CONTENT}{ROWNUM}[$i],$verbose);
+    print "=== $image_set_tmp->{IMAGE}\n    $image_set_tmp->{PSF}\n    $image_set_tmp->{MASK}\n    $image_set_tmp->{WEIGHT}\n    $image_set_tmp->{COORDINATES}\n    $image_set_tmp->{ROWNUM}\n";
+    # This appends, assuming that if we get an image, we also get the identical psf/mask/weight/etc.
+    $image_list_hash{$image_set_tmp->{IMAGE}}{IMAGE} = $image_set_tmp->{IMAGE};
+    $image_list_hash{$image_set_tmp->{IMAGE}}{PSF} = $image_set_tmp->{PSF};
+    $image_list_hash{$image_set_tmp->{IMAGE}}{MASK} = $image_set_tmp->{MASK};
+    $image_list_hash{$image_set_tmp->{IMAGE}}{WEIGHT} = $image_set_tmp->{WEIGHT};
+    $image_list_hash{$image_set_tmp->{IMAGE}}{CLASS_ID} = $image_set_tmp->{CLASS_ID};
+    $image_list_hash{$image_set_tmp->{IMAGE}}{EXTENSION_BASE} = $image_set_tmp->{EXTENSION_BASE};
+    push @{ $image_list_hash{$image_set_tmp->{IMAGE}}{COORDINATES} }, $image_set_tmp->{COORDINATES};
+    push @{ $image_list_hash{$image_set_tmp->{IMAGE}}{ROWNUM} }, $image_set_tmp->{ROWNUM};
+    $image_list_hash{$image_set_tmp->{IMAGE}}{TARGETS} = 
+	"$tmp_dir/detectability.$response{HEADER}{STAGE}[0].$response{HEADER}{FPA_ID}[0].targets";
+    $image_list_hash{$image_set_tmp->{IMAGE}}{OUTROOT} = 
+	"$tmp_dir/detectability.$response{HEADER}{STAGE}[0].$response{HEADER}{FPA_ID}[0]";
+}
+#exit(2);
+my $i = 0;
+my @image_list;
+foreach my $k (keys %image_list_hash) {
+    # write coordinates for photometry to a disk file.  Probably need to do this as a tempfile object.
+    open(T,">$image_list_hash{$k}{TARGETS}") || die "Could not open output TARGET file $image_list_hash{$k}{TARGETS}";
+    for (my $j = 0; $j <= $#{ $image_list_hash{$k}{COORDINATES} }; $j++) {
+	print T "$image_list_hash{$k}{COORDINATES}[$j]\n";
+    }
+    close(T);
+    $image_list[$i] = $image_list_hash{$k};
+    $i++;
+}
+
+#
+# Construct and run psphot
+#
+for (my $i = 0; $i <= $#image_list; $i++) {
+    my $psf_file    = $image_list[$i]->{PSF};
+    my $image_file  = $image_list[$i]->{IMAGE};
+    my $mask_file   = $image_list[$i]->{MASK};
+    my $weight_file = $image_list[$i]->{WEIGHT};
+    my $target_file = $image_list[$i]->{TARGETS};
+    my $out_root    = $image_list[$i]->{OUTROOT};
+    my $psphot_cmd = "$psphotForced -psf $psf_file ";
+    $psphot_cmd .= "-file $image_file -mask $mask_file -variance $weight_file ";
+    $psphot_cmd .= "-srctext $target_file $out_root";
+    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+	run(command => $psphot_cmd, verbose => $verbose);
+    unless ($success) {
+	warn ("Unable to perform $psphot_cmd error code: $error_code");
+    }
+}
+
+#
+# Convert psphot output to response
+#
+my @response_entries = ();
+
+foreach my $k (keys %image_list_hash) {
+    my $class_id = $image_list_hash{$k}{CLASS_ID};
+    my $cmf = $image_list_hash{$k}{OUTROOT} . ".${class_id}.cmf";
+    my ($detect_N_col,$detect_F_col,@detectability_data) = read_cmf_file($cmf,$image_list_hash{$k}{EXTENSION_BASE});
+    
+    if ($#detectability_data != $#{ $image_list_hash{$k}{ROWNUM} }) {
+	my_die("Number of measured values does not match number of requested values.");
+    }
+    for (my $i = 0; $i < $#detectability_data; $i++) {
+	push @response_entries, create_response_entry($image_list_hash{$k}{ROWNUM}[$i],
+						      $detect_N_col,$detect_F_col,
+						      $detectability_data[$i]);
+    }
+}
+
+write_response_file($output,
+		    $response{HEADER}{QUERY_ID}[0],$response{HEADER}{FPA_ID}[0],
+		    $response{HEADER}{MJD_OBS}[0],$response{HEADER}{filter}[0],
+		    $response{HEADER}{obscode}[0],
+		    @response_entries);
+
+#
+# Add to datastore
+#
+# my $finish_command = "dquery_finish.pl --req_id $req_id --req_name $req_name --product $product";
+# ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+#     run(command => $finish_command, verbose => $verbose);
+# unless ($success) {
+#     warn ("Unable to perform $finish_command error code: $error_code");
+# }
+
+#
+# Cleanup
+#
+#
+# Utilities
+#
+sub find_image_set {
+    my $FPA_ID = shift;
+    my $stage  = lc(shift);
+    my $mjd    = shift;
+    my $filter = shift;
+    my $ra     = shift;
+    my $dec    = shift;
+    my $index  = shift;
+    my $verbose = shift;
+    
+    my $option_mask |= 1;
+    $option_mask |= $PSTAMP_SELECT_IMAGE;
+    $option_mask |= $PSTAMP_SELECT_MASK;
+    $option_mask |= $PSTAMP_SELECT_VARIANCE;
+    $option_mask |= $PSTAMP_SELECT_PSF;
+    my $need_magic = 1;
+    my $mjd_min = $mjd;
+    my $mjd_max = $mjd + 1;
+
+    my @images = locate_images($ipprc,$dbname,"bycoord",$stage,
+			       undef,undef,undef,$option_mask,$need_magic,
+			       $ra,$dec,$mjd_min,$mjd_max,$filter . ".00000",undef,$verbose);  
+#     my @images = lookup($ipprc,$dbname,$req_type,$img_type,
+# 			$id,$tess_id,$componenent,$need_magic,$dateobs_begin,$dateobs_end,
+# 			$filter,$data_group,$option_mask,$verbose);
+    my @keys_i_care = ('image','mask','psf','exp_name','stage','weight','cmf','smf');
+    foreach my $i (@images) {
+	foreach my $j (@{ $i }) {
+#	    foreach my $k (keys %{ ${ $i }[$j] }) {
+	}
+    }
+
+    my %image_info  = ();
+
+    foreach my $i (@images) {
+	foreach my $j (@{ $i }) {
+	    if ($stage eq 'diff') {
+		#Diffs hold both exposure names, if they are defined (so a WW or WS diff).  
+		unless ((defined(${ $j }{exp_name_1}) && (${ $j }{exp_name_1} eq $FPA_ID))||
+			(defined(${ $j }{exp_name_2}) && (${ $j }{exp_name_2} eq $FPA_ID))) {
+		    next;
+		}
+		push @keys_i_care, 'exp_name_1','exp_name_2';
+	    }
+	    elsif ($stage eq 'stack') {
+		# Stacks do not store the exposure name information...I'll have to think up something smart here.
+		next;
+	    }
+	    else {
+		if (${ $j }{exp_name} ne $FPA_ID) {
+		    next;
+		}
+	    }
+#	    foreach my $k (@keys_i_care) {
+ 	    foreach my $k (keys %{ $j }) {
+		print "$i $j $k ${ $j }{$k}\n";
+	    }
+
+	    my $image_name = ${ $j }{image};
+	    $image_info{IMAGE} = $image_name;
+	    $image_info{PSF} = ${ $j }{psf};
+	    $image_info{MASK} = ${ $j }{mask};
+	    $image_info{WEIGHT} = ${ $j }{weight};
+	    if (exists(${ $j }{astrom})) {
+		($image_info{COORDINATES},$image_info{EXTENSION_BASE}) = 
+		    calculate_image_coordinates(${ $j }{astrom},$ra,$dec);
+		$image_info{CLASS_ID} = ${ $j }{class_id};
+	    }
+	    else {
+		($image_info{COORDINATES},$image_info{EXTENSION_BASE}) = 
+		    calculate_image_coordinates(${ $j }{cmf},$ra,$dec);
+		$image_info{CLASS_ID} = 'fpa';
+	    }
+	    $image_info{ROWNUM} = $index;
+	}
+    }
+    return(\%image_info);
+}
+# Taken largely from detect_query_read
+sub calculate_image_coordinates {
+    my $cmf = shift;
+    my $ra = shift;
+    my $dec = shift;
+    my $x = $ra;
+    my $y = $dec;
+    my $c = '';
+    
+    my $tempfile = "$tmp_dir/detect.cic.$$";
+    open(TEMP,">$tempfile");
+    print TEMP "$ra $dec\n";
+    close(TEMP);
+    my $command = "ppCoord -astrom $cmf -radec $tempfile";
+    
+    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+	run(command => $command, verbose => $verbose);
+    unless ($success) {
+	warn("Unable to perform $command error code: $error_code");
+    }
+    print ">>>$command<<<\n";
+    my @response = split /\n/, (join "", @$stdout_buf);
+    foreach my $line (@response) {
+	print ">>$line<<\n";
+	my ($r_ra,$r_dec,$trash,$r_x,$r_y,$r_chip) = split /\s+/, $line;
+	if (($r_ra == $ra)&&($r_dec == $dec)) {
+	    $x = $r_x;
+	    $y = $r_y;
+	    $c = $r_chip;
+	}
+    }
+    unlink($tempfile);
+    return("$x $y",$c);
+}
+sub read_cmf_file {
+    my $cmf_file = shift;
+    my $class_id = shift;
+    my @outdata = ();
+    my $extname = $class_id . ".psf";
+
+    my ($detect_F_col,$detect_N_col,$detect_flux_col);
+    print "RCF: $cmf_file $class_id $extname\n";
+    my $status = 0;
+    my $inFits = Astro::FITS::CFITSIO::open_file( $cmf_file, READONLY, $status ); # Open CMF
+    check_fitsio($status);
+
+    $inFits->movnam_hdu(BINARY_TBL, $extname, 0, $status) and check_fitsio($status);
+
+    # This is the only data we actually care about
+    my $column_defs = [
+	{ name => 'PSF_QF',   type => '1E', writetype => TDOUBLE },
+	{ name => 'PSF_NPIX', type => '1J', writetype => TLONG },
+	{ name => 'PSF_INST_MAG', type => '1E', writetype => TDOUBLE },
+	];
+    my @colNames; 
+    my @colTypes;
+    my @colWriteType;
+    my %colData;
+    foreach my $colSpec (@$column_defs) {
+	push @colNames, $colSpec->{name};
+	push @colTypes, $colSpec->{type};
+	push @colWriteType, $colSpec->{writetype};
+    }
+    my $numRows;
+    $inFits->get_num_rows($numRows, $status) and check_fitsio($status);
+
+    foreach my $col (@$column_defs) {
+	my ($col_num,$col_type,$col_data);
+	$inFits->get_colnum(0, $col->{name}, $col_num, $status) and check_fitsio($status);
+	$inFits->get_coltype($col_num, $col_type, undef, undef, $status) and check_fitsio($status);
+	$inFits->read_col($col_type, $col_num, 1, 1, $numRows, 0, $col_data, undef, $status) and check_fitsio($status);
+	$colData{$col->{name}} = $col_data;
+	if ($col->{name} eq 'PSF_QF') {
+	    $detect_F_col = $col_num;
+	}
+	if ($col->{name} eq 'PSF_NPIX') {
+	    $detect_N_col = $col_num;
+	}
+	if ($col->{name} eq 'PSF_INST_MAG') {
+	    $detect_flux_col = $col_num;
+	}
+    }
+    $inFits->close_file( $status ) and check_fitsio($status);    
+
+    for (my $i = 0; $i < $numRows; $i++) {
+	for (my $j = 0; $j <= $#colNames; $j++) {
+	    if ($j == $detect_flux_col) {
+		$outdata[$i][$j] = 10**(-0.4 * $colData{$colNames[$j]}->[$i]);
+	    }
+	    else {
+		$outdata[$i][$j] = $colData{$colNames[$j]}->[$i];
+	    }
+	}
+    }
+
+    return($detect_N_col,$detect_F_col,@outdata);
+}
+
+sub create_response_entry {
+    my $rownum = shift;
+    my $detect_N_col = shift;
+    my $detect_F_col = shift;
+
+    my @cmf_data = @_;
+    
+    my $entry;
+    @{ $entry } = ();
+    
+    
+    push @{ $entry }, $rownum;
+    push @{ $entry }, $cmf_data[$detect_N_col];
+    push @{ $entry }, $cmf_data[$detect_F_col];
+    
+    return($entry);
+}
+# A lot of this pulled verbatim from Bill's detect_response_create file
+sub write_response_file {
+    my $outfile = shift;
+    my $query_id = shift;
+    my $FPA_ID = shift;
+    my $MJD_OBS = shift;
+    my $filter = shift;
+    my $obscode = shift;
+    my @response_entries = @_;
+#    my ($query_id,$FPA_ID,$MJD_OBS,$filter,$obscode);
+    my $status = 0;
+    # Specification of columns to write
+    my $columns = [
+	# matching rownum from detectability original request
+        { name => 'ROWNUM',   type => '20A', writetype => TSTRING }, 
+        # number of pixels used in hypothetical PSF for the query detection
+        { name => 'DETECT_N', type => 'V',   writetype => TULONG },
+        # detectibility, indicating the fraction of PSF pixels detetable by IPP
+        { name => 'DETECT_F', type => 'D',   writetype => TDOUBLE },
+	];
+    
+    # Header translation table
+    my $headers = {
+	'QUERY_ID' => { name => 'QUERY_ID', type => TSTRING, 
+                        comment => 'MOPS Query ID for this batch query' },
+	'FPA_ID' => { name => 'FPA_ID',   type => TSTRING, 
+		      comment => 'original FPA_ID used at ingest' },
+    };
+
+    # Parse the list of columns
+    my @colNames;                   # Names of columns
+    my @colTypes;                   # Types of columns
+    my %colData;                    # Data for each column
+    my @colWriteType;                 # type to use to write
+    foreach my $colSpec ( @$columns) {
+	push @colNames, $colSpec->{name};
+	push @colTypes, $colSpec->{type};
+	push @colWriteType, $colSpec->{writetype};
+	$colData{$colSpec->{name}} = [];
+    }
+
+    my $numRows = $#response_entries + 1;
+    my $inHeader = { };
+
+    # Hack to force the data to match teh detect_response_create formats
+    $inHeader->{QUERY_ID}->{value} = $query_id;
+    $inHeader->{FPA_ID}->{value} = $FPA_ID;
+    $inHeader->{MJD_OBS}->{value} = $MJD_OBS;
+    $inHeader->{FILTER}->{value} = $filter;
+    $inHeader->{OBSCODE}->{value} = $obscode;
+    
+    for (my $i = 0; $i <= $#response_entries; $i++) {
+	push @{$colData{'ROWNUM'}}, $response_entries[$i][0];
+	push @{$colData{'DETECT_N'}}, $response_entries[$i][1];
+	push @{$colData{'DETECT_F'}}, $response_entries[$i][2];
+    }
+
+    # Back to detect_response_create:
+    my $outFits = Astro::FITS::CFITSIO::create_file( $output, $status );
+    check_fitsio( $status );
+
+    $outFits->create_img( 16, 0, undef, $status );
+    check_fitsio( $status );
+    my $EXTNAME = 'MOPS_DETECTABILITY_RESPONSE';
+    $outFits->create_tbl( BINARY_TBL(), $numRows, scalar @colNames, \@colNames, \@colTypes, undef, $EXTNAME, $status);
+    check_fitsio( $status );
+
+    # TODO: get the extension version number from somewhere common
+    $outFits->write_key( TINT, 'EXTVER', $EXTVER, 'Version of this Extension', $status );
+    check_fitsio( $status );
+
+    # Write the Extension keywords
+    foreach my $keyword ( keys %$headers ) {
+	my $value = $inHeader->{$keyword}->{value}; # Header keyword value
+	unless (defined $value) {
+	    print "Can't find header keyword $keyword\n";
+	    next;
+	}
+	$value =~ s/\'//g;
+	my $name    = $headers->{$keyword}->{name}; # New name
+	my $type    = $headers->{$keyword}->{type}; # Type
+	my $comment = $headers->{$keyword}->{comment}; # Comment
+	$outFits->write_key( $type, $name, $value, $comment, $status );
+	check_fitsio( $status );
+    }
+    for (my $i = 0; $i < scalar @colNames; $i++) {
+	my $colName = $colNames[$i];# Column name
+	my $writeType = $colWriteType[$i];
+	print "$colName $writeType $i $colData{$colName}[$i] $status $numRows\n";
+	$outFits->write_col( $writeType, $i + 1, 1, 1, $numRows, $colData{$colName}, $status );
+	print "$colName $writeType $i $colData{$colName}[$i] $status $numRows\n";
+	check_fitsio( $status );
+    }
+    $outFits->close_file( $status );
+    return($status);
+}
+
+# From Astro::FITS::CFITSIO demo
+sub check_fitsio
+{
+    my $status = shift;         # Status of FITSIO calls
+
+    if ($status != 0) {
+        my $msg;                # Message to output
+        Astro::FITS::CFITSIO::fits_get_errstatus( $status , $msg );
+        die "CFITSIO error: $msg\n";
+    }
+}