Changeset 28794 for branches/eam_branches/ipp-20100621/ippToPsps/src/ippToPspsBatchDetection.c

Timestamp:

Jul 30, 2010, 9:31:50 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

File:

• branches/eam_branches/ipp-20100621/ippToPsps/src/ippToPspsBatchDetection.c (modified) (14 diffs)

branches/eam_branches/ipp-20100621/ippToPsps/src/ippToPspsBatchDetection.c

@@ -14 +14 @@
 
 // Gets uncalibrated instrumental flux from magnitude
-static __inline bool ippToPsps_getFlux(const float zeroPoint, const float exposureTime, const float magnitude, float* flux) {
+static __inline bool ippToPsps_getFlux(const float exposureTime, const float magnitude, float* flux, const float magnitudeErr, float* fluxErr) {
 
     *flux = powf(10.0, -0.4*magnitude) / exposureTime;
-    return (!isfinite(*flux) || *flux < 0.000001) ? false : true;
+    if (!isfinite(*flux) || *flux < 0.000001) return false;
+    if (fluxErr) *fluxErr = fabsf((magnitudeErr * *flux)/1.085736);
+    //  if (fluxErr)    printf("Mag = %03.03f, Flux = %03.03f, Mag err = %03.03f, Flux Err = %03.03f\n", magnitude, *flux, magnitudeErr, *fluxErr);
+         
+    return true;
 }
 
@@ -46 +50 @@
     fits_get_hdrspace(fitsIn, &nKeys, NULL, &status);
 
-    float zptObs, zeroPoint, exposureTime;
+    float zptObs, exposureTime;
     char filterType[20];
     double obsTime;
@@ -52 +56 @@
     double expTime;
     status=0; fits_read_key(fitsIn, TFLOAT, "ZPT_OBS", &zptObs, NULL, &status);
-    status=0; fits_read_key(fitsIn, TFLOAT, "ZPT_REF", &zeroPoint, NULL, &status);
     status=0; fits_read_key(fitsIn, TFLOAT, "EXPREQ", &exposureTime, NULL, &status);
     status=0; fits_read_key(fitsIn, TSTRING, "FILTERID", filterType, NULL, &status);
     status=0; fits_read_key(fitsIn, TDOUBLE, "MJD-OBS", &expStart, NULL, &status);
     status=0; fits_read_key(fitsIn, TDOUBLE, "EXPTIME", &expTime, NULL, &status);
-    obsTime = expStart + (expTime/172800); // exp start plus half exp time (converted from secs to days)
+    obsTime = expStart + (expTime/172800.0); // exp start plus half exp time (converted from secs to days)
    
     ippToPspsConfig_writeTable(this->config, fitsIn, this->fitsOut, 1, "FrameMeta", true);
@@ -83 +86 @@
 
     // stuff to keep from psf.hdr header
-    int sourceId = -1;
-    int imageId = -1;
-    float fwhmMaj;
-    float fwhmMin;
-    float momentMaj;
-    float momentMin;
-    float psfFwhm;
-    float momentFwhm;
+    int sourceId = -1, imageId = -1;
+    float fwhmMaj, fwhmMin, momentMaj, momentMin, psfFwhm, momentFwhm;
     long pspsImageId = -1;
 
@@ -96 +93 @@
     SkyList *skylist = NULL;
     dvoDetection *dvoDetections = NULL;
-    int maxDvoDetId = -1;
-    int numDvoDetections = -1;
+    int maxDvoDetId = -1, numDvoDetections = -1;
     Image *image = NULL;
 
@@ -108 +104 @@
     strftime(timeStr, sizeof(timeStr), "%Y-%m-%d", ts);
 
-    uint32_t totalDetections = 0;
-    uint32_t s,d;
-    uint32_t invalidDvoRows;
-    uint32_t smfJumps;
-    uint32_t unmatched;
+    uint32_t s,d, invalidDvoRows, smfJumps, unmatched, totalDetections = 0;
 
     long longnull = -999;
@@ -118 +110 @@
     int anynull = 0;
 
-    char ccdNumber[3];
-    char extensionName[15];
-
+    char ccdNumber[3], extensionName[15];
+
+    // for storing FITS column data
     long ippIDet[MAXDETECT];
-    float instMag[MAXDETECT];
-    float instMagErr[MAXDETECT];
-    float peakMag[MAXDETECT];
+    float instMag[MAXDETECT], instMagErr[MAXDETECT], peakMag[MAXDETECT];
     uint64_t flags[MAXDETECT];
-    long objID[MAXDETECT];
-    long detectID[MAXDETECT];
-    long ippObjID[MAXDETECT];
-    long ippDetectID[MAXDETECT];
-    long imageID[MAXDETECT];
+    long objID[MAXDETECT], detectID[MAXDETECT], ippObjID[MAXDETECT], ippDetectID[MAXDETECT], imageID[MAXDETECT];
     double obsTimes[MAXDETECT];
-    float instFlux[MAXDETECT];
-    float instFluxErr[MAXDETECT];
-    float peakFlux[MAXDETECT];
-    int8_t filterIDs[MAXDETECT];
-    int8_t surveyIDs[MAXDETECT];
+    float instFlux[MAXDETECT], instFluxErr[MAXDETECT], peakFlux[MAXDETECT];
+    int8_t filterIDs[MAXDETECT], surveyIDs[MAXDETECT];
 
     char** assocDate = (char**)calloc(MAXDETECT, sizeof(char**));
@@ -149 +132 @@
     }
 
