Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DM-42238: Propagates list of flat pair read noises to ptcDataset #227

Merged
merged 1 commit into from
Feb 2, 2024
Merged

DM-42238: Propagates list of flat pair read noises to ptcDataset #227

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
47 changes: 20 additions & 27 deletions python/lsst/cp/pipe/ptc/cpExtractPtcTask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -434,31 +434,6 @@ def run(self, inputExp, inputDims, taskMetadata, inputPhotodiodeData=None):
pdCalib.currentScale = self.config.photodiodeCurrentScale
monitorDiodeCharge[expId] = pdCalib.integrate()

# Get read noise. Try from the exposure, then try
# taskMetadata. This adds a get() for the exposures.
readNoiseLists = {}
for pairIndex, expRefs in inputExp.items():
# This yields an index (exposure_time, seq_num, or flux)
# and a pair of references at that index.
for expRef, expId in expRefs:
# This yields an exposure ref and an exposureId.
exposureMetadata = expRef.get(component="metadata")
metadataIndex = inputDims.index(expId)
thisTaskMetadata = taskMetadata[metadataIndex]

for ampName in ampNames:
if ampName not in readNoiseLists:
readNoiseLists[ampName] = [self.getReadNoise(exposureMetadata,
thisTaskMetadata, ampName)]
else:
readNoiseLists[ampName].append(self.getReadNoise(exposureMetadata,
thisTaskMetadata, ampName))

readNoiseDict = {ampName: 0.0 for ampName in ampNames}
for ampName in ampNames:
# Take median read noise value
readNoiseDict[ampName] = np.nanmedian(readNoiseLists[ampName])

# Output list with PTC datasets.
partialPtcDatasetList = []
# The number of output references needs to match that of input
Expand Down Expand Up @@ -535,6 +510,24 @@ def run(self, inputExp, inputDims, taskMetadata, inputPhotodiodeData=None):
im1Area, im2Area, imStatsCtrl, mu1, mu2 = self.getImageAreasMasksStats(exp1, exp2,
region=region)

# Get the read noise for each exposure
readNoise1 = dict()
readNoise2 = dict()
meanReadNoise = dict()

expMetadata1 = expRef1.get(component="metadata")
metadataIndex1 = inputDims.index(expId1)
thisTaskMetadata1 = taskMetadata[metadataIndex1]

expMetadata2 = expRef2.get(component="metadata")
metadataIndex2 = inputDims.index(expId2)
thisTaskMetadata2 = taskMetadata[metadataIndex2]

readNoise1[ampName] = self.getReadNoise(expMetadata1, thisTaskMetadata1, ampName)
readNoise2[ampName] = self.getReadNoise(expMetadata2, thisTaskMetadata2, ampName)

meanReadNoise[ampName] = np.nanmean([readNoise1[ampName], readNoise2[ampName]])
Alex-Broughton marked this conversation as resolved.
Show resolved Hide resolved

# We demand that both mu1 and mu2 be finite and greater than 0.
if not np.isfinite(mu1) or not np.isfinite(mu2) \
or ((np.nan_to_num(mu1) + np.nan_to_num(mu2)/2.) <= 0.0):
Expand Down Expand Up @@ -568,7 +561,7 @@ def run(self, inputExp, inputDims, taskMetadata, inputPhotodiodeData=None):
if self.config.doGain:
gain = self.getGainFromFlatPair(im1Area, im2Area, imStatsCtrl, mu1, mu2,
correctionType=self.config.gainCorrectionType,
readNoise=readNoiseDict[ampName])
readNoise=meanReadNoise[ampName])
else:
gain = np.nan

Expand Down Expand Up @@ -667,7 +660,7 @@ def run(self, inputExp, inputDims, taskMetadata, inputPhotodiodeData=None):
covariance=covArray[0, :, :],
covSqrtWeights=covSqrtWeights[0, :, :],
gain=gain,
noise=readNoiseDict[ampName],
noise=meanReadNoise[ampName],
histVar=histVar,
histChi2Dof=histChi2Dof,
kspValue=kspValue,
Expand Down
14 changes: 14 additions & 0 deletions python/lsst/cp/pipe/ptc/cpSolvePtcTask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -334,6 +334,7 @@ def run(self, inputCovariances, camera=None, detId=0):
ptcFitType=self.config.ptcFitType,
covMatrixSide=self.config.maximumRangeCovariancesAstier,
covMatrixSideFullCovFit=self.config.maximumRangeCovariancesAstierFullCovFit)

for partialPtcDataset in inputCovariances:
# Ignore dummy datasets
if partialPtcDataset.ptcFitType == 'DUMMY':
Expand All @@ -360,6 +361,8 @@ def run(self, inputCovariances, camera=None, detId=0):
partialPtcDataset.histChi2Dofs[ampName][0])
datasetPtc.kspValues[ampName] = np.append(datasetPtc.kspValues[ampName],
partialPtcDataset.kspValues[ampName][0])
datasetPtc.noiseList[ampName] = np.append(datasetPtc.noiseList[ampName],
partialPtcDataset.noise[ampName])
datasetPtc.covariances[ampName] = np.append(
datasetPtc.covariances[ampName].ravel(),
partialPtcDataset.covariances[ampName].ravel()
Expand Down Expand Up @@ -466,6 +469,17 @@ def run(self, inputCovariances, camera=None, detId=0):
# PhotonTransferCurveDataset object.
datasetPtc = self.fitMeasurementsToModel(datasetPtc)

# Initial validation of PTC fit.
for ampName in ampNames:
noise = np.nanmedian(datasetPtc.noiseList[ampName])
noiseFitted = np.sqrt(datasetPtc.noise[ampName])

# Check if noise is close to noiseFitted
if not np.isclose(noiseFitted, noise, rtol=0.05, atol=0.0):
Alex-Broughton marked this conversation as resolved.
Show resolved Hide resolved
self.log.warning(f"Read noise from PTC fit ({noiseFitted}) is not consistent "
f"with read noise measured from overscan ({noise}) for "
f"amplifier {ampName}. Try adjusting the fit range.")
Alex-Broughton marked this conversation as resolved.
Show resolved Hide resolved

if camera:
detector = camera[detId]
else:
Expand Down
Loading