Index: /branches/pap/psModules/src/objects/Makefile.am =================================================================== --- /branches/pap/psModules/src/objects/Makefile.am (revision 25326) +++ /branches/pap/psModules/src/objects/Makefile.am (revision 25327) @@ -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: /branches/pap/psModules/src/objects/pmDetEff.c =================================================================== --- /branches/pap/psModules/src/objects/pmDetEff.c (revision 25327) +++ /branches/pap/psModules/src/objects/pmDetEff.c (revision 25327) @@ -0,0 +1,139 @@ +#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]); + psMetadataAddF32(row, PS_LIST_TAIL, "COUNTS", 0, "Number of sources recovered", + de->magOffsets->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "DIFF.MEAN", 0, "Mean magnitude difference", + de->magOffsets->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "DIFF.STDEV", 0, "Stdev magnitude difference", + de->magOffsets->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "ERR.MEAN", 0, "Mean magnitude error", + de->magOffsets->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; +} + + +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; +} + Index: /branches/pap/psModules/src/objects/pmDetEff.h =================================================================== --- /branches/pap/psModules/src/objects/pmDetEff.h (revision 25327) +++ /branches/pap/psModules/src/objects/pmDetEff.h (revision 25327) @@ -0,0 +1,64 @@ +#ifndef PM_DET_EFF_H +#define PM_DET_EFF_H + +#include +#include + +#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 + ); + +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 + ); + +pmDetEff *pmDetEffRead(psFits *fits, // FITS file from which to read + 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_F32, 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