code: def _fill_masked_intensity_resid(intensity_resid, bright_pixel_mask):
          
          filled_intensity = np.zeros((intensity_resid.data.shape))
          nebins = len(intensity_resid.data)
          for i in range(nebins):
              masked = bright_pixel_mask.data[i]
              unmasked = np.invert(masked)
              mean_intensity = intensity_resid.data[i][unmasked].mean()
              filled_intensity[i] = np.where(masked, mean_intensity, intensity_resid.data[i])
          return HpxMap(filled_intensity, intensity_resid.hpx)
docs: Fill the pixels used to compute the effective area correction with the mean intensity