Add unit test for cosmic ray detection

Author: isullivan

tests/test_detectAndMeasure.py

@@ -414,6 +414,46 @@ def _detection_wrapper(positive=True):
         _detection_wrapper(positive=True)
         _detection_wrapper(positive=False)
 
+    def test_mask_cosmic_rays(self):
+        """Run detection on a difference image containing a cosmic ray.
+        """
+        # Set up the simulated images
+        noiseLevel = 1.
+        staticSeed = 1
+        fluxLevel = 500
+        xSize = 400
+        ySize = 400
+        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "xSize": xSize, "ySize": ySize}
+        science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
+        matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
+        crMask = science.mask.getPlaneBitMask("CR")
+
+        # Configure the detection Task
+        detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
+
+        # Test that no CRs are detected
+        transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
+        science.maskedImage += transients.maskedImage
+        difference = science.clone()
+        difference.maskedImage -= matchedTemplate.maskedImage
+        output = detectionTask.run(science, matchedTemplate, difference, sources)
+        crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
+        self.assertTrue(np.all(crMaskSet == 0))
+
+        crX0 = round(sources.getX()[0] - science.getX0())
+        crY0 = round(sources.getY()[0] - science.getY0())
+        crX1 = round(sources.getX()[1] - science.getX0())
+        crY1 = round(sources.getY()[1] - science.getY0())
+        # Add CR-like shape and check that CR is detected
+        # Pick two locations on top of sources, since that is what is likely to
+        # be missed in the first stage of CR rejection.
+        difference.image.array[crX0:crX0+1, crY0:crY0+5] += 1234  # Arbitrary but large flux for the CRs
+        difference.image.array[crX1:crX1+5, crY1:crY1+1] += 1234
+        output = detectionTask.run(science, matchedTemplate, difference)
+        crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
+        self.assertTrue(np.all(crMaskSet[crX0:crX0+1, crY0:crY0+5]))
+        self.assertTrue(np.all(crMaskSet[crX1:crX1+5, crY1:crY1+1]))
+
     def test_missing_mask_planes(self):
         """Check that detection runs with missing mask planes.
         """