Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 25256) +++ trunk/psModules/src/objects/Makefile.am (revision 25383) @@ -53,5 +53,6 @@ pmGrowthCurveGenerate.c \ pmGrowthCurve.c \ - pmSourceMatch.c + pmSourceMatch.c \ + pmDetEff.c EXTRA_DIST = \ @@ -90,5 +91,6 @@ pmTrend2D.h \ pmGrowthCurve.h \ - pmSourceMatch.h + pmSourceMatch.h \ + pmDetEff.h CLEANFILES = *~ Index: trunk/psModules/src/objects/pmDetEff.c =================================================================== --- trunk/psModules/src/objects/pmDetEff.c (revision 25383) +++ trunk/psModules/src/objects/pmDetEff.c (revision 25383) @@ -0,0 +1,168 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#include "pmDetEff.h" + + +static void detEffFree(pmDetEff *de) +{ + psFree(de->magOffsets); + psFree(de->counts); + psFree(de->magDiffMean); + psFree(de->magDiffStdev); + psFree(de->magErrMean); +} + + +pmDetEff *pmDetEffAlloc(float magRef, int numSources, int numBins) +{ + pmDetEff *de = psAlloc(sizeof(pmDetEff)); // Detection efficiency, to return + psMemSetDeallocator(de, (psFreeFunc)detEffFree); + + de->magRef = magRef; + de->numSources = numSources; + de->numBins = numBins; + + de->magOffsets = NULL; + de->counts = NULL; + de->magDiffMean = NULL; + de->magDiffStdev = NULL; + de->magErrMean = NULL; + + return de; +} + + +bool pmDetEffWrite(psFits *fits, pmDetEff *de, const psMetadata *header, const char *extname) +{ + PS_ASSERT_FITS_NON_NULL(fits, false); + PS_ASSERT_FITS_WRITABLE(fits, false); + PM_ASSERT_DETEFF_RESULTS(de, false); + + psArray *table = psArrayAlloc(de->numBins); // Table to write + for (int i = 0; i < de->numBins; i++) { + psMetadata *row = psMetadataAlloc(); // Table row + psMetadataAddF32(row, PS_LIST_TAIL, "OFFSET", 0, "Magnitude offset from reference", + de->magOffsets->data.F32[i]); + psMetadataAddS32(row, PS_LIST_TAIL, "COUNTS", 0, "Number of sources recovered", + de->counts->data.S32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "DIFF.MEAN", 0, "Mean magnitude difference", + de->magDiffMean->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "DIFF.STDEV", 0, "Stdev magnitude difference", + de->magDiffStdev->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "ERR.MEAN", 0, "Mean magnitude error", + de->magErrMean->data.F32[i]); + table->data[i] = row; + } + + psMetadata *deHeader = psMetadataCopy(NULL, header); // Header for detection efficiency + psMetadataAddF32(deHeader, PS_LIST_TAIL, "DETEFF.MAGREF", PS_META_REPLACE, "Magnitude reference", + de->magRef); + psMetadataAddS32(deHeader, PS_LIST_TAIL, "DETEFF.NUM", PS_META_REPLACE, "Number of fake sources injected", + de->numSources); + + if (!psFitsWriteTable(fits, deHeader, table, extname)) { + psError(PS_ERR_IO, false, "Unable to write detection efficiency table."); + psFree(table); + psFree(deHeader); + return false; + } + + psFree(table); + psFree(deHeader); + + return true; +} + +bool pmReadoutWriteDetEff(psFits *fits, const pmReadout *readout, + const psMetadata *header, const char *extname) +{ + PM_ASSERT_READOUT_NON_NULL(readout, false); + + bool mdok; // Status of MD lookup + pmDetEff *de = psMetadataLookupPtr(&mdok, readout->analysis, PM_DETEFF_ANALYSIS); // Detection efficiency + if (!mdok || !de) { + // Wrote everything there was to write + return true; + } + return pmDetEffWrite(fits, de, header, extname); +} + + +pmDetEff *pmDetEffRead(psFits *fits, const char *extname) +{ + PS_ASSERT_FITS_NON_NULL(fits, false); + PS_ASSERT_STRING_NON_EMPTY(extname, false); + + if (!psFitsMoveExtNameClean(fits, extname)) { + // Nothing to read + return NULL; + } + + psMetadata *header = psFitsReadHeader(NULL, fits); // Header for table + if (!header) { + psError(PS_ERR_IO, false, "Unable to read FITS header"); + return NULL; + } + + int numBins = psFitsTableSize(fits);// Size of table + bool mdok; // Status of MD lookup + int numSources = psMetadataLookupS32(&mdok, header, "DETEFF.NUM"); // Number of fake sources injected + if (!mdok) { + psError(PS_ERR_UNKNOWN, false, "Unable to find number of sources"); + psFree(header); + return NULL; + } + float magRef = psMetadataLookupF32(&mdok, header, "DETEFF.MAGREF"); // Magnitude reference + if (!mdok) { + psError(PS_ERR_UNKNOWN, false, "Unable to find magnitude reference"); + psFree(header); + return NULL; + } + psFree(header); + + pmDetEff *de = pmDetEffAlloc(magRef, numSources, numBins); // Detection efficiency + de->magOffsets = psVectorAlloc(numBins, PS_TYPE_F32); + de->counts = psVectorAlloc(numBins, PS_TYPE_S32); + de->magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); + de->magDiffStdev = psVectorAlloc(numBins, PS_TYPE_F32); + de->magErrMean = psVectorAlloc(numBins, PS_TYPE_F32); + + psArray *table = psFitsReadTable(fits); // FITS table + if (!table) { + psError(PS_ERR_IO, false, "Unable to read detection efficiency table."); + psFree(de); + return false; + } + + for (int i = 0; i < numBins; i++) { + psMetadata *row = table->data[i]; // Table row + de->magOffsets->data.F32[i] = psMetadataLookupF32(NULL, row, "OFFSET"); + de->counts->data.S32[i] = psMetadataLookupS32(NULL, row, "COUNTS"); + de->magDiffMean->data.F32[i] = psMetadataLookupF32(NULL, row, "DIFF.MEAN"); + de->magDiffStdev->data.F32[i] = psMetadataLookupF32(NULL, row, "DIFF.STDEV"); + de->magErrMean->data.F32[i] = psMetadataLookupF32(NULL, row, "ERR.MEAN"); + } + + psFree(table); + return de; +} + +bool pmReadoutReadDetEff(psFits *fits, const pmReadout *readout, const char *extname) +{ + PM_ASSERT_READOUT_NON_NULL(readout, false); + + pmDetEff *de = pmDetEffRead(fits, extname); + if (!de) { + if (psErrorCodeLast() != PS_ERR_NONE) { + return false; + } + return true; + } + + return psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, PM_DETEFF_ANALYSIS, + PS_META_REPLACE | PS_DATA_UNKNOWN, "Detection efficiency", de); +} Index: trunk/psModules/src/objects/pmDetEff.h =================================================================== --- trunk/psModules/src/objects/pmDetEff.h (revision 25383) +++ trunk/psModules/src/objects/pmDetEff.h (revision 25383) @@ -0,0 +1,81 @@ +#ifndef PM_DET_EFF_H +#define PM_DET_EFF_H + +#include +#include + +#include "pmFPA.h" + +#define PM_DETEFF_ANALYSIS "DETEFF" // Location of detection efficiency on pmReadout.analysis + +// Detection efficiency characterisation +typedef struct { + float magRef; // Reference magnitude + int numSources; // Number of sources + int numBins; // Number of bins + psVector *magOffsets; // Magnitude offsets for each bin + psVector *counts; // Counts of sources recovered for each bin + psVector *magDiffMean; // Mean magnitude difference for each bin + psVector *magDiffStdev; // Stdev of magnitude difference for each bin + psVector *magErrMean; // Mean magnitude error for each bin +} pmDetEff; + + +/// Allocator +pmDetEff *pmDetEffAlloc(float magRef, // Reference magnitude + int numSources, // Number of sources + int numBins // Number of bins + ); + +/// Write detection efficiency to FITS file +bool pmDetEffWrite(psFits *fits, // FITS file to which to write + pmDetEff *deteff, // Detection efficiency to write + const psMetadata *header, // Header to write + const char *extname // Extension name + ); + +/// Write detection efficiency from a readout to a FITS file +bool pmReadoutWriteDetEff(psFits *fits,// FITS file to which to write + const pmReadout *readout, // Readout with detection efficiency + const psMetadata *header, // Header to write + const char *extname // Extension name + ); + +/// Read detection efficiency +pmDetEff *pmDetEffRead(psFits *fits, // FITS file from which to read + const char *extname // Extension name + ); + +/// Read detection efficiency into a readout +bool pmReadoutReadDetEff(psFits *fits,// FITS file to which to write + const pmReadout *readout, // Readout with detection efficiency + const char *extname // Extension name + ); + +#define PM_ASSERT_DETEFF_NON_NULL(DE, RETURN) { \ + if (!(DE)) { \ + psError(PS_ERR_UNEXPECTED_NULL, true, "Detection efficiency %s is NULL", #DE); \ + return RETURN; \ + } \ +} + +#define PM_ASSERT_DETEFF_RESULTS(DE, RETURN) { \ + PM_ASSERT_DETEFF_NON_NULL(DE, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->magOffsets, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->magOffsets, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->magOffsets, PS_TYPE_F32, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->counts, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->counts, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->counts, PS_TYPE_S32, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->magDiffMean, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->magDiffMean, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->magDiffMean, PS_TYPE_F32, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->magDiffStdev, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->magDiffStdev, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->magDiffStdev, PS_TYPE_F32, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->magErrMean, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->magErrMean, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->magErrMean, PS_TYPE_F32, RETURN); \ +} + +#endif Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 25256) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 25383) @@ -42,8 +42,82 @@ #include "pmSource.h" #include "pmModelClass.h" +#include "pmDetEff.h" #include "pmSourceIO.h" #define BLANK_HEADERS "BLANK.HEADERS" // Name of metadata in camera configuration containing header names // for putting values into a blank PHU + +// lookup the EXTNAME values used for table data and image header segments +static bool sourceExtensions(psString *headname, // Extension name for header + psString *dataname, // Extension name for data + psString *deteffname, // Extension name for detection efficiency + psString *xsrcname, // Extension name for extended sources + psString *xfitname, // Extension name for extended fits + const pmFPAfile *file, // File of interest + const pmFPAview *view // View to level of interest + ) +{ + bool status; // Status of MD lookup + + // Menu of EXTNAME rules + psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); + if (!menu) { + psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); + return false; + } + + // EXTNAME for image header + if (headname) { + const char *rule = psMetadataLookupStr(&status, menu, "CMF.HEAD"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config"); + return false; + } + *headname = pmFPAfileNameFromRule(rule, file, view); + } + + // EXTNAME for table data + if (dataname) { + const char *rule = psMetadataLookupStr(&status, menu, "CMF.DATA"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config"); + return false; + } + *dataname = pmFPAfileNameFromRule(rule, file, view); + } + + // EXTNAME for detection efficiency + if (deteffname) { + const char *rule = psMetadataLookupStr(&status, menu, "CMF.DETEFF"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DETEFF in EXTNAME.RULES in camera.config"); + return false; + } + *deteffname = pmFPAfileNameFromRule(rule, file, view); + } + + // EXTNAME for extended source data table + if (xsrcname) { + const char *rule = psMetadataLookupStr(&status, menu, "CMF.XSRC"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config"); + return false; + } + *xsrcname = pmFPAfileNameFromRule (rule, file, view); + } + + if (xfitname) { + // EXTNAME for extended source data table + const char *rule = psMetadataLookupStr(&status, menu, "CMF.XFIT"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config"); + return false; + } + *xfitname = pmFPAfileNameFromRule (rule, file, view); + } + + return true; +} + // translations between psphot object types and dophot object types @@ -271,8 +345,4 @@ char *exttype = NULL; - char *dataname = NULL; - char *xsrcname = NULL; - char *xfitname = NULL; - char *headname = NULL; // if sources is NULL, write out an empty table @@ -354,49 +424,12 @@ // define the EXTNAME values for the different data segments: - { - // lookup the EXTNAME values used for table data and image header segments - char *rule = NULL; - - // Menu of EXTNAME rules - psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); - if (!menu) { - psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); - return false; - } - - // EXTNAME for image header - rule = psMetadataLookupStr(&status, menu, "CMF.HEAD"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config"); - return false; - } - headname = pmFPAfileNameFromRule (rule, file, view); - - // EXTNAME for table data - rule = psMetadataLookupStr(&status, menu, "CMF.DATA"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config"); - return false; - } - dataname = pmFPAfileNameFromRule (rule, file, view); - - if (XSRC_OUTPUT) { - // EXTNAME for extended source data table - rule = psMetadataLookupStr(&status, menu, "CMF.XSRC"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config"); - return false; - } - xsrcname = pmFPAfileNameFromRule (rule, file, view); - } - if (XFIT_OUTPUT) { - // EXTNAME for extended source data table - rule = psMetadataLookupStr(&status, menu, "CMF.XFIT"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config"); - return false; - } - xfitname = pmFPAfileNameFromRule (rule, file, view); - } + psString headname = NULL; + psString dataname = NULL; + psString deteffname = NULL; + psString xsrcname = NULL; + psString xfitname = NULL; + if (!sourceExtensions(&headname, &dataname, &deteffname, XSRC_OUTPUT ? &xsrcname : NULL, + XFIT_OUTPUT ? &xfitname : NULL, file, view)) { + return false; } @@ -480,49 +513,54 @@ // XXX these are case-sensitive since the EXTYPE is case-sensitive - status = false; + status = true; if (!strcmp (exttype, "SMPDATA")) { - status = pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname); + status &= pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname); } if (!strcmp (exttype, "PS1_DEV_0")) { - status = pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname); + status &= pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname); } if (!strcmp (exttype, "PS1_DEV_1")) { - status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname); + status &= pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname); } if (!strcmp (exttype, "PS1_CAL_0")) { - status = pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname); + status &= pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname); } if (!strcmp (exttype, "PS1_V1")) { - status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname); + status &= pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname); } if (!strcmp (exttype, "PS1_V2")) { - status = pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname); - } + status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname); + } + + if (deteffname) { + status &= pmReadoutWriteDetEff(file->fits, readout, outhead, deteffname); + } + if (xsrcname) { if (!strcmp (exttype, "PS1_DEV_1")) { - status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe); } if (!strcmp (exttype, "PS1_CAL_0")) { - status = pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); + status &= pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } if (!strcmp (exttype, "PS1_V1")) { - status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe); } if (!strcmp (exttype, "PS1_V2")) { - status = pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe); } } if (xfitname) { if (!strcmp (exttype, "PS1_DEV_1")) { - status = pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname); } if (!strcmp (exttype, "PS1_CAL_0")) { - status = pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname); + status &= pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname); } if (!strcmp (exttype, "PS1_V1")) { - status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname); } if (!strcmp (exttype, "PS1_V2")) { - status = pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname); } } @@ -572,6 +610,6 @@ // not needed if only one chip if (file->fpa->chips->n == 1) { - pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config); - return true; + pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config); + return true; } @@ -885,26 +923,11 @@ hdu = pmFPAviewThisHDU (view, file->fpa); - // lookup the EXTNAME values used for table data and image header segments - char *rule = NULL; - // Menu of EXTNAME rules - psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); - if (!menu) { - psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); - return false; - } - // EXTNAME for image header - rule = psMetadataLookupStr(&status, menu, "CMF.HEAD"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config"); - return false; - } - char *headname = pmFPAfileNameFromRule (rule, file, view); - // EXTNAME for table data - rule = psMetadataLookupStr(&status, menu, "CMF.DATA"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config"); - return false; - } - char *dataname = pmFPAfileNameFromRule (rule, file, view); + // define the EXTNAME values for the different data segments: + psString headname = NULL; + psString dataname = NULL; + psString deteffname = NULL; + if (!sourceExtensions(&headname, &dataname, &deteffname, NULL, NULL, file, view)) { + return false; + } // advance to the IMAGE HEADER extension @@ -958,4 +981,9 @@ sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header); } + + if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) { + psError(PS_ERR_IO, false, "Unable to read detection efficiency"); + return false; + } } @@ -1070,3 +1098,3 @@ } - +