IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 10, 2020, 9:40:15 AM (6 years ago)
Author:
eugene
Message:

gDiff now runs and does not crash on forced photometry

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/gDiff.20200202/src/gDiffLoop.c

    r41260 r41413  
     1/** @file gDiffLoop.c
     2 *  @brief: image subtraction based on matched Gaussians profiles
     3 *  @author Eugene Magnier @ IfA
     4 *  @version $Revision: 1.12 $
     5 *  @date $Date: 2020-01-25 $
     6 *  Copyright 2020 Institute for Astronomy, University of Hawaii
     7 */
     8
    19#include "gDiff.h"
     10
     11# define HDUCHECK(A)                                                    \
     12    fprintf (stderr, "%s\n", A);                                        \
     13    listhdu (config, "GDIFF.OUTPUT"); listhdu (config, "PSPHOT.INPUT");
     14
     15// internal functions:
     16bool dumpout(pmConfig *config, char *name);
     17bool listhdu(pmConfig *config, char *name);
     18
     19bool gDiffLoop(gDiffData *data)
     20{
     21    bool mdok = false;
     22    bool success = true;
     23
     24    psAssert(data, "Require processing data");
     25    pmConfig *config = data->config;    // Configuration
     26    psAssert(config, "Require configuration.");
     27
     28    HDUCHECK ("0a: start gDiffLoop");
     29
     30    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim
     31    psAssert(recipe, "We checked this earlier, so it should be here.");
     32    bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images.
     33
     34    pmConfigCamerasCull(config, NULL);
     35    pmConfigRecipesCull(config, "GDIFF,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
     36
     37    pmFPAfile *input     = psMetadataLookupPtr(NULL, config->files, "GDIFF.INPUT");
     38    pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "GDIFF.REF");
     39    pmFPAfile *output    = psMetadataLookupPtr(NULL, config->files, "GDIFF.OUTPUT");
     40    psAssert(input && reference && output, "Require files");
     41
     42    // load input files : we are required to load a complete set of files that describe the input.
     43    // we are not able to process subdivisions of an image (e.g., cells vs chips)
     44    if (!gDiffFilesIterateDown(config, GDIFF_FILES_INPUT | GDIFF_FILES_CONV)) {
     45        psError(GDIFF_ERR_IO, false, "Unable to load files.");
     46        return false;
     47    }
     48
     49    psTimerStart("GDIFF_MATCH");
     50
     51    HDUCHECK ("0b: gDiffSetMasks");
     52
     53    if (!gDiffSetMasks(config)) {
     54        psError(psErrorCodeLast(), false, "Unable to set masks.");
     55        return false;
     56    }
     57
     58    HDUCHECK ("0b.1: gDiffInputDetections");
     59    if (data->forcedPhot1) {
     60        bool foundDetections = false;
     61        if (!gDiffInputDetections(&foundDetections, "GDIFF.POS1.SOURCES", "GDIFF.INPUT", data)) {
     62            psError(psErrorCodeLast(), false, "Unable to measure positive detections (1)");
     63            success = false;
     64            goto ESCAPE;
     65        }
     66        // if nothing was found, don't bother doing the forced photometry below
     67        if (!foundDetections) {
     68            psWarning ("no sources found in positive image 1, skipping forced photometry");
     69            data->forcedPhot1 = false;
     70        }
     71    }
     72    HDUCHECK ("0b.2: gDiffInputDetections");
     73    if (data->forcedPhot2) {
     74        // XXX why is this only done for Phot2 and not Phot1 ??
     75
     76        // Change the recipe to use a higher nsigma limit and quit after pass1
     77        psMetadata *psphotRecipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
     78
     79        psF32 nsigma_peak_save = psMetadataLookupF32 (NULL, psphotRecipe, "PEAKS_NSIGMA_LIMIT");
     80        char *breakPt_save =  psMetadataLookupStr (NULL, psphotRecipe, "BREAK_POINT");
     81
     82        psF32 pos2_nsigma_peak = psMetadataLookupF32 (&mdok, psphotRecipe, "PEAKS_POS2_NSIGMA_LIMIT");
     83        if (!mdok) {
     84            psWarning("PEAKS_POS2_NSIGMA_LIMIT not found in psphotRecipe. Will use 25.\n");
     85            pos2_nsigma_peak = 25.;
     86        }
     87        psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", pos2_nsigma_peak);
     88        psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", "PASS1");
     89
     90        bool foundDetections = false;
     91        if (!gDiffInputDetections(&foundDetections, "GDIFF.POS2.SOURCES", "GDIFF.REF", data)) {
     92            psError(psErrorCodeLast(), false, "Unable to measure positive detections (2)");
     93            psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save);
     94            psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save);
     95            success = false;
     96            goto ESCAPE;
     97        }
     98        psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save);
     99        psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save);
     100        // if nothing was found, don't bother doing the forced photometry below
     101        if (!foundDetections) {
     102            psWarning ("no sources found in positive image 2, skipping forced photometry");
     103            data->forcedPhot2 = false;
     104        }
     105    }
     106
     107    // XXX if it exists, use the POS1, POS2 successs for the FWHMs
     108
     109    /**** The PSF matching process is performed in this function below. ****
     110     **** This is the core of the PSF-matched image subtraction code.   ****/
     111
     112    HDUCHECK ("1: before gDiffMatchPSFs");
     113
     114    if (!gDiffMatchPSFs(data)) {
     115        psError(psErrorCodeLast(), false, "Unable to match PSFs.");
     116        success = false;
     117        goto ESCAPE;
     118    }
     119    if (data->quality) {
     120        // Can't do anything at all
     121        success = false;
     122        goto ESCAPE;
     123    }
     124
     125    // the functions in this block are not critical for initial testing, so block them out for now
     126
     127    // generate the residual stamp grid for visualization
     128    if (!gDiffResidualSampleJpeg(config)) {
     129        psError(psErrorCodeLast(), false, "Unable to update.");
     130        success = false;
     131        goto ESCAPE;
     132    }
     133
     134    // XXX add in a positive image detection step here (if needed)
     135    // XXX EAM what does this mean???
     136
     137    psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", psTimerClear("GDIFF_MATCH"));
     138
     139    // Close input files (freeing up space) : for the noConvolve case, we cannot close the
     140    // inputs since we will use them for subtraction below.  wait until later to do the work
     141    if (!noConvolve) {
     142        if (!gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) {
     143            psError(GDIFF_ERR_IO, false, "Unable to close input files.");
     144            success = false;
     145            goto ESCAPE;
     146        }
     147    }
     148    if (!gDiffLowThreshold(data)) {
     149        psError(psErrorCodeLast(), false, "Unable to threshold images.");
     150        success = false;
     151        goto ESCAPE;
     152    }
     153
     154    HDUCHECK ("2: after gDiffLowThreshold");
     155
     156    // Set up subtraction files
     157    if (!gDiffFilesIterateDown(config, GDIFF_FILES_SUB)) {
     158        psError(GDIFF_ERR_IO, false, "Unable to set up subtraction files.");
     159        success = false;
     160        goto ESCAPE;
     161    }
     162
     163    if (!gDiffDefineOutput("GDIFF.OUTPUT", config)) {
     164        psError(psErrorCodeLast(), false, "Unable to define output.");
     165        success = false;
     166        goto ESCAPE;
     167    }
     168
     169    HDUCHECK ("2a: gDiffMakePSF");
     170
     171    // the functions in this block are not critical for initial testing, so block them out for now
     172    if (!data->quality && !gDiffMakePSF(data)) {
     173        psError(psErrorCodeLast(), false, "Unable to generate PSF.");
     174        success = false;
     175        goto ESCAPE;
     176    }
     177
     178    // Now we've got a PSF, blow away detections in case they're confused with real output detections
     179    {
     180        pmFPAview *view = gDiffViewReadout(); // View to readout
     181        pmReadout *out = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT");
     182        if (psMetadataLookup(out->analysis, "PSPHOT.DETECTIONS")) {
     183            psMetadataRemoveKey(out->analysis, "PSPHOT.DETECTIONS");
     184        }
     185        psFree(view);
     186    }
     187
     188    HDUCHECK ("3: before gDiffReadoutSubtract");
     189
     190    if (!gDiffReadoutSubtract(config)) {
     191        psError(psErrorCodeLast(), false, "Unable to subtract images.");
     192        success = false;
     193        goto ESCAPE;
     194    }
     195
     196    // Close convolved files
     197    if (!gDiffFilesIterateUp(config, GDIFF_FILES_PSF | GDIFF_FILES_CONV)) {
     198        psError(GDIFF_ERR_IO, false, "Unable to close input files.");
     199        success = false;
     200        goto ESCAPE;
     201    }
     202    // Close unconvolved input files if we did the noConvolve subtraction
     203    // (otherwise they are closed above)
     204    if (noConvolve && !gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) {
     205        psError(GDIFF_ERR_IO, false, "Unable to close input files.");
     206        success = false;
     207        goto ESCAPE;
     208    }
     209   
     210    HDUCHECK ("4: before gDiffBackground");
     211
     212    // Higher order background subtraction using psphot
     213    if (!gDiffBackground(config)) {
     214        psError(psErrorCodeLast(), false, "Unable to subtract background.");
     215        success = false;
     216        goto ESCAPE;
     217    }
     218    // dumpout(config, "diff.2b.fits");
     219
     220    HDUCHECK ("5: before gDiffVarianceRescale");
     221
     222    // Perform Variance correction (rescale within a modest range)
     223    if (!gDiffVarianceRescale(config, data)) {
     224        psError(psErrorCodeLast(), false, "Unable to rescale variance.");
     225        success = false;
     226        goto ESCAPE;
     227    }
     228    // dumpout(config, "diff.2c.fits");
     229
     230    if (data->quality) {
     231        // Done all we can do up to this point
     232        success = false;
     233        goto ESCAPE;
     234    }
     235
     236    if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_SUB)) {
     237        psError(GDIFF_ERR_IO, false, "Unable to set up photometry files.");
     238        success = false;
     239        goto ESCAPE;
     240    }
     241
     242    HDUCHECK ("6: before gDiffReadoutPhotometry");
     243
     244    if (!data->quality && !gDiffReadoutPhotometry("GDIFF.OUTPUT", data)) {
     245        psError(psErrorCodeLast(), false, "Unable to perform photometry.");
     246        success = false;
     247        goto ESCAPE;
     248    }
     249    // dumpout(config, "diff.3.fits");
     250
     251    HDUCHECK ("7: before gDiffReadoutForcedPhot (1)");
     252
     253    // forced photometry for positive image 1
     254    if (data->forcedPhot1 && !data->quality) {
     255        if (!gDiffReadoutForcedPhot("GDIFF.FORCED1.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS1.SOURCES", data)) {
     256            psError(psErrorCodeLast(), false, "Unable to perform photometry.");
     257            success = false;
     258            goto ESCAPE;
     259        }
     260    }
     261
     262    // forced photometry for positive image 2
     263    if (data->forcedPhot2 && !data->quality) {
     264        if (!gDiffReadoutForcedPhot("GDIFF.FORCED2.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS2.SOURCES", data)) {
     265            psError(psErrorCodeLast(), false, "Unable to perform photometry.");
     266            success = false;
     267            goto ESCAPE;
     268        }
     269    }
     270
     271    if (!gDiffFilesIterateUp(config, GDIFF_FILES_PHOT_SUB)) {
     272        psError(GDIFF_ERR_IO, false, "Unable to set up photometry files.");
     273        success = false;
     274        goto ESCAPE;
     275    }
     276
     277    // Perform statistics on the cell
     278    if (!gDiffReadoutStats(data)) {
     279        psError(psErrorCodeLast(), false, "Unable to collect statistics");
     280        success = false;
     281        goto ESCAPE;
     282    }
     283    // Do Mask Stats
     284    {
     285        pmFPAview *view = gDiffViewReadout(); // View to readout
     286        if (!gDiffMaskStats(config, view,data->stats)) {
     287            psError(psErrorCodeLast(), false, "Unable to generate mask statistics");
     288            success = false;
     289            psFree(view);
     290            goto ESCAPE;
     291        }
     292        psFree(view);
     293    }
     294    // dumpout(config, "diff.4.fits");
     295   
     296    // generate the binned image used to write the jpeg
     297    if (!gDiffReadoutJpeg(config)) {
     298        psError(psErrorCodeLast(), false, "Unable to update.");
     299        success = false;
     300        goto ESCAPE;
     301    }
     302
     303    if (data->inverse) {
     304        // Set up inverse subtraction files
     305        if (!gDiffFilesIterateDown(config, GDIFF_FILES_INV)) {
     306            psError(GDIFF_ERR_IO, false, "Unable to set up inverse files.");
     307            success = false;
     308            goto ESCAPE;
     309        }
     310
     311        if (data->inverse && !gDiffDefineOutput("GDIFF.INVERSE", config)) {
     312            psError(psErrorCodeLast(), false, "Unable to define inverse.");
     313            success = false;
     314            goto ESCAPE;
     315        }
     316
     317        if (!gDiffReadoutInverse(config)) {
     318            psError(psErrorCodeLast(), false, "Unable to invert images.");
     319            success = false;
     320            goto ESCAPE;
     321        }
     322
     323        //MEH -- need to also clear out detections or inv.cmf corrupted on update
     324        pmFPAview *viewinv = gDiffViewReadout(); // View to readout
     325        pmReadout *inv = pmFPAfileThisReadout(config->files, viewinv, "GDIFF.INVERSE");
     326        if (psMetadataLookup(inv->analysis, "PSPHOT.DETECTIONS")) {
     327            psMetadataRemoveKey(inv->analysis, "PSPHOT.DETECTIONS");
     328        }
     329        psFree(viewinv);
     330
     331        // Close subtraction files and open inverse photometry files
     332        if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) {
     333            psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
     334            success = false;
     335            goto ESCAPE;
     336        }
     337
     338        if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_INV)) {
     339            psError(GDIFF_ERR_IO, false, "Unable to set up inverse files.");
     340            success = false;
     341            goto ESCAPE;
     342        }
     343
     344        if (!data->quality && !gDiffReadoutPhotometry("GDIFF.INVERSE", data)) {
     345            psError(psErrorCodeLast(), false, "Unable to perform photometry.");
     346            success = false;
     347            goto ESCAPE;
     348        }
     349
     350        // Close inverse subtraction files
     351        if (!gDiffFilesIterateUp(config, GDIFF_FILES_INV | GDIFF_FILES_PHOT_INV)) {
     352            psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
     353            success = false;
     354            goto ESCAPE;
     355        }
     356    } else {
     357        // Close subtraction files
     358        // dumpout(config, "diff.5.fits");
     359        if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) {
     360            psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
     361            success = false;
     362            goto ESCAPE;
     363        }
     364    }
     365
     366 ESCAPE:
     367    pmFPAfileDropInternal(config->files, "PSPHOT.BACKGND");
     368    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
     369    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
     370
     371    return success;
     372}
     373
     374/**** internal functions ****/
    2375
    3376bool dumpout(pmConfig *config, char *name)
    … …  
    10383}
    11384
    12 
    13 bool gDiffLoop(gDiffData *data)
     385bool listhdu(pmConfig *config, char *name)
    14386{
    15     bool mdok = false;
    16     bool success = true;
    17 
    18     psAssert(data, "Require processing data");
    19     pmConfig *config = data->config;    // Configuration
    20     psAssert(config, "Require configuration.");
    21 
    22     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim
    23     psAssert(recipe, "We checked this earlier, so it should be here.");
    24     bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images.
    25 
    26     pmConfigCamerasCull(config, NULL);
    27     pmConfigRecipesCull(config, "GDIFF,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
    28 
    29     pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "GDIFF.INPUT");
    30     pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "GDIFF.REF");
    31     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "GDIFF.OUTPUT");
    32     psAssert(input && reference && output, "Require files");
    33 
    34     if (!gDiffFilesIterateDown(config, GDIFF_FILES_INPUT | GDIFF_FILES_CONV)) {
    35         psError(GDIFF_ERR_IO, false, "Unable to load files.");
    36         return false;
    37     }
    38 
    39     psTimerStart("GDIFF_MATCH");
    40 
    41     if (!gDiffSetMasks(config)) {
    42         psError(psErrorCodeLast(), false, "Unable to set masks.");
    43         return false;
    44     }
    45 
    46     if (data->forcedPhot1) {
    47         bool foundDetections = false;
    48         if (!gDiffInputDetections(&foundDetections, "GDIFF.POS1.SOURCES", "GDIFF.INPUT", data)) {
    49             psError(psErrorCodeLast(), false, "Unable to measure positive detections (1)");
    50             success = false;
    51             goto ESCAPE;
    52         }
    53         // if nothing was found, don't bother doing the forced photometry below
    54         if (!foundDetections) {
    55             psWarning ("no sources found in positive image 1, skipping forced photometry");
    56             data->forcedPhot1 = false;
    57         }
    58     }
    59     if (data->forcedPhot2) {
    60         // Change the recipe to use a higher nsigma limit and quit after pass1
    61         psMetadata *psphotRecipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
    62 
    63         psF32 nsigma_peak_save = psMetadataLookupF32 (NULL, psphotRecipe, "PEAKS_NSIGMA_LIMIT");
    64         char *breakPt_save =  psMetadataLookupStr (NULL, psphotRecipe, "BREAK_POINT");
    65 
    66         psF32 pos2_nsigma_peak = psMetadataLookupF32 (&mdok, psphotRecipe, "PEAKS_POS2_NSIGMA_LIMIT");
    67         if (!mdok) {
    68             psWarning("PEAKS_POS2_NSIGMA_LIMIT not found in psphotRecipe. Will use 25.\n");
    69             pos2_nsigma_peak = 25.;
    70         }
    71         psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", pos2_nsigma_peak);
    72         psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", "PASS1");
    73 
    74         bool foundDetections = false;
    75         if (!gDiffInputDetections(&foundDetections, "GDIFF.POS2.SOURCES", "GDIFF.REF", data)) {
    76             psError(psErrorCodeLast(), false, "Unable to measure positive detections (2)");
    77             psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save);
    78             psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save);
    79             success = false;
    80             goto ESCAPE;
    81         }
    82         psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save);
    83         psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save);
    84         // if nothing was found, don't bother doing the forced photometry below
    85         if (!foundDetections) {
    86             psWarning ("no sources found in positive image 2, skipping forced photometry");
    87             data->forcedPhot2 = false;
    88         }
    89     }
    90 
    91     // XXX if it exists, use the POS1, POS2 successs for the FWHMs
    92     if (!gDiffMatchPSFs(data)) {
    93         psError(psErrorCodeLast(), false, "Unable to match PSFs.");
    94         success = false;
    95         goto ESCAPE;
    96     }
    97 
    98     if (data->quality) {
    99         // Can't do anything at all
    100         success = false;
    101         goto ESCAPE;
    102     }
    103     // generate the residual stamp grid for visualization
    104     if (!gDiffResidualSampleJpeg(config)) {
    105         psError(psErrorCodeLast(), false, "Unable to update.");
    106         success = false;
    107         goto ESCAPE;
    108     }
    109 
    110     // XXX add in a positive image detection step here (if needed)
    111    
    112 
    113     psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs",
    114                      psTimerClear("GDIFF_MATCH"));
    115 
    116     // Close input files (freeing up space) : for the noConvolve case, we cannot close the
    117     // inputs since we will use them for subtraction below.  wait until later to do the work
    118     if (!noConvolve) {
    119         if (!gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) {
    120             psError(GDIFF_ERR_IO, false, "Unable to close input files.");
    121             success = false;
    122             goto ESCAPE;
    123         }
    124     }
    125     if (!gDiffLowThreshold(data)) {
    126         psError(psErrorCodeLast(), false, "Unable to threshold images.");
    127         success = false;
    128         goto ESCAPE;
    129     }
    130 
    131     // Set up subtraction files
    132     if (!gDiffFilesIterateDown(config, GDIFF_FILES_SUB)) {
    133         psError(GDIFF_ERR_IO, false, "Unable to set up subtraction files.");
    134         success = false;
    135         goto ESCAPE;
    136     }
    137 
    138     if (!gDiffDefineOutput("GDIFF.OUTPUT", config)) {
    139         psError(psErrorCodeLast(), false, "Unable to define output.");
    140         success = false;
    141         goto ESCAPE;
    142     }
    143 
    144     if (!data->quality && !gDiffMakePSF(data)) {
    145         psError(psErrorCodeLast(), false, "Unable to generate PSF.");
    146         success = false;
    147         goto ESCAPE;
    148     }
    149 
    150     // Now we've got a PSF, blow away detections in case they're confused with real output detections
    151     {
    152         pmFPAview *view = gDiffViewReadout(); // View to readout
    153         pmReadout *out = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT");
    154         if (psMetadataLookup(out->analysis, "PSPHOT.DETECTIONS")) {
    155             psMetadataRemoveKey(out->analysis, "PSPHOT.DETECTIONS");
    156         }
    157         psFree(view);
    158     }
    159 
    160     if (!gDiffReadoutSubtract(config)) {
    161         psError(psErrorCodeLast(), false, "Unable to subtract images.");
    162         success = false;
    163         goto ESCAPE;
    164     }
    165     // dumpout(config, "diff.1.fits");
    166 
    167     // Close convolved files
    168     if (!gDiffFilesIterateUp(config, GDIFF_FILES_PSF | GDIFF_FILES_CONV)) {
    169         psError(GDIFF_ERR_IO, false, "Unable to close input files.");
    170         success = false;
    171         goto ESCAPE;
    172     }
    173     // dumpout(config, "diff.2a.fits");
    174     if (noConvolve) {
    175         if (!gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) {
    176             psError(GDIFF_ERR_IO, false, "Unable to close input files.");
    177             success = false;
    178             goto ESCAPE;
    179         }
    180     }
    181    
    182     // Higher order background subtraction using psphot
    183     if (!gDiffBackground(config)) {
    184         psError(psErrorCodeLast(), false, "Unable to subtract background.");
    185         success = false;
    186         goto ESCAPE;
    187     }
    188     // dumpout(config, "diff.2b.fits");
    189 
    190     // Perform Variance correction (rescale within a modest range)
    191     if (!gDiffVarianceRescale(config, data)) {
    192         psError(psErrorCodeLast(), false, "Unable to rescale variance.");
    193         success = false;
    194         goto ESCAPE;
    195     }
    196     // dumpout(config, "diff.2c.fits");
    197 
    198     if (data->quality) {
    199         // Done all we can do up to this point
    200         success = false;
    201         goto ESCAPE;
    202     }
    203 
    204     if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_SUB)) {
    205         psError(GDIFF_ERR_IO, false, "Unable to set up photometry files.");
    206         success = false;
    207         goto ESCAPE;
    208     }
    209 
    210     if (!data->quality && !gDiffReadoutPhotometry("GDIFF.OUTPUT", data)) {
    211         psError(psErrorCodeLast(), false, "Unable to perform photometry.");
    212         success = false;
    213         goto ESCAPE;
    214     }
    215     // dumpout(config, "diff.3.fits");
    216 
    217     // forced photometry for positive image 1
    218     if (data->forcedPhot1 && !data->quality) {
    219         if (!gDiffReadoutForcedPhot("GDIFF.FORCED1.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS1.SOURCES", data)) {
    220             psError(psErrorCodeLast(), false, "Unable to perform photometry.");
    221             success = false;
    222             goto ESCAPE;
    223         }
    224     }
    225 
    226     // forced photometry for positive image 2
    227     if (data->forcedPhot2 && !data->quality) {
    228         if (!gDiffReadoutForcedPhot("GDIFF.FORCED2.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS2.SOURCES", data)) {
    229             psError(psErrorCodeLast(), false, "Unable to perform photometry.");
    230             success = false;
    231             goto ESCAPE;
    232         }
    233     }
    234 
    235     if (!gDiffFilesIterateUp(config, GDIFF_FILES_PHOT_SUB)) {
    236         psError(GDIFF_ERR_IO, false, "Unable to set up photometry files.");
    237         success = false;
    238         goto ESCAPE;
    239     }
    240 
    241     // Perform statistics on the cell
    242     if (!gDiffReadoutStats(data)) {
    243         psError(psErrorCodeLast(), false, "Unable to collect statistics");
    244         success = false;
    245         goto ESCAPE;
    246     }
    247     // Do Mask Stats
    248     {
    249         pmFPAview *view = gDiffViewReadout(); // View to readout
    250         if (!gDiffMaskStats(config, view,data->stats)) {
    251             psError(psErrorCodeLast(), false, "Unable to generate mask statistics");
    252             success = false;
    253             psFree(view);
    254             goto ESCAPE;
    255         }
    256         psFree(view);
    257     }
    258     // dumpout(config, "diff.4.fits");
    259    
    260     // generate the binned image used to write the jpeg
    261     if (!gDiffReadoutJpeg(config)) {
    262         psError(psErrorCodeLast(), false, "Unable to update.");
    263         success = false;
    264         goto ESCAPE;
    265     }
    266 
    267     if (data->inverse) {
    268         // Set up inverse subtraction files
    269         if (!gDiffFilesIterateDown(config, GDIFF_FILES_INV)) {
    270             psError(GDIFF_ERR_IO, false, "Unable to set up inverse files.");
    271             success = false;
    272             goto ESCAPE;
    273         }
    274 
    275         if (data->inverse && !gDiffDefineOutput("GDIFF.INVERSE", config)) {
    276             psError(psErrorCodeLast(), false, "Unable to define inverse.");
    277             success = false;
    278             goto ESCAPE;
    279         }
    280 
    281         if (!gDiffReadoutInverse(config)) {
    282             psError(psErrorCodeLast(), false, "Unable to invert images.");
    283             success = false;
    284             goto ESCAPE;
    285         }
    286 
    287         //MEH -- need to also clear out detections or inv.cmf corrupted on update
    288         pmFPAview *viewinv = gDiffViewReadout(); // View to readout
    289         pmReadout *inv = pmFPAfileThisReadout(config->files, viewinv, "GDIFF.INVERSE");
    290         if (psMetadataLookup(inv->analysis, "PSPHOT.DETECTIONS")) {
    291             psMetadataRemoveKey(inv->analysis, "PSPHOT.DETECTIONS");
    292         }
    293         psFree(viewinv);
    294 
    295         // Close subtraction files and open inverse photometry files
    296         if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) {
    297             psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
    298             success = false;
    299             goto ESCAPE;
    300         }
    301 
    302         if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_INV)) {
    303             psError(GDIFF_ERR_IO, false, "Unable to set up inverse files.");
    304             success = false;
    305             goto ESCAPE;
    306         }
    307 
    308         if (!data->quality && !gDiffReadoutPhotometry("GDIFF.INVERSE", data)) {
    309             psError(psErrorCodeLast(), false, "Unable to perform photometry.");
    310             success = false;
    311             goto ESCAPE;
    312         }
    313 
    314         // Close inverse subtraction files
    315         if (!gDiffFilesIterateUp(config, GDIFF_FILES_INV | GDIFF_FILES_PHOT_INV)) {
    316             psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
    317             success = false;
    318             goto ESCAPE;
    319         }
    320     } else {
    321         // Close subtraction files
    322         // dumpout(config, "diff.5.fits");
    323         if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) {
    324             psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
    325             success = false;
    326             goto ESCAPE;
    327         }
    328     }
    329 
    330  ESCAPE:
    331     pmFPAfileDropInternal(config->files, "PSPHOT.BACKGND");
    332     pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
    333     pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
    334 
    335     return success;
     387    pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, name);
     388    fprintf (stderr, "HDU: %s, %llx\n", name, (long long unsigned int) file->fpa->hdu);
     389    return true;
    336390}
Note: See TracChangeset for help on using the changeset viewer.