-    long maxObjID = LONG_MIN; 
-    long minObjID = LONG_MAX;
+    long maxObjID = LONG_MIN, minObjID = LONG_MAX;
     uint64_t imageFlags;
     short nOta = 0;
     long i;
     bool isDuplicate;
-    uint32_t numOfDuplicates;
-    uint32_t numInvalidFlux;
-    long numDetectionsOut;
-    long totalDetectionsOut = 0;
+    uint32_t numOfDuplicates, numInvalidFlux;
+    long numDetectionsOut, totalDetectionsOut = 0, numPhotoRef, totalNumPhotoRef = 0;
     long removeList[MAXDETECT];
     long thisObjId;
+    bool firstTimeIn = true, peakFluxOk, instFluxOk;
+    int ippIDetNum, instMagNum, instMagErrNum, peakMagNum;
 
     // loop round all 60 chips
@@ -197 +179 @@
             status=0; fits_read_key(fitsIn, TLONG, "IMAGEID", &imageId, NULL, &status);
             status=0; fits_read_key(fitsIn, TLONG, "SOURCEID", &sourceId, NULL, &status);
+            status=0; fits_read_key(fitsIn, TLONG, "NASTRO", &numPhotoRef, NULL, &status);
+            totalNumPhotoRef += numPhotoRef; // total up for storing in FrameMeta
 
             // access DVO database
