NOTE!! This is now out of date - I haven't got around to writing a new one.
subroutine lccorr
foreach (srcTS, bkgTS) {
# Do basic existence and sense checks on the datasets.
RATE = COUNTS / bin duration
ERROR = ERROR / bin duration
# Retrieve region from DSS and make list of all pixels within it.
# Calculate the total area of the region.
if (use measured spectrum) {
# Rebin spectrum if specified by user.
} else {
# Make a 'natural' or model spectrum.
}
# Obtain the effective area as a function of energy.
# Calculate the correction for GTIs as a function of time.
tdx = calcGtiCorrection()
# Calculate the other time-dependent exposure factors (cosmic rays and
# dead time).
tdx = tdx * calcOtherTimeDepCorrns()
# Now calculate spatial contributions.
if (numCcdsInRegion == 1) {
if (withspecset) {
# Routine 00: 1 ccd/node, spectrum supplied.
(numerator, denominator) = calcSpatial00(ccd, node, imageModel, etc)
} else {
# Routine 01: 1 ccd/node, model spectrum used.
(numerator, denominator) = calcSpatial01(ccd, node, imageModel, etc)
}
FRACEXP = tdx * numerator / denominator
} else {
if (withspecset) {
# Routine 10: many ccds/nodes, spectrum supplied.
(numerator, denominator) = calcSpatial10(imageModel, etc)
FRACEXP = tdx * numerator / denominator
} else {
# Routine 11: many ccds/nodes, model spectrum used.
FRACEXP = calcFracexp11(tdx, imageModel, etc)
}
}
# Correct for external-source OOTEs.
# Correct for binning noise.
if (applyCorrections) {
foreach (timebin) {
if (FRACEXP(timebin)) > 0) {
RATE[timebin] = RATE[timebin] / FRACEXP(t)
ERROR[timebin] = ERROR[timebin] / FRACEXP(t)
} else {
RATE[timebin] = setToNull
ERROR[timebin] = setToNull
}
}
}
} # end foreach (srcTS, bkgTS)
# Normalise bkgTS to src collection area.
if (subtractBkg) {
foreach (timebin) {
RATE[timebin] = RATE[timebin] - BACKV[timebin]
ERROR[timebin] = sqrt(ERROR[timebin] * ERROR[timebin]
+ bkgERROR[timebin] * bkgERROR[timebin])
}
}
# Delete bkg FRACEXP column from output file.
# Release output file.
end subroutine lccorr
The workings of the four `spatial' routines are very similar to each other. To serve as an example, the first is laid out below.
subroutine calcSpatial00(ccd, node, imageModel, etc)
# startEbin = first Ebin with centre energy greater than TS minimum E.
# stopEbin = last Ebin with centre energy less than TS maximum E.
numerator = 0.0
denominator = 0.0
if (imageModel eq 'point') { # eqn 33, 34 of SSC-LUX-TN-0053
foreach Ebin (startEbin .. stopEbin) {
call CAL_getPSFmap(PSFmap)
expr = 0.0
foreach pixel (listOfPixels) {
call CAL_getFilterTransmission(filterTrans)
totalQe = 0.0
foreach patternID (0 .. 31) {
qeOfPattern = CAL_getQuantumEfficiency(patternID etc)
totalQe = totalQe + qeOfPattern
} # next patternID
expr = expr + filterTrans * totalQe * psf(srcX, srcY, pixelX, pixelY, PSFmap)
} # next pixel
call CAL_getVignettingFactor(vigFactor)
expr = expr * pixelArea * vigFactor
denominator = denominator + spectrum(Ebin) * dE / expr
numerator = numerator + spectrum(Ebin) * dE
} # next Ebin
} elsif (imageModel eq 'uniform') { # eqn 35, 36 of SSC-LUX-TN-0053
foreach Ebin (startEbin..stopEbin) {
expr = 0.0
foreach pixel (listOfPixels) {
call CAL_getFilterTransmission(filterTrans)
totalQe = 0.0
foreach patternID (0 .. 31) {
qeOfPattern = CAL_getQuantumEfficiency(patternID etc)
totalQe = totalQe + qeOfPattern
} # next patternID
call CAL_getVignettingFactor(vigFactor)
expr = expr + filterTrans * totalQe * vigFactor
} # next pixel
denominator = denominator + spectrum(Ebin) * dE / expr
numerator = numerator + spectrum(Ebin) * dE
} # next Ebin
denominator = denominator * regionArea
}
end subroutine calcSpatial00