Most of the code for this ticket was ported/adapted from
which was used to implement the models in Astier+19 (1905.08677, "The Shape of the Photon Transfer Curve in CCDs") to fit for the covariances from flat fields, in particular, the full model of Equation 20 (ptc.py already implemented the approximation in Equation 14 of that paper, to fit the model of that equation to the PTC curve (Cov00 (var) vs signal)). Pierre Astier kindly provided some slides to explain the structure of the code above: pipeline_astier_shape_of_ptc.pdf
For this ticket, if doCovariancesAstier=True an array with tags (numpy.recarray, referred to as "tuple" in Pierre's code) is calculated and returned by computeCovariancesAstier from the input flat pairs. The array contains the covariances at all signal levels and all lags up to a maximum lag/range. Then, the function fitData uses this array to fit the data to the full model:
The function plotCovariancesAstier reproduces some of the plots in Astier+19 (the style of a few of them changed to match the plotting style already in ptc.py). The plotting routines are in the file astierCovPtcPlots.py, fitData, in astierCovPtcFit.py, and the file astierCovFitParameters.py contains external code that aids in the fitting process (even external to Pierre, as evidenced by the disclaimer at the top of that file, which I chose to leave so that we can decide what to do with this external package). The file astierCovPtcUtils.py has some helper routines.
I implemented some unit tests in test/test_ptc.py under test_covAstier. The tests check that the gain is the same as the input gain from mock data, that the covariances via FFT (as it is in MeasurePhotonTransferCurveTask when doCovariancesAstier=True) are the same as calculated in real space, and that Cov[0, 0] (i.e., the variances) are consistent with the variances calculated with the standard method (when doCovariancesAstier=false).
I created some mock data to run the code (adding some brighter-fatter via the power-law model in GalSim); the resulting plots are here (the input gain of 0.75 e/ADU is recovered; same data for all channels): PTC_COVS_ASTIER_det0.pdf
There are some other things that still could be done (e.g., testing the code with on-sky or lab data, refactoring the plotting code of ptc.py and astierCovPtcPlots.py to put them in the same place, as well as the code in utils.py and astierCovPtcUtils.py, implementing the correction for the bias induced in the variance due to sigma-clipping (discussed above and mentioned in the paper in section 5.2, but we have decided that we will do this later:
DM-24762), etc), but with John Swinbank we decided that the code should be reviewed at the current state (as it is already a lot of code) and that further work should be tackled in other tickets.