IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 8, 2010, 5:53:28 PM (16 years ago)
Author:
watersc1
Message:

More robust and tested detectability response creator. Seems to run well, and create all the necessary table entries. Now to tie it in with the rest of the postage stamp server.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/pstamp/scripts/detectability_respond.pl

    r27640 r27642  
    99use warnings;
    1010
     11use Sys::Hostname;
     12use Carp;
     13use File::Basename;
     14use File::Copy;
     15use IPC::Cmd 0.36 qw( can_run run );
    1116use Getopt::Long qw( GetOptions );
    1217use Pod::Usage qw( pod2usage );
    13 
     18use File::Temp qw( tempfile tempdir);
     19
     20use PS::IPP::PStamp::RequestFile qw( :standard );
     21use PS::IPP::PStamp::Job qw( :standard );
    1422use PS::IPP::Config qw($PS_EXIT_SUCCESS
    1523                       $PS_EXIT_UNKNOWN_ERROR
     
    2331                       caturi
    2432                       );
    25 
    26 use Sys::Hostname;
    27 use Carp;
    28 use File::Basename;
    29 use File::Copy;
    30 use PS::IPP::PStamp::RequestFile qw( :standard );
    31 use PS::IPP::PStamp::Job qw( :standard );
    3233use Astro::FITS::CFITSIO qw( :constants );
    3334Astro::FITS::CFITSIO::PerlyUnpacking(1);
    3435
    35 use IPC::Cmd 0.36 qw( can_run run );
    36 
    3736#
    3837# Set up
     
    4140
    4241my $EXTVER = 1.0;
    43 my ($missing_tools,$request_file,$tmp_dir);
    44 my ($req_id,$req_name,$product);
    45 my ($dbname,$need_magic);
    46 my ($input,$output,$workdir,$verbose,$save_temps);
     42my $EXTNAME = 'MOPS_DETECTABILITY_RESPONSE';
     43my ($req_id,$req_name,$product,$need_magic,$missing_tools);
     44my ($request_file,$output,$dbname,$verbose,$save_temps);
    4745GetOptions(
    4846    'input=s'          =>     \$request_file,
    4947    'output=s'        =>      \$output,
    5048    'dbname=s'        =>      \$dbname,
    51     'workdir=s'       =>      \$workdir,
    5249    'verbose'         =>      \$verbose,
    5350    'save-temps'      =>      \$save_temps,
     
    6158
    6259my $detect_query_read = can_run('detect_query_read') or (warn "Can't find detect_query_read" and $missing_tools = 1);
    63 my $psphotForced = can_run('psphotForced') or (warn "Can't find psphotForced" and $missing_tools = 1);
     60my $psphotForced      = can_run('psphotForced') or (warn "Can't find psphotForced" and $missing_tools = 1);
     61my $dquery_finish     = can_run('dquery_finish.pl') or (warn "Can't find dquery_finish.pl" and $missing_tools = 1);
     62my $ppCoord           = can_run('ppCoord') or (warn "Can't find ppCoord" and $missing_tools = 1);
    6463if ($missing_tools) {
    6564    warn("Can't find required tools.");
     
    7170}
    7271
    73 
    7472my $ipprc = PS::IPP::Config->new();
    75 $tmp_dir = "/data/${host}.0/tmp/";
    76 #
    77 # Parse Input
     73#my $tmp_dir = "/data/${host}.0/tmp/";
     74
     75#
     76# Parse input request file using detect_query_read (as it's already written).
    7877#
    7978
     
    8483    warn("Unable to perform $dqr_command error code: $error_code");
    8584}
    86 
    87 my @column_names = ();
    88 my %response = ();
    89 my $section = '';
     85my %query = ();
    9086my $Nrows = 0;
    91 foreach my $entry (split /\n/, (join "", @$stdout_buf)) {
    92     if ($entry =~ /^#/) {
    93         @column_names = split /\s+/, $entry;
    94         shift(@column_names);  # Dump the hash sign.
    95         if ($section eq 'HEADER') {
    96             $section = 'CONTENT';
     87{
     88    my @column_names = ();
     89    my $section = '';
     90    foreach my $entry (split /\n/, (join "", @$stdout_buf)) {
     91        if ($entry =~ /^#/) {
     92            @column_names = split /\s+/, $entry;
     93            shift(@column_names);  # Dump the hash sign.
     94            if ($section eq 'HEADER') {
     95                $section = 'CONTENT';
     96            }
     97            else {
     98                $section = 'HEADER';
     99            }
    97100        }
    98101        else {
    99             $section = 'HEADER';
    100         }
    101     }
    102     else {
    103         # HEADER:
    104         # QUERY_ID FPA_ID MJD_OBS FILTER OBSCODE STAGE
    105         # CONTENT:
    106         # ROWNUM RA1_DEG DEC1_DEG RA2_DEG DEC2_DEG MAG
    107         my @columns = split /\s+/, $entry;
    108         for (my $i = 0; $i <= $#columns; $i++) {
    109             $response{$section}{$column_names[$i]}[$Nrows] = $columns[$i];
    110             print "$section $column_names[$i] $Nrows $columns[$i]\n";
    111         }
    112         $Nrows++;
     102            # HEADER:
     103            # QUERY_ID FPA_ID MJD_OBS FILTER OBSCODE STAGE
     104            # CONTENT:
     105            # ROWNUM RA1_DEG DEC1_DEG RA2_DEG DEC2_DEG MAG
     106            my @columns = split /\s+/, $entry;
     107            for (my $i = 0; $i <= $#columns; $i++) {
     108                $query{$section}{$column_names[$i]}[$Nrows] = $columns[$i];
     109#               print "$section $column_names[$i] $Nrows $columns[$i]\n";
     110            }
     111            $Nrows++;
     112        }
    113113    }
    114114}
     
    119119my %image_list_hash;
    120120for (my $i = 1; $i < $Nrows; $i++) {
    121     print "$i $Nrows $response{CONTENT}{RA1_DEG}[$i] $response{CONTENT}{DEC1_DEG}[$i]\n";
    122 
    123     my $image_set_tmp  = find_image_set($response{HEADER}{FPA_ID}[0],$response{HEADER}{STAGE}[0],
    124                                         $response{HEADER}{MJD_OBS}[0],$response{HEADER}{FILTER}[0],
    125                                         $response{CONTENT}{RA1_DEG}[$i],
    126                                         $response{CONTENT}{DEC1_DEG}[$i],
    127                                         $response{CONTENT}{ROWNUM}[$i],$verbose);
    128     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";
     121#    print "$i $Nrows $query{CONTENT}{RA1_DEG}[$i] $query{CONTENT}{DEC1_DEG}[$i]\n";
     122
     123    my $image_set_tmp  = find_image_set($query{HEADER}{FPA_ID}[0],$query{HEADER}{STAGE}[0],
     124                                        $query{HEADER}{MJD_OBS}[0],$query{HEADER}{FILTER}[0],
     125                                        $query{CONTENT}{RA1_DEG}[$i],$query{CONTENT}{DEC1_DEG}[$i],
     126                                        $query{CONTENT}{ROWNUM}[$i],$verbose);
     127    print "=== $image_set_tmp->{IMAGE}\n    $image_set_tmp->{PSF}\n    $image_set_tmp->{MASK}\n    $image_set_tmp->{WEIGHT}\n    $image_set_tmp->{SKY_COORDINATES}\n    $image_set_tmp->{ROWNUM}\n";
    129128    # This appends, assuming that if we get an image, we also get the identical psf/mask/weight/etc.
    130     $image_list_hash{$image_set_tmp->{IMAGE}}{IMAGE} = $image_set_tmp->{IMAGE};
    131     $image_list_hash{$image_set_tmp->{IMAGE}}{PSF} = $image_set_tmp->{PSF};
    132     $image_list_hash{$image_set_tmp->{IMAGE}}{MASK} = $image_set_tmp->{MASK};
    133     $image_list_hash{$image_set_tmp->{IMAGE}}{WEIGHT} = $image_set_tmp->{WEIGHT};
     129    $image_list_hash{$image_set_tmp->{IMAGE}}{IMAGE}    = $image_set_tmp->{IMAGE};
     130    $image_list_hash{$image_set_tmp->{IMAGE}}{PSF}      = $image_set_tmp->{PSF};
     131    $image_list_hash{$image_set_tmp->{IMAGE}}{MASK}     = $image_set_tmp->{MASK};
     132    $image_list_hash{$image_set_tmp->{IMAGE}}{WEIGHT}   = $image_set_tmp->{WEIGHT};
     133    $image_list_hash{$image_set_tmp->{IMAGE}}{CATALOG}  = $image_set_tmp->{CATALOG};
    134134    $image_list_hash{$image_set_tmp->{IMAGE}}{CLASS_ID} = $image_set_tmp->{CLASS_ID};
    135     $image_list_hash{$image_set_tmp->{IMAGE}}{EXTENSION_BASE} = $image_set_tmp->{EXTENSION_BASE};
    136     push @{ $image_list_hash{$image_set_tmp->{IMAGE}}{COORDINATES} }, $image_set_tmp->{COORDINATES};
     135    push @{ $image_list_hash{$image_set_tmp->{IMAGE}}{SKY_COORDINATES} }, $image_set_tmp->{SKY_COORDINATES};
    137136    push @{ $image_list_hash{$image_set_tmp->{IMAGE}}{ROWNUM} }, $image_set_tmp->{ROWNUM};
    138     $image_list_hash{$image_set_tmp->{IMAGE}}{TARGETS} =
    139         "$tmp_dir/detectability.$response{HEADER}{STAGE}[0].$response{HEADER}{FPA_ID}[0].targets";
    140     $image_list_hash{$image_set_tmp->{IMAGE}}{OUTROOT} =
    141         "$tmp_dir/detectability.$response{HEADER}{STAGE}[0].$response{HEADER}{FPA_ID}[0]";
    142 }
    143 #exit(2);
     137}
    144138my $i = 0;
    145 my @image_list;
     139
    146140foreach my $k (keys %image_list_hash) {
    147     # write coordinates for photometry to a disk file.  Probably need to do this as a tempfile object.
    148     open(T,">$image_list_hash{$k}{TARGETS}") || die "Could not open output TARGET file $image_list_hash{$k}{TARGETS}";
    149     for (my $j = 0; $j <= $#{ $image_list_hash{$k}{COORDINATES} }; $j++) {
    150         print T "$image_list_hash{$k}{COORDINATES}[$j]\n";
    151     }
    152     close(T);
    153     $image_list[$i] = $image_list_hash{$k};
    154     $i++;
    155 }
    156 
    157 #
    158 # Construct and run psphot
    159 #
    160 for (my $i = 0; $i <= $#image_list; $i++) {
    161     my $psf_file    = $image_list[$i]->{PSF};
    162     my $image_file  = $image_list[$i]->{IMAGE};
    163     my $mask_file   = $image_list[$i]->{MASK};
    164     my $weight_file = $image_list[$i]->{WEIGHT};
    165     my $target_file = $image_list[$i]->{TARGETS};
    166     my $out_root    = $image_list[$i]->{OUTROOT};
    167     my $psphot_cmd = "$psphotForced -psf $psf_file ";
    168     $psphot_cmd .= "-file $image_file -mask $mask_file -variance $weight_file ";
    169     $psphot_cmd .= "-srctext $target_file $out_root";
     141    # Write coordinates of the requested targets to a file.
     142    my ($coordfile,$coordname) = tempfile("/tmp/detect.coords.$i.XXXX",
     143                                            UNLINK => !$save_temps);
     144    my ($targetfile,$targetname) = tempfile("/tmp/detect.targets.$i.XXXX",
     145                                            UNLINK => !$save_temps);
     146
     147    for (my $j = 0; $j <= $#{ $image_list_hash{$k}{SKY_COORDINATES} }; $j++) {
     148        print $coordfile "$image_list_hash{$k}{SKY_COORDINATES}[$j]\n";
     149    }
     150
     151    # Convert the sky coordinates to image coordinates with ppCoord.
     152    my $command = "ppCoord -astrom $image_list_hash{$k}{CATALOG} -radec $coordname";
    170153    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
     154        run(command => $command, verbose => $verbose);
     155    unless ($success) {
     156        warn("Unable to perform $command error code: $error_code");
     157    }
     158    my @response = split /\n/, (join "", @$stdout_buf);
     159    foreach my $line (@response) {
     160        my ($r_ra,$r_dec,$trash,$r_x,$r_y,$r_chip) = split /\s+/, $line;
     161        print $targetfile "$r_x $r_y\n";
     162        $image_list_hash{$k}{EXTENSION_BASE} = $r_chip;
     163    }
     164
     165   
     166    # Run psphotForced on the target list.
     167    my $tmpdir  = tempdir("detect.$i.XXXX", DIR => "/tmp/", CLEANUP => !$save_temps);
     168    $image_list_hash{$k}{OUTROOT} = "$tmpdir/detectability.$query{HEADER}{STAGE}[0].$query{HEADER}{FPA_ID}[0]";
     169   
     170    my $psphot_cmd = "$psphotForced -psf $image_list_hash{$k}{PSF} ";
     171    $psphot_cmd .= "-file $image_list_hash{$k}{IMAGE} ";
     172    $psphot_cmd .= "-mask $image_list_hash{$k}{MASK} ";
     173    $psphot_cmd .= "-variance $image_list_hash{$k}{WEIGHT} ";
     174    $psphot_cmd .= "-srctext $targetname $image_list_hash{$k}{OUTROOT}";
     175
     176    ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
    171177        run(command => $psphot_cmd, verbose => $verbose);
    172178    unless ($success) {
     
    178184# Convert psphot output to response
    179185#
    180 my @response_entries = ();
    181 
     186my @rownums = ();
     187my @psphot_Npix = ();
     188my @psphot_Qfact= ();
     189my @psphot_flux = ();
    182190foreach my $k (keys %image_list_hash) {
    183     my $class_id = $image_list_hash{$k}{CLASS_ID};
    184     my $cmf = $image_list_hash{$k}{OUTROOT} . ".${class_id}.cmf";
    185     my ($detect_N_col,$detect_F_col,@detectability_data) = read_cmf_file($cmf,$image_list_hash{$k}{EXTENSION_BASE});
    186    
    187     if ($#detectability_data != $#{ $image_list_hash{$k}{ROWNUM} }) {
    188         my_die("Number of measured values does not match number of requested values.");
    189     }
    190     for (my $i = 0; $i < $#detectability_data; $i++) {
    191         push @response_entries, create_response_entry($image_list_hash{$k}{ROWNUM}[$i],
    192                                                       $detect_N_col,$detect_F_col,
    193                                                       $detectability_data[$i]);
    194     }
     191    my $cmf = "$image_list_hash{$k}{OUTROOT}.$image_list_hash{$k}{CLASS_ID}.cmf";
     192   
     193    my ($tmp_Npix,$tmp_Qfact,$tmp_flux) = read_cmf_file($cmf,$image_list_hash{$k}{EXTENSION_BASE});
     194
     195    push @rownums,        @{ $image_list_hash{$k}{ROWNUM} };
     196    push @psphot_Npix,    @{ $tmp_Npix };
     197    push @psphot_Qfact,   @{ $tmp_Qfact };
     198    push @psphot_flux,    @{ $tmp_flux };
    195199}
    196200
    197201write_response_file($output,
    198                     $response{HEADER}{QUERY_ID}[0],$response{HEADER}{FPA_ID}[0],
    199                     $response{HEADER}{MJD_OBS}[0],$response{HEADER}{filter}[0],
    200                     $response{HEADER}{obscode}[0],
    201                     @response_entries);
     202                    $query{HEADER}{QUERY_ID}[0],$query{HEADER}{FPA_ID}[0],
     203                    $query{HEADER}{MJD_OBS}[0],$query{HEADER}{filter}[0],
     204                    $query{HEADER}{obscode}[0],
     205                    \@rownums, \@psphot_Npix, \@psphot_Qfact, \@psphot_flux);
    202206
    203207#
     
    214218# Cleanup
    215219#
     220# Since everything is written to temporary files, there should be nothing to cleanup.
    216221#
    217222# Utilities
     
    236241    my $mjd_max = $mjd + 1;
    237242
     243    # Call the PStamp code to find the images that contain the target on the given MJD in the specified filter.
    238244    my @images = locate_images($ipprc,$dbname,"bycoord",$stage,
    239245                               undef,undef,undef,$option_mask,$need_magic,
    240246                               $ra,$dec,$mjd_min,$mjd_max,$filter . ".00000",undef,$verbose); 
    241 #     my @images = lookup($ipprc,$dbname,$req_type,$img_type,
    242 #                       $id,$tess_id,$componenent,$need_magic,$dateobs_begin,$dateobs_end,
    243 #                       $filter,$data_group,$option_mask,$verbose);
    244     my @keys_i_care = ('image','mask','psf','exp_name','stage','weight','cmf','smf');
    245     foreach my $i (@images) {
    246         foreach my $j (@{ $i }) {
    247 #           foreach my $k (keys %{ ${ $i }[$j] }) {
    248         }
    249     }
    250247
    251248    my %image_info  = ();
    252 
    253     foreach my $i (@images) {
    254         foreach my $j (@{ $i }) {
     249    foreach my $i (@images) {           # Scan over each result
     250        foreach my $j (@{ $i }) {       # Scan over each image in the result
     251            # We only care about an image if it matches the FPA_ID in the request
    255252            if ($stage eq 'diff') {
    256                 #Diffs hold both exposure names, if they are defined (so a WW or WS diff). 
     253                # Diffs match if either exposure name is defined and matches the FPA_ID
    257254                unless ((defined(${ $j }{exp_name_1}) && (${ $j }{exp_name_1} eq $FPA_ID))||
    258255                        (defined(${ $j }{exp_name_2}) && (${ $j }{exp_name_2} eq $FPA_ID))) {
    259256                    next;
    260257                }
    261                 push @keys_i_care, 'exp_name_1','exp_name_2';
    262258            }
    263259            elsif ($stage eq 'stack') {
    264                 # Stacks do not store the exposure name information...I'll have to think up something smart here.
    265                 next;
     260                # Stacks hide the exposure name very well, so I'm
     261                # choosing the stack_id as the FPA_ID to match. This
     262                # probably needs to be approved somehow.
     263                if (${ $j }{stage_id} ne $FPA_ID) {
     264                    next;
     265                }
    266266            }
    267267            else {
     268                # For all the other stages (warp and chip are the ones I've tested), we can simply
     269                # directly match the exposure name to the FPA_ID
    268270                if (${ $j }{exp_name} ne $FPA_ID) {
    269271                    next;
    270272                }
    271273            }
    272 #           foreach my $k (@keys_i_care) {
    273             foreach my $k (keys %{ $j }) {
    274                 print "$i $j $k ${ $j }{$k}\n";
    275             }
    276 
    277             my $image_name = ${ $j }{image};
    278             $image_info{IMAGE} = $image_name;
    279             $image_info{PSF} = ${ $j }{psf};
    280             $image_info{MASK} = ${ $j }{mask};
     274            # Debug prints of all the components of this image
     275#           foreach my $k (keys %{ $j }) {
     276#               print "$i $j $k ${ $j }{$k}\n";
     277#           }
     278           
     279            # This image matches, so we want to save the information into our output structure
     280            $image_info{ROWNUM} = $index;
     281            $image_info{IMAGE}  = ${ $j }{image};
     282            $image_info{PSF}    = ${ $j }{psf};
     283            $image_info{MASK}   = ${ $j }{mask};
    281284            $image_info{WEIGHT} = ${ $j }{weight};
     285            $image_info{SKY_COORDINATES} = "$ra $dec";
     286            # To do sky->image coordinate transformations, we need to use the cmf/smf file. If
     287            # an astrom reference (the camera stage smf file) exists, then use that, as we're dealing with
     288            # with the chip stage. Otherwise, use the stage-dependent cmf (and set the class_id to fpa).
     289            # The EXTENSION_BASE stores the basename of the extension that will be generated by psphotForced.
    282290            if (exists(${ $j }{astrom})) {
    283                 ($image_info{COORDINATES},$image_info{EXTENSION_BASE}) =
    284                     calculate_image_coordinates(${ $j }{astrom},$ra,$dec);
     291                $image_info{CATALOG} = ${ $j }{astrom};
    285292                $image_info{CLASS_ID} = ${ $j }{class_id};
     293
    286294            }
    287295            else {
    288                 ($image_info{COORDINATES},$image_info{EXTENSION_BASE}) =
    289                     calculate_image_coordinates(${ $j }{cmf},$ra,$dec);
     296                $image_info{CATALOG} = ${ $j }{cmf};
    290297                $image_info{CLASS_ID} = 'fpa';
    291             }
    292             $image_info{ROWNUM} = $index;
     298
     299            }
     300
    293301        }
    294302    }
    295303    return(\%image_info);
    296304}
     305
    297306# Taken largely from detect_query_read
    298 sub calculate_image_coordinates {
    299     my $cmf = shift;
    300     my $ra = shift;
    301     my $dec = shift;
    302     my $x = $ra;
    303     my $y = $dec;
    304     my $c = '';
    305    
    306     my $tempfile = "$tmp_dir/detect.cic.$$";
    307     open(TEMP,">$tempfile");
    308     print TEMP "$ra $dec\n";
    309     close(TEMP);
    310     my $command = "ppCoord -astrom $cmf -radec $tempfile";
    311    
    312     my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
    313         run(command => $command, verbose => $verbose);
    314     unless ($success) {
    315         warn("Unable to perform $command error code: $error_code");
    316     }
    317     print ">>>$command<<<\n";
    318     my @response = split /\n/, (join "", @$stdout_buf);
    319     foreach my $line (@response) {
    320         print ">>$line<<\n";
    321         my ($r_ra,$r_dec,$trash,$r_x,$r_y,$r_chip) = split /\s+/, $line;
    322         if (($r_ra == $ra)&&($r_dec == $dec)) {
    323             $x = $r_x;
    324             $y = $r_y;
    325             $c = $r_chip;
    326         }
    327     }
    328     unlink($tempfile);
    329     return("$x $y",$c);
    330 }
    331307sub read_cmf_file {
    332308    my $cmf_file = shift;
    333309    my $class_id = shift;
    334     my @outdata = ();
    335310    my $extname = $class_id . ".psf";
    336311
    337     my ($detect_F_col,$detect_N_col,$detect_flux_col);
    338     print "RCF: $cmf_file $class_id $extname\n";
     312    my @tmp_Npix = ();
     313    my @tmp_Qfact = ();
     314    my @tmp_flux = ();
     315
     316    my ($detect_F_col,$detect_N_col,$detect_flux_col,$detect_mag_col) = (-1, -1, -1, -1);
     317
     318   
    339319    my $status = 0;
    340320    my $inFits = Astro::FITS::CFITSIO::open_file( $cmf_file, READONLY, $status ); # Open CMF
    341321    check_fitsio($status);
    342 
    343322    $inFits->movnam_hdu(BINARY_TBL, $extname, 0, $status) and check_fitsio($status);
    344323
     324    my $inHeader = $inFits->read_header();
     325
     326    my $CMFversion = $inHeader->{EXTTYPE};
     327    $CMFversion =~ s/\'//g;
     328    $CMFversion =~ s/\s+//g;
     329
    345330    # This is the only data we actually care about
    346     my $column_defs = [
    347         { name => 'PSF_QF',   type => '1E', writetype => TDOUBLE },
    348         { name => 'PSF_NPIX', type => '1J', writetype => TLONG },
    349         { name => 'PSF_INST_MAG', type => '1E', writetype => TDOUBLE },
    350         ];
     331    my $column_defs;
     332    if ($CMFversion eq 'PS1_V2') {
     333        $column_defs = [
     334            { name => 'PSF_QF',   type => '1E', writetype => TDOUBLE },
     335            { name => 'PSF_NPIX', type => '1J', writetype => TLONG },
     336            { name => 'PSF_INST_MAG', type => '1E', writetype => TDOUBLE },
     337            ];
     338    }
     339    elsif ($CMFversion eq 'PS1_DV2') {
     340        $column_defs = [
     341            { name => 'PSF_QF',   type => '1E', writetype => TDOUBLE },
     342            { name => 'PSF_NPIX', type => '1J', writetype => TLONG },
     343            { name => 'PSF_INST_FLUX', type => '1E', writetype => TDOUBLE },
     344            ];
     345    }
     346
    351347    my @colNames;
    352348    my @colTypes;
     
    364360        my ($col_num,$col_type,$col_data);
    365361        $inFits->get_colnum(0, $col->{name}, $col_num, $status) and check_fitsio($status);
     362        if ($status != 0) {
     363            $status = 0;
     364            next;
     365        }
    366366        $inFits->get_coltype($col_num, $col_type, undef, undef, $status) and check_fitsio($status);
    367367        $inFits->read_col($col_type, $col_num, 1, 1, $numRows, 0, $col_data, undef, $status) and check_fitsio($status);
    368         $colData{$col->{name}} = $col_data;
     368       
    369369        if ($col->{name} eq 'PSF_QF') {
    370             $detect_F_col = $col_num;
    371         }
    372         if ($col->{name} eq 'PSF_NPIX') {
    373             $detect_N_col = $col_num;
    374         }
    375         if ($col->{name} eq 'PSF_INST_MAG') {
    376             $detect_flux_col = $col_num;
     370            @tmp_Qfact = @{ $col_data };
     371        }
     372        elsif ($col->{name} eq 'PSF_NPIX') {
     373            @tmp_Npix = @{ $col_data };
     374        }
     375        elsif ($col->{name} eq 'PSF_INST_MAG') {
     376            @tmp_flux = map { $_ = 10**(-0.4 * $_) } @{ $col_data };
     377        }
     378        elsif ($col->{name} eq 'PSF_INST_FLUX') {
     379            @tmp_flux = @{ $col_data };
    377380        }
    378381    }
    379382    $inFits->close_file( $status ) and check_fitsio($status);   
    380 
    381     for (my $i = 0; $i < $numRows; $i++) {
    382         for (my $j = 0; $j <= $#colNames; $j++) {
    383             if ($j == $detect_flux_col) {
    384                 $outdata[$i][$j] = 10**(-0.4 * $colData{$colNames[$j]}->[$i]);
    385             }
    386             else {
    387                 $outdata[$i][$j] = $colData{$colNames[$j]}->[$i];
    388             }
    389         }
    390     }
    391 
    392     return($detect_N_col,$detect_F_col,@outdata);
    393 }
    394 
    395 sub create_response_entry {
    396     my $rownum = shift;
    397     my $detect_N_col = shift;
    398     my $detect_F_col = shift;
    399 
    400     my @cmf_data = @_;
    401    
    402     my $entry;
    403     @{ $entry } = ();
    404    
    405    
    406     push @{ $entry }, $rownum;
    407     push @{ $entry }, $cmf_data[$detect_N_col];
    408     push @{ $entry }, $cmf_data[$detect_F_col];
    409    
    410     return($entry);
    411 }
     383    return(\@tmp_Npix, \@tmp_Qfact, \@tmp_flux);
     384}
     385
    412386# A lot of this pulled verbatim from Bill's detect_response_create file
    413387sub write_response_file {
     
    418392    my $filter = shift;
    419393    my $obscode = shift;
    420     my @response_entries = @_;
    421 #    my ($query_id,$FPA_ID,$MJD_OBS,$filter,$obscode);
     394    my $rownum_ref = shift;
     395    my $psphot_Npix_ref = shift;
     396    my $psphot_Qfact_ref = shift;
     397    my $psphot_flux_ref = shift;
    422398    my $status = 0;
     399
    423400    # Specification of columns to write
    424401    my $columns = [
     
    429406        # detectibility, indicating the fraction of PSF pixels detetable by IPP
    430407        { name => 'DETECT_F', type => 'D',   writetype => TDOUBLE },
     408        # flux of the target source
     409        { name => 'TARGET_FLUX', type => 'D', writetype => TDOUBLE },
    431410        ];
    432411   
     
    451430    }
    452431
    453     my $numRows = $#response_entries + 1;
     432    my $numRows = $#{ $rownum_ref } + 1;
    454433    my $inHeader = { };
    455434
    456     # Hack to force the data to match teh detect_response_create formats
     435    # Hack to force the data to match the detect_response_create formats
    457436    $inHeader->{QUERY_ID}->{value} = $query_id;
    458437    $inHeader->{FPA_ID}->{value} = $FPA_ID;
     
    461440    $inHeader->{OBSCODE}->{value} = $obscode;
    462441   
    463     for (my $i = 0; $i <= $#response_entries; $i++) {
    464         push @{$colData{'ROWNUM'}}, $response_entries[$i][0];
    465         push @{$colData{'DETECT_N'}}, $response_entries[$i][1];
    466         push @{$colData{'DETECT_F'}}, $response_entries[$i][2];
     442    # Fill the table columns with the data, making sure the flux is defined
     443    for (my $i = 0; $i < $numRows; $i++) {
     444        push @{$colData{'ROWNUM'}},      ${ $rownum_ref }[$i];
     445        push @{$colData{'DETECT_N'}},    ${ $psphot_Npix_ref }[$i];
     446        push @{$colData{'DETECT_F'}},    ${ $psphot_Qfact_ref }[$i];
     447        push @{$colData{'TARGET_FLUX'}}, ${ $psphot_flux_ref }[$i];
     448        unless (defined(${ $colData{'TARGET_FLUX'}}[-1])) {
     449            $colData{'TARGET_FLUX'}[-1] = 0.0;
     450        }
    467451    }
    468452
     
    470454    my $outFits = Astro::FITS::CFITSIO::create_file( $output, $status );
    471455    check_fitsio( $status );
    472 
    473456    $outFits->create_img( 16, 0, undef, $status );
    474457    check_fitsio( $status );
     
    498481        my $colName = $colNames[$i];# Column name
    499482        my $writeType = $colWriteType[$i];
    500         print "$colName $writeType $i $colData{$colName}[$i] $status $numRows\n";
    501483        $outFits->write_col( $writeType, $i + 1, 1, 1, $numRows, $colData{$colName}, $status );
    502         print "$colName $writeType $i $colData{$colName}[$i] $status $numRows\n";
    503484        check_fitsio( $status );
    504485    }
Note: See TracChangeset for help on using the changeset viewer.