Changeset 42885

Timestamp:

May 29, 2025, 3:28:25 PM (14 months ago)

Author:

tdeboer

Message:

build in safeties to prevent bad fits leading to random pattern row correction and bad cells

File:

• tags/ipp-ops-20220906/psModules/src/detrend/pmPattern.c (modified) (6 diffs)

tags/ipp-ops-20220906/psModules/src/detrend/pmPattern.c

@@ -239 +239 @@
 # define READNOISE 10 /* arbitrary number */
     float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    if(isnan(sigma)) sigma=20.;
     float delta = PS_MIN (thresh * sigma, 2*nominalAmplitude); 
     float lower = stats->robustMedian - delta; // Lower bound for data
@@ -246 +247 @@
     // if the noise from the background is too large 
     float significance = nominalAmplitude / sigma;
+    significance = abs(significance);
    
     // XXX EAM : arbitrary number
     if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      return true; 
+    }
+    //Add some safeties in case of bad fits
+    if (isnan(significance) || (background < 0.0)) {
       psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
       return true; 
@@ -554 +561 @@
 # define READNOISE 10 /* arbitrary number */
     float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    if(isnan(sigma)) sigma=20.; // XXX EAM : somewhat arbitrary number
     float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude); 
     float lower = stats->robustMedian - delta; // Lower bound for data
@@ -561 +569 @@
     // if the noise from the background is too large 
     float significance = nominalAmplitude / sigma;
-   
+    significance = abs(significance);
+      
     // XXX EAM : arbitrary number
     if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
@@ -569 +578 @@
       return true; 
     }
+    //Add some safeties in case of bad fits
+    if (isnan(significance) || (background < 0.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction due to bad fit for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      psFree(stats);
+      psFree(rng);
+      return true; 
+    }
+    
     psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
 # endif    
@@ -1005 +1022 @@
     psFree(meanMask);
 
+    return true;
+}
+
+// Compare the backs (cellBackgrounds) vector to each of the patterns loaded for this chip.
+// Choose the one that matches best & mask as appropriate
+// can we pass the backgrounds in using an image (8 x 8)
+bool pmPatternDeadCells(pmChip *chip, const psVector *backs, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+    PS_ASSERT_VECTOR_NON_NULL(backs, false);
+    PS_ASSERT_VECTOR_SIZE(backs, chip->cells->n, false);
+    PS_ASSERT_VECTOR_TYPE(backs, PS_TYPE_F32, false);
+
+    // Look for the dead cell pattern cube in the analysis metadata.
+    // NOTE: not all chips have the pattern, skip if not found.
+    bool mdok;
+    psImage *deadCellPattern = (psImage *) psMetadataLookupPtr (&mdok, chip->analysis, "PTN.DEAD.CELL");
+    if (!mdok) {
+        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "No DEAD CELL pattern for chip, skipping\n");
+        return true;
+    }
+
+    int numCells = backs->n;            // Number of cells
+    int numColsD = deadCellPattern->numCols;
+    int numRowsD = deadCellPattern->numRows;
+
+    assert (backs->n == numRowsD);
+
+    // The background values need to be normalized (divide by the median).
+    // First extract the valid data values
+    psVector *valid = psVectorAllocEmpty(backs->n, PS_TYPE_F32); // Mean for each cell
+    for (int i = 0; i < backs->n; i++) {
+        float value = backs->data.F32[i];
+        if (!isfinite(value)) { backs->data.F32[i] = NAN; continue; }
+        if (value > 50000.0)  { backs->data.F32[i] = NAN; continue; } // XXX warning: hard-wired value
+        psVectorAppend (valid, value);
+    }
+    if (valid->n == 0) {
+        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "No valid backgrounds, skipping\n");
+        psFree (valid);
+        return true;
+    }
+
+    // Second, calculate the median
+    psVectorSortInPlace (valid);
+    int midPt = valid->n / 2.0;
+    float median = valid->n % 2 ? valid->data.F32[midPt] : 0.5*(valid->data.F32[midPt] + valid->data.F32[midPt-1]);
+    psFree (valid);
+
+    // Finally, renormalize:
+    for (int i = 0; i < backs->n; i++) {
+        if (!isfinite(backs->data.F32[i])) { continue; }
+        backs->data.F32[i] /= median;
+    }
+
+    // there are 2 columns (value & mask) for each mode, plus the constant mode
+    int nModes = numColsD / 2 + 1;
+    
+    psVector *stdev = psVectorAlloc(nModes, PS_TYPE_F32); // Mean for each cell
+
+    // measure stdev for each comparison
+    for (int i = 0; i < nModes; i++) {
+        float Sum1 = 0.0;
+        float Sum2 = 0.0;
+        int   Npts = 0;
+
+        // choose the column with the background values for this mode
+        int nModeColumn = 2*(i - 1);
+
+        for (int j = 0; j < backs->n; j++) {
+
+            float valObs = backs->data.F32[j];
+            float valRef = (i == 0) ? 0.0 : deadCellPattern->data.F32[j][nModeColumn];
+
+            if (!isfinite(valObs)) continue;
+            if (!isfinite(valRef)) continue;
+
+            float dS = valObs - valRef;
+            Sum1 += dS;
+            Sum2 += dS*dS;
+            Npts ++;
+        }
+        if (Npts == 0) {
+            // is this is an error?
+            // no valid data, ignore this test
+            psLogMsg("psModules.detrend", PS_LOG_DETAIL, "No valid backgrounds, skipping\n");
+            psFree (stdev);
+            return true;
+        }
+
+        float mean = Sum1 / Npts;
+        float sigma = sqrt(Sum2 / Npts - mean*mean);
+        stdev->data.F32[i] = sigma;
+        fprintf (stderr, "mode: %d, stdev: %f\n", i, sigma);
+    }
+
+    // loop over stdev and choose the lowest one
+    int   minI = 0;
+    float minV = stdev->data.F32[minI];
+
+    for (int i = 1; i < nModes; i++) {
+        if (stdev->data.F32[i] < minV) {
+            minI = i;
+            minV = stdev->data.F32[i];
+        }
+    }
+    psFree (stdev);
+
+    if (minI == 0) return true;
+
+    int nModeColumnMask = 2*(minI - 1) + 1;
+
+    // now mask bad cells
+    for (int i = 0; i < numCells; i++) {
+        if (deadCellPattern->data.F32[i][nModeColumnMask] > 0) continue;
+
+        pmCell *cell = chip->cells->data[i]; // Cell of interest
+        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+
+        psImage *mask = ro->mask; // mask of interest
+        int numCols = mask->numCols, numRows = mask->numRows; // Size of image
+        
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal;
+            }
+        }
+    }
     return true;
 }