@@ -202 +186 @@
             if (skylist == NULL) {
                 psError(PS_ERR_IO, false, 
-                        "DVO: can't find SkyList for sourceId='%d' imageId='%d' (CCD = XY%s): skipping\n", 
+                        "DVO: can't find SkyList for sourceId='%d' imageId='%d' (CCD = XY%s): skipping this chip\n",
                         sourceId, imageId, ccdNumber);
                 continue;
@@ -251 +235 @@
             s = d = totalDetections = invalidDvoRows = smfJumps = unmatched = 0;
 
+            // determine column numbers for certain IPP detection columns
+            if (firstTimeIn) {
+
+                status=0;fits_get_colnum(fitsIn, CASESEN, "IPP_IDET", &ippIDetNum, &status);
+                if (status) psError(PS_ERR_IO, false, "Unable to read col num for IPP_IDET");
+                status=0;fits_get_colnum(fitsIn, CASESEN, "PSF_INST_MAG", &instMagNum, &status);
+                if (status) psError(PS_ERR_IO, false, "Unable to read col num for PSF_INST_MAG");
+                status=0;fits_get_colnum(fitsIn, CASESEN, "PSF_INST_MAG_SIG", &instMagErrNum, &status);
+                if (status) psError(PS_ERR_IO, false, "Unable to read col num for PSF_INST_MAG_SIG");
+                status=0;fits_get_colnum(fitsIn, CASESEN, "PEAK_FLUX_AS_MAG", &peakMagNum, &status);
+                if (status) psError(PS_ERR_IO, false, "Unable to read col num for PEAK_FLUX_AS_MAG");
+
+                firstTimeIn=false;
+            }
+
             anynull = 0;
-            fits_read_col(fitsIn, TLONG, 1, 1, 1, nDet, &longnull, ippIDet, &anynull, &status);
-            fits_read_col(fitsIn, TFLOAT, 8, 1, 1, nDet, &floatnull, instMag, &anynull, &status);
-            fits_read_col(fitsIn, TFLOAT, 9, 1, 1, nDet, &floatnull, instMagErr, &anynull, &status);
-            fits_read_col(fitsIn, TFLOAT, 12, 1, 1, nDet, &floatnull, peakMag, &anynull, &status);
+            fits_read_col(fitsIn, TLONG, ippIDetNum, 1, 1, nDet, &longnull, ippIDet, &anynull, &status);
+            fits_read_col(fitsIn, TFLOAT, instMagNum, 1, 1, nDet, &floatnull, instMag, &anynull, &status);
+            fits_read_col(fitsIn, TFLOAT, instMagErrNum, 1, 1, nDet, &floatnull, instMagErr, &anynull, &status);
+            fits_read_col(fitsIn, TFLOAT, peakMagNum, 1, 1, nDet, &floatnull, peakMag, &anynull, &status);
 
             // DVO detections are ordered by IPP_IDET, which increments from 0 in SMF table
@@ -297 +296 @@
                     obsTimes[s] = obsTime;
 
-                    ippToPsps_getFlux(zeroPoint, exposureTime, instMagErr[s], &instFluxErr[s]);
+                    // calculate flux and flux errors from magnitudes
+                    peakFluxOk = ippToPsps_getFlux(exposureTime, peakMag[s], &peakFlux[s], 0.0, NULL);
+                    instFluxOk = ippToPsps_getFlux(exposureTime, instMag[s], &instFlux[s], instMagErr[s], &instFluxErr[s]);
 
                     // check for invalid flux values (if not already labelled as a duplicate)
-                    if ((!ippToPsps_getFlux(zeroPoint, exposureTime, peakMag[s], &peakFlux[s]) ||
-                                !ippToPsps_getFlux(zeroPoint, exposureTime, instMag[s], &instFlux[s])) && !isDuplicate) {
+                    if ( (!peakFluxOk || !instFluxOk) && !isDuplicate) {
                         removeList[numOfDuplicates+numInvalidFlux] = s+1;
                         numInvalidFlux++;
@@ -379 +379 @@
     if (fits_movnam_hdu(this->fitsOut, BINARY_TBL, "FrameMeta", 0, &status))
         psError(PS_ERR_IO, false, "Can't move back to FrameMeta extension to write number of OTAs\n");
-    else 
+    else { 
         fits_write_col(this->fitsOut, TSHORT, FRAMEMETA_NOTA, 1, 1, 1, &nOta, &status);
-
+        fits_write_col(this->fitsOut, TSHORT, FRAMEMETA_NUMPHOTOREF, 1, 1, 1, &totalNumPhotoRef, &status);
+    }
     status=0;
     if (fits_close_file(fitsIn, &status)) fits_report_error(stderr, status);
 
     // write results
-    if (this->resultsFile) {
-
-        if (fprintf(this->resultsFile, "%ld\n", minObjID) < 1 ) 
-            psError(PS_ERR_IO, false, "Unable to write min Object ID to'%s'\n", this->resultsPath);
-        if (fprintf(this->resultsFile, "%ld\n", maxObjID) < 1 ) 
-            psError(PS_ERR_IO, false, "Unable to write max Object ID to'%s'\n", this->resultsPath);
-        if (fprintf(this->resultsFile, "%ld\n", totalDetectionsOut) < 1 ) 
-            psError(PS_ERR_IO, false, "Unable to write totalDetectionsOut to'%s'\n", this->resultsPath);
+    if (this->resultsXmlDoc) {
+
+        xmlNodePtr rootNode = xmlDocGetRootElement(this->resultsXmlDoc);
+        char tmp[100];
+        sprintf(tmp, "%ld", minObjID);
+        xmlNewChild(rootNode, NULL, BAD_CAST "minObjID", BAD_CAST tmp);
+        sprintf(tmp, "%ld", maxObjID);
+        xmlNewChild(rootNode, NULL, BAD_CAST "maxObjID", BAD_CAST tmp);
+        sprintf(tmp, "%ld", totalDetectionsOut);
+        xmlNewChild(rootNode, NULL, BAD_CAST "totalDetections", BAD_CAST tmp);
     }
 
@@ -403 +406 @@
 }
 
-