From 5dd6d09b42471c30a6fa4580ae944d81baabcab7 Mon Sep 17 00:00:00 2001 From: mdlpstsci Date: Thu, 21 Nov 2024 13:00:54 -0500 Subject: [PATCH] HLA-1133/1241/1288 Implementation of ACS Serial CTE correction (#609) * Layered all the ACS Serial CTE changes onto "The Reunion" HSTCAL. This was done manually to view all the differences between the structure of the previous and "The Reunion" versions of HSTCAL. * Clean up and streamline of routines in dopcte-gen2.c. * Enable the parallel CTE correction * Added the constraint the serial CTE correction only applies to post-SM4 data. * Update date associated with version of code. * Checking in cleaned up files to support serial CTE correction. Renamed dopcte-gen2.c to dpcte-gen3.c. * Moved a message regarding no serial CTE correction for subarrays to a more appropriate location in the code. * Added the comment to the Dates file that the serial CTE implementation is NOT backwards compatible. WFC RAW files will need to have the PCTETAB keyword in the primary header updated with the name of the update calibration file which supports both parallel and serial CTE corrections. * Small clean-up * Fixed some text alignment issues and typos. --- ctegen2/ctegen2.h | 7 +- ctegen2/ctehelpers.c | 492 ++++++++++++++--------- pkg/acs/Dates | 11 + pkg/acs/History | 10 + pkg/acs/Updates | 11 + pkg/acs/include/acs.h | 2 +- pkg/acs/include/acsversion.h | 14 +- pkg/acs/lib/CMakeLists.txt | 2 +- pkg/acs/lib/acscte/docte.c | 501 ++++++++++++++++-------- pkg/acs/lib/acscte/dopcte-gen2.c | 498 ------------------------ pkg/acs/lib/acscte/dopcte-gen3.c | 649 +++++++++++++++++++++++++++++++ pkg/acs/lib/acscte/dopcte.c | 230 +---------- 12 files changed, 1354 insertions(+), 1073 deletions(-) delete mode 100644 pkg/acs/lib/acscte/dopcte-gen2.c create mode 100644 pkg/acs/lib/acscte/dopcte-gen3.c diff --git a/ctegen2/ctegen2.h b/ctegen2/ctegen2.h index 214a8f59d..34cf52f58 100644 --- a/ctegen2/ctegen2.h +++ b/ctegen2/ctegen2.h @@ -75,6 +75,9 @@ typedef struct { char descrip2[256]; /*descrip from table row, not read in for cte purposes*/ char cte_name[256]; /*name of cte algorithm */ char cte_ver[256]; /*version of algorithm */ + + // Serial CTE FITS extension information only + char ccdamp[2]; /* ID of specific amp for the serial CTE correction */ } CTEParamsFast; int inverseCTEBlurWithRowMajorIput(const SingleGroup * rsz, SingleGroup * rsc, const SingleGroup * trapPixelMap, CTEParamsFast * cte); @@ -102,9 +105,9 @@ void initCTEParamsFast(CTEParamsFast * pars, const unsigned _nTraps, const unsig int allocateCTEParamsFast(CTEParamsFast * pars); void freeCTEParamsFast(CTEParamsFast * pars); -int populateImageFileWithCTEKeywordValues(SingleGroup *group, CTEParamsFast *pars); +int populateImageFileWithCTEKeywordValues(SingleGroup *group, CTEParamsFast *pars, char * corrType); int getCTEParsFromImageHeader(SingleGroup * input, CTEParamsFast * params); -int loadPCTETAB(char *filename, CTEParamsFast * params); +int loadPCTETAB(char *filename, CTEParamsFast * params, int extn, Bool skipLoadPrimary); void ctewarn (char *message); void cteerror (char *message); void ctemessage (char *message); diff --git a/ctegen2/ctehelpers.c b/ctegen2/ctehelpers.c index e19dcabdd..63ac99af8 100644 --- a/ctegen2/ctehelpers.c +++ b/ctegen2/ctehelpers.c @@ -70,12 +70,11 @@ void initCTEParamsFast(CTEParamsFast * pars, const unsigned _nTraps, pars->qlevq_data = NULL; pars->dpdew_data = NULL; - pars->rprof = NULL; /*differential trail profile as image*/ - pars->cprof = NULL; /*cummulative trail profile as image*/ + pars->rprof = NULL; /* differential trail profile as image */ + pars->cprof = NULL; /* cummulative trail profile as image */ *pars->cte_name='\0'; *pars->cte_ver='\0'; - *pars->descrip2='\0'; } int allocateCTEParamsFast(CTEParamsFast * pars) @@ -177,185 +176,266 @@ MLS 2015: read in the CTE parameters from the PCTETAB file Jan 19, 2016: MLS Updated to check for existence of PCTENSMD in raw science header +MDD Sept 2023: Implementation to support a PCTETAB which is an + update to the Generation 2 CTE correction. This file will + now contain both parallel(y) and serial(x) CTE information. + Extensions 1-4 contain the parallel CTE correction data. The + remaining extensions contain the serial CTE correction, one + full set of corrections (QPROF, SCLBYCOL, RPROF, and CPROF) + per amplifier. */ -int loadPCTETAB (char *filename, CTEParamsFast *pars) { - /* Read the cte parameters from the reference table PCTETAB - - These are taken from the PCTETAB global header: - CTE_NAME - name of cte algorithm - CTE_VER - version number of cte algorithm - CTEDATE0 - date of instrument installation in HST, in fractional years - CTEDATE1 - reference date of CTE model pinning, in fractional years - - PCTETLEN - max length of CTE trail - PCTERNOI - readnoise amplitude and clipping level - PCTESMIT - number of iterations used in CTE forward modeling - PCTESHFT - number of iterations used in the parallel transfer - PCTENSMD - readnoise mitigation algorithm - PCTETRSH - over-subtraction threshold - - The table has 4 extensions: - -Filename: wfc3_cte.fits -No. Name Type Cards Dimensions Format -0 PRIMARY PrimaryHDU 21 () -1 QPROF BinTableHDU 16 nTraps>R x 3C ['i', 'i', 'i'] -2 SCLBYCOL BinTableHDU 20 nScaleTableColumns> x 5C ['i', 'e', 'e', 'e', 'e'] -3 RPROF ImageHDU 12 (nTraps>, 100) float32 -4 CPROF ImageHDU 12 (nTraps>, 100) float32 - - */ +int loadPCTETAB (char *filename, CTEParamsFast *pars, int extn, Bool skipLoadPrimary) { +/* Read the cte parameters from the reference table PCTETAB + + These are taken from the PCTETAB global header (aka primary header): + CTE_NAME - name of cte algorithm + CTE_VER - version number of cte algorithm + PCTERNOI - readnoise amplitude and clipping level + PCTENSMD - readnoise mitigation algorithm + PCTETRSH - over-subtraction threshold + FIXROCR - account for cosmic ray over-subtraction? + + These values are now amp-dependent for the *serial* CTE correction. This + means that they are read from the QPROF extension header of the amp being + processed. In contrast, there is only one set for the *parallel* CTE + correction which applies to all amps. The parallel values are also read + from the QPROF extension header where the parallel correction is defined + as discussed below. + CTEDATE0 - date of instrument installation in HST, in fractional years (MJD) + CTEDATE1 - reference date of CTE model pinning, in fractional years (MJD) + PCTENFOR - number of iterations used in CTE forward modeling + PCTENPAR - number of iterations used in the parallel transfer + PCTETLEN - max length of CTE trail + + The parallel/serial PCTETAB which supports the latest CTE algorithm + CTE_VER = '3.0 ' + CTE_NAME= 'Par/Serial PixelCTE 2023' + + The updated PCTETAB has a primary HDU and 20 extensions. Extensions 1-4 + apply to the parallel CTE and the remaining extensions apply to the serial CTE. + +No. Name Ver Type Cards Dimensions Format + 0 PRIMARY 1 PrimaryHDU 75 () + 1 QPROF 1 BinTableHDU 21 9999R x 4C [I, J, E, 20A] + 2 SCLBYCOL 1 BinTableHDU 25 8192R x 6C [I, E, E, E, E, 20A] + 3 RPROF 1 ImageHDU 21 (9999, 100) float32 + 4 CPROF 1 ImageHDU 21 (9999, 100) float32 + 5 QPROF 2 BinTableHDU 22 9999R x 4C [I, J, E, 20A] + 6 SCLBYCOL 2 BinTableHDU 26 8192R x 6C [I, E, E, E, E, 20A] + 7 RPROF 2 ImageHDU 22 (9999, 100) float32 + 8 CPROF 2 ImageHDU 22 (9999, 100) float32 + 9 QPROF 3 BinTableHDU 22 9999R x 4C [I, J, E, 20A] + 10 SCLBYCOL 3 BinTableHDU 26 8192R x 6C [I, E, E, E, E, 20A] + 11 RPROF 3 ImageHDU 22 (9999, 100) float32 + 12 CPROF 3 ImageHDU 22 (9999, 100) float32 + 13 QPROF 4 BinTableHDU 22 9999R x 4C [I, J, E, 20A] + 14 SCLBYCOL 4 BinTableHDU 26 8192R x 6C [I, E, E, E, E, 20A] + 15 RPROF 4 ImageHDU 22 (9999, 100) float32 + 16 CPROF 4 ImageHDU 22 (9999, 100) float32 + 17 QPROF 5 BinTableHDU 22 9999R x 4C [I, J, E, 20A] + 18 SCLBYCOL 5 BinTableHDU 26 8192R x 6C [I, E, E, E, E, 20A] + 19 RPROF 5 ImageHDU 22 (9999, 100) float32 + 20 CPROF 5 ImageHDU 22 (9999, 100) float32 + + The "extn" parameter indicates the starting extension for the QPROF table to be read + as this routine will be called multiple times - first for the parallel CTE correction + and then for the multiple sets (one set per amp) of serial CTE corrections. The serial + CTE is performed first. A "set" here means QPROF, SCLBYCOL, RPROF, and CPROF extensions. + Parallel: starting extension 1 (extensions 1 - 4) + Serial Amp A: starting extension 5 (extensions 5 - 8) + Serial Amp B: starting extension 9 (extensions 9 - 12) + Serial Amp C: starting extension 13 (extensions 13 - 16) + Serial Amp D: starting extension 17 (extensions 17 - 20) +*/ extern int status; /* variable for return status */ /* VARIABLE FOR FILENAME + EXTENSION NUMBER. */ char filename_wext[strlen(filename) + 4]; - /* NAMES OF DATA COLUMNS WE WANT FROM THE FILE, DATA WILL BE STORED IN THE PARS STRUCTURE */ + /* QPROF table - the last column is a description string */ const char wcol[] = "W"; const char qlevq[] = "QLEV_Q"; const char dpdew[] = "DPDE_W"; + + /* SCLBYCOL table - the last column is a description string */ const char iz[] = "IZ"; const char sens512[]= "SENS_0512"; const char sens1024[] = "SENS_1024"; const char sens1536[] = "SENS_1536"; const char sens2048[] = "SENS_2048"; - /* HSTIO VARIABLES */ - Hdr hdr_ptr; - initHdr(&hdr_ptr); - /* LOAD PRIMARY HEADER */ - if (LoadHdr(filename, &hdr_ptr)) { - sprintf(MsgText,"(pctecorr) Error loading header from %s",filename); - cteerror(MsgText); - status = OPEN_FAILED; - return status; - } - - /* GET CTE_NAME KEYWORD */ - if (GetKeyStr (&hdr_ptr, "CTE_NAME", NO_DEFAULT, "", pars->cte_name, SZ_CBUF)) { - cteerror("(pctecorr) Error reading CTE_NAME keyword from PCTETAB"); - status = KEYWORD_MISSING; - return status; - } - sprintf(MsgText,"\nCTE_NAME: %s",pars->cte_name); - trlmessage(MsgText); - - /* GET VERSION NUMBER */ - if (GetKeyStr(&hdr_ptr, "CTE_VER", NO_DEFAULT, "", pars->cte_ver, SZ_CBUF)) { - cteerror("(pctecorr) Error reading CTE_VER keyword from PCTETAB"); - status = KEYWORD_MISSING; - return status; - } - sprintf(MsgText,"CTE_VER: %s",pars->cte_ver); - trlmessage(MsgText); - - /* GET DATE OF INSTRUMENT INSTALLATION IN HST */ - if (GetKeyDbl(&hdr_ptr, "CTEDATE0", NO_DEFAULT, -999, &pars->cte_date0)) { - cteerror("(pctecorr) Error reading CTEDATE0 keyword from PCTETAB"); - status = KEYWORD_MISSING; - return status; - } - - sprintf(MsgText,"CTEDATE0: %g",pars->cte_date0); - trlmessage(MsgText); + /* Read in the primary header keywords */ + if (!skipLoadPrimary) + { + sprintf(MsgText, "(ctehelpers) Reading PRIMARY. cte_name: %s skipLoadPrimary: %c\n", pars->cte_name, skipLoadPrimary); + /* HSTIO VARIABLES */ + Hdr hdr_ptr; + initHdr(&hdr_ptr); + /* LOAD PRIMARY HEADER */ + if (LoadHdr(filename, &hdr_ptr)) { + sprintf(MsgText,"(pctecorr) Error loading header from %s",filename); + cteerror(MsgText); + status = OPEN_FAILED; + return status; + } - /* GET REFRENCE DATE OF CTE MODEL PINNING */ - if (GetKeyDbl(&hdr_ptr, "CTEDATE1", NO_DEFAULT, -999, &pars->cte_date1)) { - cteerror("(pctecorr) Error reading CTEDATE1 keyword from PCTETAB"); - status = KEYWORD_MISSING; - return status; + /* GET CTE_NAME KEYWORD */ + if (GetKeyStr (&hdr_ptr, "CTE_NAME", NO_DEFAULT, "", pars->cte_name, SZ_CBUF)) { + cteerror("(pctecorr) Error reading CTE_NAME keyword from PCTETAB"); + status = KEYWORD_MISSING; + return status; + } + sprintf(MsgText,"\nCTE_NAME: %s",pars->cte_name); + trlmessage(MsgText); + + /* GET VERSION NUMBER */ + if (GetKeyStr(&hdr_ptr, "CTE_VER", NO_DEFAULT, "", pars->cte_ver, SZ_CBUF)) { + cteerror("(pctecorr) Error reading CTE_VER keyword from PCTETAB"); + status = KEYWORD_MISSING; + return status; + } + sprintf(MsgText,"CTE_VER: %s",pars->cte_ver); + trlmessage(MsgText); + + /* GET READ NOISE CLIPPING LEVEL */ + if (GetKeyDbl(&hdr_ptr, "PCTERNOI", NO_DEFAULT, -999, &pars->rn_amp)) { + cteerror("(pctecorr) Error reading PCTERNOI keyword from PCTETAB"); + status = KEYWORD_MISSING; + return status; + } + sprintf(MsgText,"PCTERNOI: %f",pars->rn_amp); + trlmessage(MsgText); + + /* GET READ NOISE MITIGATION ALGORITHM*/ + if (GetKeyInt(&hdr_ptr, "PCTENSMD", NO_DEFAULT, -999, &pars->noise_mit)) { + cteerror("(pctecorr) Error reading PCTENSMD keyword from PCTETAB"); + status = KEYWORD_MISSING; + return status; + } + sprintf(MsgText,"PCTENSMD: %d",pars->noise_mit); + trlmessage(MsgText); + + /* GET OVER SUBTRACTION THRESHOLD */ + if (GetKeyDbl(&hdr_ptr, "PCTETRSH", NO_DEFAULT, -999, &pars->thresh)) { + cteerror("(pctecorr) Error readsng PCTETRSH keyword from PCTETAB"); + status = KEYWORD_MISSING; + return status; + } + sprintf(MsgText,"PCTETRSH: %g",pars->thresh); + trlmessage(MsgText); + + /* FIX THE READOUT CR'S? */ + if (GetKeyInt(&hdr_ptr, "FIXROCR", NO_DEFAULT, -999, &pars->fix_rocr)){ + cteerror("(pctecorr) Error reading FIXROCR keyword from PCTETAB"); + status = KEYWORD_MISSING; + return status; + } + sprintf(MsgText,"FIXROCR: %d",pars->fix_rocr); + trlmessage(MsgText); + + /* DONE READING STUFF FROM THE PRIMARY HEADER */ + freeHdr(&hdr_ptr); } - sprintf(MsgText,"CTEDATE1: %g",pars->cte_date1); - trlmessage(MsgText); - - /* READ MAX LENGTH OF CTE TRAIL */ - if (GetKeyInt(&hdr_ptr, "PCTETLEN", NO_DEFAULT, -999, &pars->cte_len)) { - cteerror("(pctecorr) Error reading PCTETLEN keyword from PCTETAB"); - status = KEYWORD_MISSING; - return status; - } + /* + Read in the remaining keywords necessary for proper processing and are + amp-dependent from the associated "extn". + + The variable extn contains the number of the FITS extension which is the starting + extension for the "set" of qprof/sclbycol/rprof/cprof data. The values for the + starting extension of a new set are: 5, 9, 13, 17, and 1. For example, the starting + serial CTE correct for Amp A as seen from the description of the PCTETAB earlier + in this file: + Extension EXTNAME EXTVER + 5 QPROF 2 + 6 SCLBYCOL 2 + 7 RPROF 2 + 8 CPROF 2 + */ - sprintf(MsgText,"PCTETLEN: %d",pars->cte_len); - trlmessage(MsgText); + /*****************************************************************************/ + /* READ DATA FROM THE SPECIFIED QPROF TABLE - EXTENSIONS 5, 9, 13, 17, and 1 */ + sprintf(filename_wext, "%s[%i]", filename, extn); - /* GET READ NOISE CLIPPING LEVEL */ - if (GetKeyDbl(&hdr_ptr, "PCTERNOI", NO_DEFAULT, -999, &pars->rn_amp)) { - cteerror("(pctecorr) Error reading PCTERNOI keyword from PCTETAB"); - status = KEYWORD_MISSING; + sprintf(MsgText,"Opening %s to read QPROF table.",filename_wext); + ctemessage(MsgText); + /* OPEN PARAMETERS FILE TO QPROF EXTENSION */ + IRAFPointer tbl_ptr = c_tbtopn(filename_wext, IRAF_READ_ONLY, 0); // xtables table pointer + if (c_iraferr()) { + sprintf(MsgText,"(pctecorr) Error opening %s with xtables",filename_wext); + cteerror(MsgText); + status = OPEN_FAILED; + c_tbtclo(tbl_ptr); return status; } - sprintf(MsgText,"PCTERNOI: %f",pars->rn_amp); - trlmessage(MsgText); - - /* GET NUMBER OF ITERATIONS USED IN FORWARD MODEL */ - if (GetKeyInt(&hdr_ptr, "PCTENFOR", NO_DEFAULT, -999, &pars->n_forward)) { - cteerror("(pctecorr) Error reading PCTENFOR keyword from PCTETAB"); - status = KEYWORD_MISSING; - return status; - } - sprintf(MsgText,"PCTERNFOR: %d",pars->n_forward); - trlmessage(MsgText); + /* + First get the keywords necessary for processing from the extension header + All the extension headers pertaining to an amp correction have the same + values for these keywords, so they only have to be read from the QPROF + extension. + */ - /* GET NUMBER OF ITERATIONS USED IN PARALLEL TRANSFER*/ - if (GetKeyInt(&hdr_ptr, "PCTENPAR", NO_DEFAULT, -999, &pars->n_par)) { - cteerror("(pctecorr) Error reading PCTENPAR keyword from PCTETAB"); - status = KEYWORD_MISSING; + /* GET DATE SERIAL CTE BECAME RELEVANT */ + pars->cte_date0 = c_tbhgtd(tbl_ptr, "CTEDATE0"); + if (status = c_iraferr()) { + sprintf(MsgText,"(pctecorr) Error reading CTEDATE0 from extension %d.", extn); + cteerror(MsgText); + c_tbtclo(tbl_ptr); return status; } + sprintf(MsgText,"CTEDATE0: %g",pars->cte_date0); + trlmessage(MsgText); - sprintf(MsgText,"PCTERNPAR: %d",pars->n_par); - trlmessage(MsgText); - - /* GET READ NOISE MITIGATION ALGORITHM*/ - if (GetKeyInt(&hdr_ptr, "PCTENSMD", NO_DEFAULT, -999, &pars->noise_mit)) { - cteerror("(pctecorr) Error reading PCTENSMD keyword from PCTETAB"); - status = KEYWORD_MISSING; + /* GET REFERENCE DATE OF CTE MODEL PINNING */ + pars->cte_date1 = c_tbhgtd(tbl_ptr, "CTEDATE1"); + if (status = c_iraferr()) { + sprintf(MsgText,"(pctecorr) Error reading CTEDATE1 from extension %d.", extn); + cteerror(MsgText); + c_tbtclo(tbl_ptr); return status; } - sprintf(MsgText,"PCTENSMD: %d",pars->noise_mit); - trlmessage(MsgText); + sprintf(MsgText,"CTEDATE1: %g",pars->cte_date1); + trlmessage(MsgText); - /* GET OVER SUBTRACTION THRESHOLD */ - if (GetKeyDbl(&hdr_ptr, "PCTETRSH", NO_DEFAULT, -999, &pars->thresh)) { - cteerror("(pctecorr) Error reading PCTETRSH keyword from PCTETAB"); - status = KEYWORD_MISSING; + /* READ MAX LENGTH OF CTE TRAIL */ + pars->cte_len = c_tbhgti(tbl_ptr, "PCTETLEN"); + if (status = c_iraferr()) { + sprintf(MsgText,"(pctecorr) Error reading PCTETLEN from extension %d.", extn); + cteerror(MsgText); + c_tbtclo(tbl_ptr); return status; } - sprintf(MsgText,"PCTETRSH: %g",pars->thresh); - trlmessage(MsgText); + sprintf(MsgText,"PCTETLEN: %d",pars->cte_len); + trlmessage(MsgText); - /*FIX THE READOUT CR'S? */ - if (GetKeyInt(&hdr_ptr, "FIXROCR", NO_DEFAULT, -999, &pars->fix_rocr)){ - cteerror("(pctecorr) Error reading FIXROCR keyword from PCTETAB"); - status = KEYWORD_MISSING; + /* GET NUMBER OF ITERATIONS USED IN FORWARD MODEL */ + pars->n_forward = c_tbhgti(tbl_ptr, "PCTENFOR"); + if (status = c_iraferr()) { + sprintf(MsgText,"(pctecorr) Error reading PCTENFOR from extension %d.", extn); + cteerror(MsgText); + c_tbtclo(tbl_ptr); return status; } - sprintf(MsgText,"FIXROCR: %d",pars->fix_rocr); - trlmessage(MsgText); - - - /* DONE READING STUFF FROM THE PRIMARY HEADER */ - freeHdr(&hdr_ptr); - /****************************************************************************/ - /* READ DATA FROM FIRST TABLE EXTENSIONS */ - sprintf(filename_wext, "%s[%i]", filename, 1); + sprintf(MsgText,"PCTENFOR: %d",pars->n_forward); + trlmessage(MsgText); - /* OPEN PARAMETERS FILE TO EXTENSION NUMBER 1 */ - IRAFPointer tbl_ptr = c_tbtopn(filename_wext, IRAF_READ_ONLY, 0); // xtables table pointer - if (c_iraferr()) { - sprintf(MsgText,"(pctecorr) Error opening %s with xtables",filename_wext); + /* GET NUMBER OF ITERATIONS USED IN TRANSFER */ + pars->n_par = c_tbhgti(tbl_ptr, "PCTENPAR"); + if (status = c_iraferr()) { + sprintf(MsgText,"(pctecorr) Error reading PCTENPAR from extension %d.", extn); cteerror(MsgText); - status = OPEN_FAILED; c_tbtclo(tbl_ptr); return status; } + sprintf(MsgText,"PCTENPAR: %d",pars->n_par); + trlmessage(MsgText); /* READ DATA FROM TABLE */ + /* get column pointer for w */ IRAFPointer w_ptr = c_tbcfnd1_retPtr(tbl_ptr, wcol); if (c_iraferr() || !w_ptr) { @@ -383,7 +463,6 @@ No. Name Type Cards Dimensions Format return status; } - // LOOP OVER TABLE ROWS UP TO SIZE TRAPS int ctraps = 0; // actual usable traps, i.e. see if more traps were added to reference file {unsigned j; @@ -433,14 +512,17 @@ No. Name Type Cards Dimensions Format trlmessage(MsgText); */ - /* CLOSE CTE PARAMETERS FILE FOR EXTENSION 1*/ + /* CLOSE CTE PARAMETERS FILE FOR QPROF EXTENSION */ c_tbtClose((void*)&tbl_ptr); assert(!tbl_ptr); + /****************************************************************************/ /****************************************************************************/ - /* READ CTE SCALING DATA FROM SECOND TABLE EXTENSION */ - sprintf(filename_wext, "%s[%i]", filename, 2); + /* READ CTE SCALING DATA FROM THE SPECIFIED SCLBYCOL TABLE - EXTENSIONS 6, 10, 14, 18, and 2 */ + sprintf(filename_wext, "%s[%i]", filename, extn + 1); + sprintf(MsgText,"Opening %s to read SCLBYCOL table.",filename_wext); + ctemessage(MsgText); tbl_ptr = c_tbtopn(filename_wext, IRAF_READ_ONLY, 0); if (c_iraferr()) { sprintf(MsgText,"(pctecorr) Error opening %s with xtables",filename_wext); @@ -493,7 +575,6 @@ No. Name Type Cards Dimensions Format return status; } - /* read data from table */ /* loop over table rows */ {unsigned j; @@ -541,13 +622,19 @@ No. Name Type Cards Dimensions Format } */ - /* close CTE parameters file for extension 2*/ + /* close CTE parameters file for SCLBYCOL extension */ c_tbtClose((void*)&tbl_ptr); assert(!tbl_ptr); - /****************************************************************************/ - /* extension 3: differential trail profile as image */ - ctemessage("Reading in image from extension 3"); + /*********************************************************************************/ + /* extensions 7, 11, 15, 19, and 3 - RPROF : differential trail profile as image */ + /* Images are read by EXTNAM and EXTVER and not extension number */ + + /* Determine the extension version number for the RPROF/CPROF images */ + int extver = (extn / 4) + 1; + + sprintf(MsgText,"Reading in image from RPROF EXTVER %d.", extver); + ctemessage(MsgText); /* Get the coefficient images from the PCTETAB */ pars->rprof = malloc(sizeof(*pars->rprof)); @@ -556,16 +643,20 @@ No. Name Type Cards Dimensions Format trlerror (MsgText); return (status = 1); } + initFloatHdrData(pars->rprof); pars->rprof->data.storageOrder = COLUMNMAJOR; - if (getFloatHD (filename, "RPROF", 1, pars->rprof)){ + if (getFloatHD (filename, "RPROF", extver, pars->rprof)){ return (status=1); } - /****************************************************************************/ - /* ext number 4 : cummulative trail profile as image */ - ctemessage("Reading in image from extension 4"); + /*********************************************************************************/ + /* extensions 8, 12, 16, 20, and 4 - CPROF : cummulative trail profile as image */ + /* This image has the same EXTVER as the RPROF. */ + + sprintf(MsgText,"Reading in image from CPROF EXTVER %d.", extver); + ctemessage(MsgText); pars->cprof = malloc(sizeof(*pars->cprof)); if (pars->cprof == NULL){ @@ -577,7 +668,7 @@ No. Name Type Cards Dimensions Format /* Get the coefficient images from the PCTETAB */ initFloatHdrData (pars->cprof); pars->cprof->data.storageOrder = COLUMNMAJOR; - if (getFloatHD (filename, "CPROF", 1, pars->cprof)){ + if (getFloatHD (filename, "CPROF", extver, pars->cprof)){ return (status=1); } @@ -587,38 +678,76 @@ No. Name Type Cards Dimensions Format /* Some CTE parameters will likely end up in the image headers so users can tune them. This will first check whether these parameters are set in the image header. If they - are then the values found there will be used instead of the values read from + are, then the values found there will be used instead of the values read from the PCTETAB. If the values are not set they will be populated with the values from the PCTETAB, namely these: - 'CTE_NAME':'pixelCTE 2012', #name of cte algorithm - 'CTE_VER':'1.0' , #version number of algorithm - 'CTEDATE0':54962.0, #date of instrument installation in HST in MJD - 'CTEDATE1':56173.0, #reference date of cte model pinning in MJd - 'PCTETLEN':60, #max length of CTE trail - 'PCTERNOI':3.25, #read noise amplitude, clipping limit - 'PCTENFOR':5 ,#number of iterations used in cte forward modeling - 'PCTENPAR':7 ,#number of iterations used in parallel transfer - 'PCTENSMD':0 ,#read noise mitigation algorithm + 'CTE_NAME':'Par/Serial pixelCTE 2023', #name of cte algorithm + 'CTE_VER':'3.0' , #version number of algorithm + 'CTEDATE0':nnnnn.0, #date of instrument installation in HST in MJD + 'CTEDATE1':nnnnn.0, #reference date of cte model pinning in MJd + 'PCTETLEN':nnn, #max length of CTE trail + 'PCTERNOI':n.nn, #read noise amplitude, clipping limit + 'PCTENFOR':n ,#number of iterations used in cte forward modeling + 'PCTENPAR':n ,#number of iterations used in parallel transfer + 'PCTENSMD':n ,#read noise mitigation algorithm 'PCTETRSH':-10.0 ,#over subtraction threshold, always use reference value 'FIXROCR' : 1, #set to 1 for true, fix the readout cr's */ -int populateImageFileWithCTEKeywordValues(SingleGroup *group, CTEParamsFast *pars) +int populateImageFileWithCTEKeywordValues(SingleGroup *group, CTEParamsFast *pars, char *corrType) { extern int status; - if ((status = updateKeyOrAddAsHistKeyStr(group->globalhdr,"CTE_NAME", pars->cte_name, "CTE algorithm name"))) - { - trlerror("(pctecorr) failed to update (or add as history) CTE_NAME keyword in image header"); - return status; + /* + HISTORY information written only one time + */ + if (strncmp(corrType, "parallel", 6) == 0) { + + if ((status = updateKeyOrAddAsHistKeyStr(group->globalhdr,"CTE_NAME", pars->cte_name, "CTE algorithm name"))) + { + trlerror("(pctecorr) failed to update (or add as history) CTE_NAME keyword in image header"); + return status; + } + + if ((status = updateKeyOrAddAsHistKeyStr(group->globalhdr,"CTE_VER", pars->cte_ver, "CTE algorithm version"))) + { + trlerror("(pctecorr) failed to update (or add as history) CTE_VER keyword in image header"); + return status; + } + + if ((status = updateKeyOrAddAsHistKeyStr(group->globalhdr,"CORRTYPE", corrType, "CTE correction type"))) + { + trlerror("(pctecorr) failed to update (or add as history) CORRTYPE keyword in image header"); + return status; + } + + if ((status = updateKeyOrAddAsHistKeyInt(group->globalhdr,"FIXROCR",pars->fix_rocr,"fix readout cosmic rays"))) + { + trlerror("(pctecorr) failed to update (or add as history) FIXROCR keyword in header"); + return status; + } + + if ((status = updateKeyOrAddAsHistKeyDouble(group->globalhdr, "PCTETRSH", pars->thresh,"CTE over subtraction threshold"))) + { + trlerror("(pctecorr) failed to update (or add as history) PCTETRSH keyword in header"); + return status; + } + + if ((status = updateKeyOrAddAsHistKeyDouble(group->globalhdr, "PCTEFRAC", pars->scale_frac, "CTE scaling factor"))) + { + trlerror("(pctecorr) failed to update (or add as history) PCTEFRAC to image header"); + return status; + } } - if ((status = updateKeyOrAddAsHistKeyStr(group->globalhdr,"CTE_VER", pars->cte_ver, "CTE algorithm version"))) - { - trlerror("(pctecorr) failed to update (or add as history) CTE_VER keyword in image header"); - return status; + if (strncmp(corrType, "serial", 6) == 0) { + if ((status = updateKeyOrAddAsHistKeyStr(group->globalhdr,"CORRTYPE", corrType, "CTE correction type"))) + { + trlerror("(pctecorr) failed to update (or add as history) CORRTYPE keyword in image header"); + return status; + } } if ((status = updateKeyOrAddAsHistKeyDouble(group->globalhdr, "CTEDATE0", pars->cte_date0,"Date of instrument installation"))) @@ -627,25 +756,12 @@ int populateImageFileWithCTEKeywordValues(SingleGroup *group, CTEParamsFast *par return status; } - if ((status = updateKeyOrAddAsHistKeyDouble(group->globalhdr, "PCTETRSH", pars->thresh,"CTE over subtraction threshold"))) - { - trlerror("(pctecorr) failed to update (or add as history) PCTETRSH keyword in header"); - return status; - } - - if ((status = updateKeyOrAddAsHistKeyDouble(group->globalhdr, "CTEDATE1", pars->cte_date1, "Date of CTE model pinning"))) { trlerror("(pctecorr) failed to update (or add as history) CTEDATE1 keyword in header"); return status; } - if ((status = updateKeyOrAddAsHistKeyDouble(group->globalhdr, "PCTEFRAC", pars->scale_frac, "CTE scaling factor"))) - { - trlerror("(pctecorr) failed to update (or add as history) PCTEFRAC to image header"); - return status; - } - if ((status = updateKeyOrAddAsHistKeyInt(group->globalhdr,"PCTETLEN",pars->cte_len,"max length of CTE trail"))) { trlerror("(pctecorr) failed to update (or add as history) PCTETLEN in header"); @@ -672,12 +788,6 @@ int populateImageFileWithCTEKeywordValues(SingleGroup *group, CTEParamsFast *par return status; } - if ((status = updateKeyOrAddAsHistKeyInt(group->globalhdr,"FIXROCR",pars->fix_rocr,"fix readout cosmic rays"))) - { - trlerror("(pctecorr) failed to update (or add as history) FIXROCR keyword in header"); - return status; - } - return HSTCAL_OK; } @@ -692,6 +802,8 @@ int getCTEParsFromImageHeader(SingleGroup *group, CTEParamsFast *pars) { int fix_rocr = 0; int noise_mit = 0; + trlmessage("Begin attempt to read CTE keywords from input image header. Missing CTE keywords"); + trlmessage("from input image header are not a fatal error as the values from PCTETAB will be used."); int tempStatus = HSTCAL_OK; /*check the PCTENSMD keyword in the header*/ tempStatus = GetKeyInt(group->globalhdr, "PCTENSMD", NO_DEFAULT, -999, &noise_mit); @@ -788,6 +900,8 @@ int getCTEParsFromImageHeader(SingleGroup *group, CTEParamsFast *pars) { else if (tempStatus == KEYWORD_MISSING) status = HSTCAL_OK; + trlmessage("End attempt to read CTE keywords from input image header."); + return HSTCAL_OK; } diff --git a/pkg/acs/Dates b/pkg/acs/Dates index 058017c35..6a9832033 100644 --- a/pkg/acs/Dates +++ b/pkg/acs/Dates @@ -1,3 +1,14 @@ + 07-May-2024 CALACS 10.4.0 Implementation of the "Parallel and Serial CTE correction". This + version of the CTE correction (Generation 3) is based upon the + algorithm of Generation 2, but includes the additional correction + for the serial direction which is amp-dependent. Previous + generations of CTE correction are no longer supported, and + obsolete code has been removed. Note the PCTETAB has changed + substantially to accommodate the amp-dependent serial corrections. + The serial CTE correction only applies to full-frame, post-SM4 data. + NOTE: This is NOT a backwards compatible change. WFC RAW files + will need to have the PCTETAB keyword updated with the Generation 3 + Parallel/Serial Pixel CTE 2023 file. 08-Feb-2022 CALACS 10.3.5 Update to the cosmic ray rejection algorithm as to the way the output ERR extension is computed for the CRJ file. The output ERR is now propagated from the usable input ERR extensions diff --git a/pkg/acs/History b/pkg/acs/History index 6e6aeea74..ec58f51d5 100644 --- a/pkg/acs/History +++ b/pkg/acs/History @@ -1,3 +1,13 @@ +### 07-May-2024 - MDD - Version 10.4.0 + Implementation of the "Parallel and Serial CTE correction". This + version of the CTE correction (Generation 3) is based upon the + algorithm of Generation 2, but includes the additional correction + for the serial direction which is amp-dependent. Previous + generations of CTE correction are no longer supported, and obsolete + code has been removed. Note the PCTETAB has changed substantially + to accommodate the amp-dependent serial corrections. + The serial CTE correction only applies to full-frame, post-SM4 data. + ### 08-Feb-2022 - MDD - Version 10.3.5 Update to the cosmic ray rejection algorithm as to the way the output ERR extension is computed for the CRJ file. The diff --git a/pkg/acs/Updates b/pkg/acs/Updates index 99f18530c..823c2b0e5 100644 --- a/pkg/acs/Updates +++ b/pkg/acs/Updates @@ -1,3 +1,14 @@ +Update for version 10.4.0 - 07-May-2024 (MDD) + pkg/acs/Dates + pkg/acs/History + pkg/acs/Updates + pkg/acs/include/acs.h + pkg/acs/lib/CMakeLists.txt + pkg/acs/lib/acscte/docte.c + pkg/acs/lib/acscte/dopcte-gen2.c -> pkg/acs/calace/acscte/dopcte-gen3.c + ctegen2/ctegen2.h + ctegen2/ctehelpers.c + Update for version 10.3.5 - 08-Feb-2022 (MDD) pkg/acs/calacs/Dates pkg/acs/calacs/History diff --git a/pkg/acs/include/acs.h b/pkg/acs/include/acs.h index 9b421ee4a..abbdc6507 100644 --- a/pkg/acs/include/acs.h +++ b/pkg/acs/include/acs.h @@ -83,7 +83,7 @@ void errchk (); /* HSTIO error check */ # define DEFAULT_OFFSET 3 -# define SM4MJD 54967 +# define SM4MJD 54962 /* May 11, 2009 - beginning of SM4 */ /* October 01, 2016: Date of first observation in Cycle 24. The new ACS subarray configurations validated for Cycle 24. */ # define CYCLE24 57662 diff --git a/pkg/acs/include/acsversion.h b/pkg/acs/include/acsversion.h index c624c55e4..7686d446f 100644 --- a/pkg/acs/include/acsversion.h +++ b/pkg/acs/include/acsversion.h @@ -2,15 +2,25 @@ #define INCL_ACSVERSION_H /* This string is written to the output primary header as CAL_VER. */ -#define ACS_CAL_VER "10.3.5 (08-Feb-2022)" -#define ACS_CAL_VER_NUM "10.3.5" +#define ACS_CAL_VER "10.4.0 (07-May-2024)" +#define ACS_CAL_VER_NUM "10.4.0" +/* This Generation of the CTE algorithm is obsolete. These strings + are only maintained in case a user has an older version of the + PCTETAB. */ /* name and version number of the CTE correction algorithm */ #define ACS_GEN1_CTE_NAME "PixelCTE 2012" #define ACS_GEN1_CTE_VER "3.3" +/* This Generation of the CTE algorithm is obsolete. These strings + are only maintained in case a user has an older version of the + PCTETAB. */ /* name and version number of generation 2 CTE correction algorithm */ #define ACS_GEN2_CTE_NAME "PixelCTE 2017" #define ACS_GEN2_CTE_VER "2.0" +/* name and version number of generation 3 CTE correction algorithm */ +#define ACS_GEN3_CTE_NAME "Par/Serial PixelCTE 2023" +#define ACS_GEN3_CTE_VER "3.0" + #endif /* INCL_ACSVERSION_H */ diff --git a/pkg/acs/lib/CMakeLists.txt b/pkg/acs/lib/CMakeLists.txt index e812fdcf4..7ee926f04 100644 --- a/pkg/acs/lib/CMakeLists.txt +++ b/pkg/acs/lib/CMakeLists.txt @@ -33,7 +33,7 @@ add_library(${PROJECT_NAME} SHARED acscte/acscte.c acscte/docte.c acscte/domaincte.c - acscte/dopcte-gen2.c + acscte/dopcte-gen3.c acscte/dopcte.c acscte/getcteflag.c acscte/getctesw.c diff --git a/pkg/acs/lib/acscte/docte.c b/pkg/acs/lib/acscte/docte.c index f85e5984a..f2936b369 100644 --- a/pkg/acs/lib/acscte/docte.c +++ b/pkg/acs/lib/acscte/docte.c @@ -1,8 +1,11 @@ -/* Do CTE loss correction for the current set of extensions. +/* The primary header will be updated with history information, and the calibration switches in the header will be reset from "PERFORM" to "COMPLETE". 12-Aug-2012 PLL: Separated PCTECORR from ACSCCD. + 14-Sep-2023 MDD: Updated to accommodate only the "Parallel/Serial PixelCTE 2023" + (aka Generation 3) correction. Code now applies a serial CTE correction for + full-frame data. */ # include @@ -27,7 +30,32 @@ static void PCTEMsg (ACSInfo *, int); static int OscnTrimmed (Hdr*, Hdr *); -#define SZ_CBUF 24 +/* + Typical order of processing is Chip 2 (Amps C and D) and then + Chip 1 (Amps A and B). See ctehelpers.c for full description of + PCTETAB which clarifies the order of the extensions for processing. + + Background information added here in August 2024 MDD: + Since February 2018 according to the ACS keyword rules which are + used to assign the default values of PERFORM or OMIT to the calibration + steps (e.g., DARKCORR, PCTECORR), no subarrays are to be processed with + the CTE correction during pipeline processing. However, the ACS team has + allowed for the possibility a user would want to apply the CTE correction + to subarray data by setting PCTECORR = PERFORM. In this case, ONLY the + parallel CTE correction will be applied to the data. + + At one point for post-SM4 data, the parallel CTE correction was allowed + to be applied to certain subarray modes. However, in February 2018, this + was changed to PCTECORR=OMIT for post-SM4 where SUBARRAY=T. As + PCTECORR=OMIT by default, except for specified modes, pre-SM4 subarrays + are also not processed with parallel CTE correction. +*/ + +/* CTE calibration data as stored in the PCTETAB */ +int SET_TO_PROCESS[] = {13, 17, 5, 9}; // Starting extension number of a set of qprof/sclbycol/rprof/cprof +# define AMPCALIBORDER "CDAB" + +#define SZ_CBUF 30 int getCTE_NAME(char * filename, char * cteName, int cteNameBufferLength); int DoCTE (ACSInfo *acs_info, const bool forwardModelOnly) { @@ -38,15 +66,18 @@ int DoCTE (ACSInfo *acs_info, const bool forwardModelOnly) { extern int status; - SingleGroup * x; /* used for both input and output */ + SingleGroup * x; /* used for both input and output */ ACSInfo * acs; /* hold a copy of the acs_info struct for each extension */ int option = 0; + int ampID; /* index amp A:0, B:1, etc. */ + int ampIDInCalib; /* index amp C:0, D:1, etc. in calibration file */ + char * amploc; /* pointer to amp character in AMPSORDER */ + char * amplocInCalib; /* pointer to amp character in calibration file */ + int numamps; - int i; /* loop index */ Bool subarray; int CCDHistory (ACSInfo *, Hdr *); - int doPCTEGen1 (ACSInfo *, SingleGroup *); - int doPCTEGen2 (ACSInfo *, CTEParamsFast * pars, SingleGroup *, const bool forwardModelOnly); + int doPCTEGen3 (ACSInfo *, CTEParamsFast * pars, SingleGroup *, const bool forwardModelOnly, char * corrType, char *ccdamp, int nthAmp, char *amploc, int ampID); int pcteHistory (ACSInfo *, Hdr *); int GetACSGrp (ACSInfo *, Hdr *); int OmitStep (int); @@ -60,6 +91,8 @@ int DoCTE (ACSInfo *acs_info, const bool forwardModelOnly) { void UFilename (char *, Hdr *); int GetCCDTab (ACSInfo *, int, int); int GetKeyBool (Hdr *, char *, int, Bool, Bool *); + int GetKeyStr (Hdr *, char *, int, char *, char *, int); + void parseWFCamps (char *acsamps, int chip, char *ccdamp); /*========================Start Code here =======================*/ @@ -67,9 +100,11 @@ int DoCTE (ACSInfo *acs_info, const bool forwardModelOnly) { x = (SingleGroup *) malloc(acs_info->nimsets * sizeof(SingleGroup)); acs = (ACSInfo *) malloc(acs_info->nimsets * sizeof(ACSInfo)); - for (i = 0; i < acs_info->nimsets; i++) { - initSingleGroup(&x[i]); - acs[i] = *acs_info; + { unsigned i; + for (i = 0; i < acs_info->nimsets; i++) { + initSingleGroup(&x[i]); + acs[i] = *acs_info; + } } if (acs_info->printtime) @@ -81,12 +116,14 @@ int DoCTE (ACSInfo *acs_info, const bool forwardModelOnly) { processing step functions and pass along modified input image from one step to the next... */ - for (i = 0; i < acs_info->nimsets; i++) { - getSingleGroup(acs[i].input, i+1, &x[i]); + { unsigned i; + for (i = 0; i < acs_info->nimsets; i++) { + getSingleGroup(acs[i].input, i+1, &x[i]); - if (hstio_err()) { - freeSingleGroup(&x[i]); - return (status = OPEN_FAILED); + if (hstio_err()) { + freeSingleGroup(&x[i]); + return (status = OPEN_FAILED); + } } } @@ -99,35 +136,37 @@ int DoCTE (ACSInfo *acs_info, const bool forwardModelOnly) { Probably redundant here but does not hurt. */ - for (i = 0; i < acs_info->nimsets; i++) { - if (GetACSGrp(&acs[i], &x[i].sci.hdr)) { - freeSingleGroup(&x[i]); - return (status); - } - - /* Has the CCD image been overscan trimmed? */ - if (OscnTrimmed (x[i].globalhdr, &x[i].sci.hdr)) { - freeSingleGroup (&x[i]); - return (status); - } + { unsigned i; + for (i = 0; i < acs_info->nimsets; i++) { + if (GetACSGrp(&acs[i], &x[i].sci.hdr)) { + freeSingleGroup(&x[i]); + return (status); + } - /* Read in keywords from primary header... */ + /* Has the CCD image been overscan trimmed? */ + if (OscnTrimmed (x[i].globalhdr, &x[i].sci.hdr)) { + freeSingleGroup (&x[i]); + return (status); + } - if (GetKeyBool(x[i].globalhdr, "SUBARRAY", NO_DEFAULT, 0, &subarray)) - return (status); + /* Read in keywords from primary header... */ - if (subarray) - acs[i].subarray = YES; - else - acs[i].subarray = NO; + if (GetKeyBool(x[i].globalhdr, "SUBARRAY", NO_DEFAULT, 0, &subarray)) + return (status); - /* For the CCD, update primary header keywords. - Also reset CRCORR if there's only one image set. - */ - /* Get values from tables, using same function used in ACSCTE. */ - if (GetCCDTab(&acs[i], x[i].sci.data.nx, x[i].sci.data.ny)) { - freeSingleGroup(&x[i]); - return (status); + if (subarray) + acs[i].subarray = YES; + else + acs[i].subarray = NO; + + /* For the CCD, update primary header keywords. + Also reset CRCORR if there's only one image set. + */ + /* Get values from tables, using same function used in ACSCTE. */ + if (GetCCDTab(&acs[i], x[i].sci.data.nx, x[i].sci.data.ny)) { + freeSingleGroup(&x[i]); + return (status); + } } } @@ -162,152 +201,307 @@ int DoCTE (ACSInfo *acs_info, const bool forwardModelOnly) { if (acs_info->pctecorr == PERFORM) { - PtrRegister ptrReg;//move this up later + PtrRegister ptrReg; initPtrRegister(&ptrReg); + + PtrRegister ptrParallelReg; + initPtrRegister(&ptrParallelReg); + + /* These variables will hold the serial(aka 'x') and the parallel(aka 'y') + CTE values. */ CTEParamsFast ctePars; + CTEParamsFast cteParallelPars; //NOTE: The char * below should be const but this would require a massive refactoring. char * cteTabFilename = acs->pcte.name; - //open PCTETAB and check for CTE algorithm version name: CTE_NAME + // Generation 1 CTE algorithm is obsolete. + // Generation 3 CTE algorithm is the Generation 2 algorithm with updates to apply + // not only a parallel correction, but also possibly an amp-based serial CTE correction. + // There is no longer any choice for the generation of the CTE algorithm. The correction + // is only "Generation 3" which is also know by its FITS keyword CTE_NAME + // "Par/Serial PixelCTE 2023" or the "Parallel and Serial CTE correction". + + // open PCTETAB and check for CTE algorithm version name: CTE_NAME char cteName[SZ_CBUF]; *cteName = '\0'; int cteAlgorithmGen = 0; int tmpStatus = getCTE_NAME(cteTabFilename, cteName, SZ_CBUF); if (tmpStatus == KEYWORD_MISSING) - cteAlgorithmGen = 1; + { + trlerror("(pctecorr) No CTE_NAME keyword found in the primary header of PCTETAB file."); + trlerror("(pctecorr) This version of PCTETAB is obsolete - use the 'Parallel and Serial PCTETAB' PCTETAB."); + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status = ERROR_RETURN); + } else if (tmpStatus) { freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); return (status = tmpStatus); } + else if (strcmp(cteName, ACS_GEN3_CTE_NAME) == 0) + { + cteAlgorithmGen = 3; + trlmessage("(pctecorr) Parallel and Serial PCTETAB file auto-detected."); + } else if (strcmp(cteName, ACS_GEN2_CTE_NAME) == 0) { - cteAlgorithmGen = 2; - trlmessage("(pctecorr) Generation 2 PCTETAB file auto-detected."); + trlerror("(pctecorr) Generation 2 PCTETAB file auto-detected."); + trlerror("(pctecorr) This version of PCTETAB is obsolete - use the 'Parallel and Serial PCTETAB' PCTETAB."); + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status = ERROR_RETURN); } else { - cteAlgorithmGen = 1; - trlmessage("(pctecorr) Generation 1 PCTETAB file auto-detected."); + trlerror("(pctecorr) Generation 1 PCTETAB file auto-detected."); + trlerror("(pctecorr) This version of PCTETAB is obsolete - use the 'Parallel and Serial PCTETAB' PCTETAB."); + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status = ERROR_RETURN); } + acs_info->cteAlgorithmGen = cteAlgorithmGen; - if (acs_info->cteAlgorithmGen && (acs_info->cteAlgorithmGen != cteAlgorithmGen)) + sprintf(MsgText, "(pctecorr) Reading CTE parameters from PCTETAB file: '%s'...", cteTabFilename); + trlmessage(MsgText); + if ((status = PutKeyStr(x[0].globalhdr, "PCTETAB", cteTabFilename, "CTE Correction Table"))) { - char msgBuffer[MSG_BUFF_LENGTH]; - sprintf(msgBuffer, "(pctecorr) Cmd line option '--ctegen %d' specified yet gen%d algorithm detected in PCTETAB (CTE_NAME).", - acs_info->cteAlgorithmGen, cteAlgorithmGen); - trlerror(msgBuffer); + trlerror("(pctecorr) failed to update PCTETAB keyword in image header"); freeOnExit(&ptrReg); - return (status = CAL_FILE_MISSING); + freeOnExit(&ptrParallelReg); + return status; } - else if (!acs_info->cteAlgorithmGen && (cteAlgorithmGen != 2)) + + unsigned nScaleTableColumns = N_COLUMNS_FOR_RAZ_CDAB_ALIGNED_IMAGE; + + /* Read the parallel CTE parameters */ + sprintf(MsgText,"(pctecorr) Collecting data for Correction Type: parallel.\n"); + trlmessage(MsgText); + + addPtr(&ptrParallelReg, &cteParallelPars, &freeCTEParamsFast); + initCTEParamsFast(&cteParallelPars, TRAPS, 0, 0, nScaleTableColumns, acs_info->nThreads); + cteParallelPars.refAndIamgeBinsIdenticle = True; + cteParallelPars.verbose = acs->verbose = 0 ? False : True; + + if ((status = allocateCTEParamsFast(&cteParallelPars))) { - char msgBuffer[MSG_BUFF_LENGTH]; - sprintf(msgBuffer, "(pctecorr) Gen%d algorithm detected in PCTETAB (CTE_NAME). Default is gen2, use '--ctegen %d' to override.", - cteAlgorithmGen, cteAlgorithmGen); - trlerror(msgBuffer); freeOnExit(&ptrReg); - return (status = CAL_FILE_MISSING); + freeOnExit(&ptrParallelReg); + return (status); } - acs_info->cteAlgorithmGen = cteAlgorithmGen; - if (acs_info->cteAlgorithmGen == 2) + int startOfSetInCalib = 1; + Bool skipLoadPrimary = False; + if ((status = loadPCTETAB(cteTabFilename, &cteParallelPars, startOfSetInCalib, skipLoadPrimary))) { - sprintf(MsgText, "(pctecorr) Reading CTE parameters from PCTETAB file: '%s'...", cteTabFilename); - trlmessage(MsgText); - if ((status = PutKeyStr(x[0].globalhdr, "PCTETAB", cteTabFilename, "CTE Correction Table"))) - { - trlerror("(pctecorr) failed to update PCTETAB keyword in image header"); - return status; - } - //Get parameters from PCTETAB reference file - addPtr(&ptrReg, &ctePars, &freeCTEParamsFast); - unsigned nScaleTableColumns = N_COLUMNS_FOR_RAZ_CDAB_ALIGNED_IMAGE; - initCTEParamsFast(&ctePars, TRAPS, 0, 0, nScaleTableColumns, acs_info->nThreads); - ctePars.refAndIamgeBinsIdenticle = True; - ctePars.verbose = acs->verbose = 0 ? False : True; - if ((status = allocateCTEParamsFast(&ctePars))) - { - freeOnExit(&ptrReg); - return (status); - } - - if ((status = loadPCTETAB(cteTabFilename, &ctePars))) - { - freeOnExit(&ptrReg); - return (status); - } + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status); + } - if ((status = getCTEParsFromImageHeader(&x[0], &ctePars))) - { - freeOnExit(&ptrReg); - return (status); - } + if ((status = getCTEParsFromImageHeader(&x[0], &cteParallelPars))) + { + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status); + } - ctePars.scale_frac = (acs->expstart - ctePars.cte_date0) / (ctePars.cte_date1 - ctePars.cte_date0); + cteParallelPars.scale_frac = (acs->expstart - cteParallelPars.cte_date0) / (cteParallelPars.cte_date1 - cteParallelPars.cte_date0); - if ((status = populateImageFileWithCTEKeywordValues(&x[0], &ctePars))) - { - freeOnExit(&ptrReg); - return (status); - } + /* + This routine writes information to the HISTORY portion of the output + primary header which includes KEYWORD values read from the PCTETAB. + Some of the KEYWORDs have different values for their use in the application + of the parallel and serial CTE corrections. Further, the serial CTE is + amp-dependent so the KEYWORDs are documented for each amp. - sprintf(MsgText, "(pctecorr) IGNORING read noise level PCTERNOI from PCTETAB: %f. Using amp dependent values from CCDTAB instead", ctePars.rn_amp); - trlwarn(MsgText); - sprintf(MsgText, "(pctecorr) Readout simulation forward modeling iterations PCTENFOR: %i", - ctePars.n_forward); - trlmessage(MsgText); - sprintf(MsgText, "(pctecorr) Number of iterations used in the parallel transfer PCTENPAR: %i", - ctePars.n_par); - trlmessage(MsgText); - sprintf(MsgText, "(pctecorr) CTE_FRAC: %f", ctePars.scale_frac); - trlmessage(MsgText); - - trlmessage("(pctecorr) PCTETAB read"); + Write the parallel HISTORY information here. + */ + if ((status = populateImageFileWithCTEKeywordValues(&x[0], &cteParallelPars, "parallel"))) + { + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status); } - for (i = 0; i < acs_info->nimsets; i++) - { - clock_t begin = (double)clock(); - if (acs_info->cteAlgorithmGen == 1)//make explicit as not using bool - { - if (forwardModelOnly) + sprintf(MsgText, "(pctecorr) The parallel CTE processing parameters have been read."); + trlmessage(MsgText); + sprintf(MsgText, "(pctecorr) IGNORING read noise level PCTERNOI from PCTETAB: %f. Using amp dependent values from CCDTAB instead", cteParallelPars.rn_amp); + trlwarn(MsgText); + sprintf(MsgText, "(pctecorr) Readout simulation forward modeling iterations PCTENFOR: %i", + cteParallelPars.n_forward); + trlmessage(MsgText); + sprintf(MsgText, "(pctecorr) Number of iterations used in the parallel transfer PCTENPAR: %i", + cteParallelPars.n_par); + trlmessage(MsgText); + sprintf(MsgText, "(pctecorr) CTE_FRAC: %f", cteParallelPars.scale_frac); + trlmessage(MsgText); + + sprintf(MsgText,"\n(pctecorr) NOTE: No serial CTE correction is done for any subarray data.\n"); + trlmessage(MsgText); + /* End read of the parallel CTE parameters */ + + /* + Loop over the imsets as the CTE is applied per amp + */ + char ccdamp[strlen(AMPSTR1)+1]; // string to hold amps on current chip + addPtr(&ptrReg, &ctePars, &freeCTEParamsFast); + { unsigned i; + for (i = 0; i < acs_info->nimsets; i++) { + + // Determine which amps are on the current chip of the input image + ccdamp[0] = '\0'; // "reset" the string for reuse + + // Amps on this chip + parseWFCamps(acs[i].ccdamp, acs[i].chip, ccdamp); + numamps = strlen(ccdamp); + + char corrType[20]; + corrType[0] = '\0'; + for (unsigned nthAmp = 0; nthAmp < numamps; ++nthAmp) { - trlerror("--forwardModelOnly NOT compatible with 1st generation CTE algorithm"); - freeOnExit(&ptrReg); - return (INVALID_VALUE); - } - if ((status = doPCTEGen1(&acs[i], &x[i]))) - return status; - } - else - { - if ((status = doPCTEGen2(&acs[i], &ctePars, &x[i], forwardModelOnly))) - return status; + /* Get the amp letter and number where A:0, B:1, etc. as defined in acs.h */ + amploc = strchr(AMPSORDER, ccdamp[nthAmp]); + ampID = amploc - AMPSORDER; + + /* + Get the amp letter and number where C:0, D:1, etc. as defined in this file. + This order is based upon the arrangement of the chips in the input science + image (chip 2 with amps C and D in the first imset, and then chip 1 with + amps A and B in the second imset. + */ + amplocInCalib = strchr(AMPCALIBORDER, ccdamp[nthAmp]); // This is a full string. + ampIDInCalib = amplocInCalib - AMPCALIBORDER; // This is a number. + + /* + Only perform the serial CTE correction for full-frame, post-SM4 data + */ + if ((acs_info->expstart >= SM4MJD) && (!acs[i].subarray)) { + + startOfSetInCalib = SET_TO_PROCESS[ampIDInCalib]; + strcpy(corrType, "serial"); + sprintf(MsgText,"(pctecorr) Collecting data for Correction Type: %s \n", corrType); + trlmessage(MsgText); + + /* + Loop to control the application of the CTE corrections - the serial + (Extensions 5-8 [A], 9-12 [B], 13-16 [C], 17-20 [D]) correction for each + Amp will be done first, then the parallel (Extensions 1-4) correction will + be done. The "Extensions" refer to the set of extensions in the calibration + reference file pertaining to a particular amp correction. + + As noted at the top of this file, the typical order of processing is + Chip 2 (Amps C and D) and then Chip 1 (Amps A and B). + */ + + if (!skipLoadPrimary) { + initCTEParamsFast(&ctePars, TRAPS, 0, 0, nScaleTableColumns, acs_info->nThreads); + + ctePars.refAndIamgeBinsIdenticle = True; + ctePars.verbose = acs->verbose = 0 ? False : True; + if ((status = allocateCTEParamsFast(&ctePars))) + { + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status); + } + } + + if ((status = loadPCTETAB(cteTabFilename, &ctePars, startOfSetInCalib, skipLoadPrimary))) + { + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status); + } + skipLoadPrimary = True; + + if ((status = getCTEParsFromImageHeader(&x[0], &ctePars))) + { + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status); + } + + ctePars.scale_frac = (acs->expstart - ctePars.cte_date0) / (ctePars.cte_date1 - ctePars.cte_date0); + + /* + Write the amp-dependent serial HISTORY information here. + */ + char amp_corrType[20] = "serial - Amp "; + strcat(amp_corrType, &ccdamp[nthAmp]); + amp_corrType[14] = '\0'; + if ((status = populateImageFileWithCTEKeywordValues(&x[0], &ctePars, amp_corrType))) + { + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return (status); + } + + sprintf(MsgText, "(pctecorr) IGNORING read noise level PCTERNOI from PCTETAB: %f. Using amp dependent values from CCDTAB instead", ctePars.rn_amp); + trlwarn(MsgText); + sprintf(MsgText, "(pctecorr) Readout simulation forward modeling iterations PCTENFOR: %i", + ctePars.n_forward); + trlmessage(MsgText); + sprintf(MsgText, "(pctecorr) Number of iterations used in the parallel transfer PCTENPAR: %i", + ctePars.n_par); + trlmessage(MsgText); + sprintf(MsgText, "(pctecorr) CTE_FRAC: %f", ctePars.scale_frac); + trlmessage(MsgText); + + sprintf(MsgText, "(pctecorr) The %s CTE processing parameters have been read.", corrType); + trlmessage(MsgText); + } // End if block for collecting and reporting serial CTE correction + + /* + The number of acs_info->nimsets represents the number of chips to process for the + input file, and each chip can contain two amps. Since the CTE correction is done on + an amp basis, only the serial correction information for that amp is available, + in addition to the parallel correction information which is the same for all amps. + */ + + clock_t begin = (double)clock(); + + /* Perform the serial CTE correction for only full-frame, post-SM4 data */ + if ((acs_info->expstart >= SM4MJD) && (!acs[i].subarray)) { + /* Serial correction */ + strcpy(corrType, "serial"); + if ((status = doPCTEGen3(&acs[i], &ctePars, &x[i], forwardModelOnly, corrType, &ccdamp, nthAmp, amploc, ampID))) + { + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return status; + } + } // End short block for applying serial CTE correction + + /* Parallel correction */ + strcpy(corrType, "parallel"); + if ((status = doPCTEGen3(&acs[i], &cteParallelPars, &x[i], forwardModelOnly, corrType, &ccdamp, nthAmp, amploc, ampID))) + { + freeOnExit(&ptrReg); + freeOnExit(&ptrParallelReg); + return status; + } + + double time_spent = ((double) clock()- begin +0.0) / CLOCKS_PER_SEC; + sprintf(MsgText,"(pctecorr) CTE run time for current chip: %.2f(s) with %i procs/threads\n", time_spent/acs_info->nThreads, acs_info->nThreads); + trlmessage(MsgText); + } // End loop to accommodate both the serial and parallel CTE corrections } - double time_spent = ((double) clock()- begin +0.0) / CLOCKS_PER_SEC; - sprintf(MsgText,"(pctecorr) CTE run time for current chip: %.2f(s) with %i procs/threads\n", time_spent/acs_info->nThreads, acs_info->nThreads); - trlmessage(MsgText); } - freeOnExit(&ptrReg); PrSwitch("pctecorr", COMPLETE); - - if (acs_info->printtime) - TimeStamp ("PCTECORR complete", acs->rootname); } + if (acs_info->printtime) + TimeStamp ("PCTECORR complete", acs->rootname); + if (!OmitStep(acs_info->pctecorr)) { if (pcteHistory(&acs[0], x[0].globalhdr)) return (status); - if (acs_info->cteAlgorithmGen == 1) - UCteVer(x[0].globalhdr, ACS_GEN1_CTE_NAME, ACS_GEN1_CTE_VER); - else if (acs_info->cteAlgorithmGen == 2) - UCteVer(x[0].globalhdr, ACS_GEN2_CTE_NAME, ACS_GEN2_CTE_VER); - else - assert(0);//unimplemented + UCteVer(x[0].globalhdr, ACS_GEN3_CTE_NAME, ACS_GEN3_CTE_VER); } else if (acs_info->pctecorr == OMIT) { /* this is just to make sure that the PCTECORR flag is set to OMIT in the @@ -326,25 +520,28 @@ int DoCTE (ACSInfo *acs_info, const bool forwardModelOnly) { UCalVer(x[0].globalhdr); UFilename(acs_info->output, x[0].globalhdr); - for (i = 0; i < acs_info->nimsets; i++) { - if (acs_info->verbose) { - sprintf(MsgText, "Outputting chip %d", acs[i].chip); - trlmessage(MsgText); - } + { unsigned i; + for (i = 0; i < acs_info->nimsets; i++) { + if (acs_info->verbose) { + sprintf(MsgText, "Outputting chip %d", acs[i].chip); + trlmessage(MsgText); + } - putSingleGroup(acs_info->output, i+1, &x[i], option); + putSingleGroup(acs_info->output, i+1, &x[i], option); - if (hstio_err()) { - sprintf(MsgText, "Couldn't write imset %d.", i+1); - trlerror(MsgText); - return (status = WRITE_FAILED); - } + if (hstio_err()) { + sprintf(MsgText, "Couldn't write imset %d.", i+1); + trlerror(MsgText); + return (status = WRITE_FAILED); + } + + if (acs_info->printtime) + TimeStamp ("Output written to disk", acs_info->rootname); - if (acs_info->printtime) - TimeStamp ("Output written to disk", acs_info->rootname); + /* Free x, which still has memory allocated. */ + freeSingleGroup (&x[i]); - /* Free x, which still has memory allocated. */ - freeSingleGroup (&x[i]); + } } free(x); diff --git a/pkg/acs/lib/acscte/dopcte-gen2.c b/pkg/acs/lib/acscte/dopcte-gen2.c deleted file mode 100644 index fc906759d..000000000 --- a/pkg/acs/lib/acscte/dopcte-gen2.c +++ /dev/null @@ -1,498 +0,0 @@ -#include -#include -#include -#include -#include - -#include "hstcal_memory.h" -#include "hstcal.h" -#include "hstio.h" - -#include "acs.h" -#include "acsinfo.h" -#include "hstcalerr.h" -#include "trlbuf.h" - -#include "pcte.h" - -# ifdef _OPENMP -# include -# endif -# include "../../../../ctegen2/ctegen2.h" -#include - -int get_amp_array_size_acs_cte(const ACSInfo *acs, SingleGroup *amp, - const int ampID, char *amploc, char *ccdamp, - int *xsize, int *ysize, int *xbeg, - int *xend, int *ybeg, int *yend); - -static int extractAmp(SingleGroup * amp, const SingleGroup * image, const unsigned ampID, CTEParamsFast * ctePars); -static int insertAmp(SingleGroup * amp, const SingleGroup * image, const unsigned ampID, CTEParamsFast * ctePars); -static int alignAmpData(FloatTwoDArray * amp, const unsigned ampID); -static int alignAmp(SingleGroup * amp, const unsigned ampID); - -int doPCTEGen2 (ACSInfo *acs, CTEParamsFast * ctePars, SingleGroup * chipImage, const bool forwardModelOnly) -{ - - /* arguments: - acs i: calibration switches, etc - x io: image to be calibrated; written to in-place - */ - - extern int status; - - if (!acs || !ctePars || !chipImage) - return (status = ALLOCATION_PROBLEM); - - PtrRegister ptrReg; - initPtrRegister(&ptrReg); -#ifdef _OPENMP - if (acs->nThreads > 1) - trlmessage("(pctecorr) Using parallel processing provided by OpenMP inside CTE routine"); -#endif - - char ccdamp[strlen(AMPSTR1)+1]; /* string to hold amps on current chip */ - int numamps; /* number of amps on chip */ - int ampID; /* index amp A:0, B:1, etc. */ - char * amploc; /* pointer to amp character in AMPSORDER */ - int nRows, nColumns; /* int for C array dimension sizes */ - int amp_xsize, amp_ysize; /* int for amp array size (x/y in CCD coords) */ - int amp_xbeg, amp_xend; /* int for beg and end of amp arrays on chip */ - int amp_ybeg, amp_yend; - - /* functions from calacs/lib */ - void parseWFCamps (char *acsamps, int chip, char *ccdamp); - void TimeStamp (char *message, char *rootname); - int PutKeyDbl (Hdr *hd, char *keyword, double value, char *comment); - - /* test whether this we've been given ACS WFC data, since the CTE algorithm - is currently only valid for WFC data. */ - if (acs->detector != WFC_CCD_DETECTOR) { - trlerror("(pctecorr) only valid for WFC CCD data, PCTECORR should be OMIT for all others."); - return (status = ERROR_RETURN); - } - - if (acs->printtime) { - TimeStamp("Starting CTE correction...",""); - } - - /* need to figure out which amps are on this chip */ - ccdamp[0] = '\0'; /* "reset" the string for reuse */ - parseWFCamps(acs->ccdamp, acs->chip, ccdamp); - - /* loop over amps on this chip and do CTE correction */ - numamps = strlen(ccdamp); - - //Now loop over each amp and compute cte correction - {unsigned nthAmp; - for (nthAmp = 0; nthAmp < numamps; ++nthAmp) - { - sprintf(MsgText, "(pctecorr) Performing CTE correction for amp %c", ccdamp[nthAmp]); - trlmessage(MsgText); - - /* get the amp letter and number where A:0, B:1, etc. */ - amploc = strchr(AMPSORDER, ccdamp[nthAmp]); - ampID = *amploc - AMPSORDER[0]; - - ctePars->rn_amp = acs->readnoise[ampID]; - sprintf(MsgText, "(pctecorr) Read noise level from CCDTAB: %f.", ctePars->rn_amp); - trlmessage(MsgText); - - /* get amp array size */ - if (get_amp_array_size_acs_cte(acs, chipImage, ampID, amploc, ccdamp, - &_xsize, &_ysize, &_xbeg, - &_xend, &_ybeg, &_yend)) - { - freeOnExit(&ptrReg); - return (status); - } - - nRows = amp_ysize; - nColumns = amp_xsize; - ctePars->nRows = nRows; - ctePars->nColumns = nColumns; - ctePars->columnOffset = amp_xbeg; - ctePars->rowOffset = amp_ybeg; - // razColumnOffset is used to align the image with the column scaling (SCLBYCOL) in the PCTETAB. - // The SCLBYCOL array is 8192 cols wide which is 2048*4 and therefore does NOT include the physical - // overscan columns (24 for ACS WFC) like this array does for calwf3. - // The raz format is all 4 amps side by side in CDAB order. - ctePars->razColumnOffset = ctePars->columnOffset; - if (ampID == AMP_A || ampID == AMP_B) - ctePars->razColumnOffset += ACS_WFC_N_COLUMNS_PER_CHIP_EXCL_OVERSCAN; - - //This is used for the final output - SingleGroup ampImage; - initSingleGroup(&Image); - addPtr(&ptrReg, &Image, &freeSingleGroup); - if (allocSingleGroupExts(&Image, nColumns, nRows, SCIEXT | ERREXT, False) != 0) - { - freeOnExit(&ptrReg); - return (status = OUT_OF_MEMORY); - } - - // read data from the SingleGroup into an array containing data from - // just one amp - if ((status = extractAmp(&Image, chipImage, ampID, ctePars))) - { - freeOnExit(&ptrReg); - return (status); - } - if ((status = alignAmp(&Image, ampID))) - { - freeOnExit(&ptrReg); - return (status); - } - - //copy to column major storage - SingleGroup columnMajorImage; - initSingleGroup(&columnMajorImage); - addPtr(&ptrReg, &columnMajorImage, &freeSingleGroup); - if (allocSingleGroupExts(&columnMajorImage, nColumns, nRows, SCIEXT, False) != 0) - { - freeOnExit(&ptrReg); - return (status = OUT_OF_MEMORY); - } - if ((status = copySingleGroup(&columnMajorImage, &Image, COLUMNMAJOR))) - { - freeOnExit(&ptrReg); - return (status = ALLOCATION_PROBLEM); - } - - if (forwardModelOnly) - { - // Not needed as we won't be smoothing. - } - else - { - //CALCULATE THE SMOOTH READNOISE IMAGE - trlmessage("(pctecorr) Calculating smooth readnoise image..."); - if (ctePars->noise_mit != 0) - { - trlerror("(pctecorr) Only noise model 0 implemented!"); - freeOnExit(&ptrReg); - return (status = ERROR_RETURN); - } - } - SingleGroup smoothedImage; - initSingleGroup(&smoothedImage); - addPtr(&ptrReg, &smoothedImage, &freeSingleGroup); - if (allocSingleGroupExts(&smoothedImage, nColumns, nRows, SCIEXT, False) != 0) - { - freeOnExit(&ptrReg); - return (status = OUT_OF_MEMORY); - } - setStorageOrder(&smoothedImage, COLUMNMAJOR); - - if (forwardModelOnly) - { - // Don't actually smooth but copy to smoothedImage - if ((status = copySingleGroup(&smoothedImage, &columnMajorImage, COLUMNMAJOR))) - { - freeOnExit(&ptrReg); - return (status = ALLOCATION_PROBLEM); - } - } - else - { - // do some smoothing on the data so we don't amplify the read noise. - if ((status = cteSmoothImage(&columnMajorImage, &smoothedImage, ctePars, ctePars->rn_amp))) - { - freeOnExit(&ptrReg); - return (status); - } - } - trlmessage("(pctecorr) ...complete."); - - trlmessage("(pctecorr) Creating charge trap image..."); - SingleGroup trapPixelMap; - initSingleGroup(&trapPixelMap); - addPtr(&ptrReg, &trapPixelMap, &freeSingleGroup); - if (allocSingleGroupExts(&trapPixelMap, nColumns, nRows, SCIEXT, False) != 0) - { - freeOnExit(&ptrReg); - return (status = OUT_OF_MEMORY); - } - setStorageOrder(&trapPixelMap, COLUMNMAJOR); - if ((status = populateTrapPixelMap(&trapPixelMap, ctePars))) - { - freeOnExit(&ptrReg); - return status; - } - trlmessage("(pctecorr) ...complete."); - - SingleGroup * cteCorrectedImage = &columnMajorImage; - if (forwardModelOnly) - { - trlmessage("(pctecorr) Running forward model simulation..."); - //perform CTE correction - if ((status = forwardModel(&smoothedImage, cteCorrectedImage, &trapPixelMap, ctePars))) - { - freeOnExit(&ptrReg); - return status; - } - } - else - { - trlmessage("(pctecorr) Running correction algorithm..."); - //perform CTE correction - if ((status = inverseCTEBlur(&smoothedImage, cteCorrectedImage, &trapPixelMap, ctePars))) - { - freeOnExit(&ptrReg); - return status; - } - } - freePtr(&ptrReg, &trapPixelMap); - trlmessage("(pctecorr) ...complete."); - - // add 10% correction to error in quadrature. - float totalCounts = 0; - float totalRawCounts = 0; -#ifdef _OPENMP - #pragma omp parallel shared(nRows, nColumns, cteCorrectedImage, smoothedImage, ampImage) -#endif - { - double temp_err; - float threadCounts = 0; - float threadRawCounts = 0; - float delta; - {unsigned k; -#ifdef _OPENMP - #pragma omp for schedule(static) -#endif - //This loop order is row major as more row major storage is accessed than column. - //MORE: look into worth splitting out ops, prob needs a order swap (copy) so perhaps not worth it. - for (k = 0; k < nRows; ++k) - { - {unsigned m; - for (m = 0; m < nColumns; ++m) - { - delta = (PixColumnMajor(cteCorrectedImage->sci.data,k,m) - PixColumnMajor(smoothedImage.sci.data,k,m)); - threadCounts += delta; - threadRawCounts += Pix(ampImage.sci.data, m, k); - temp_err = 0.1 * fabs(delta); - Pix(ampImage.sci.data, m, k) += delta; - - float err2 = Pix(ampImage.err.data, m, k); - err2 *= err2; - Pix(ampImage.err.data, m, k) = sqrt(err2 + temp_err*temp_err); - }} - }} - -#ifdef _OPENMP - #pragma omp critical(deltaAggregate) -#endif - { - totalCounts += threadCounts; - totalRawCounts += threadRawCounts; - } - }//close parallel block - sprintf(MsgText, "(pctecorr) Total count difference (corrected-raw) incurred from correction: %f (%f%%)", totalCounts, totalCounts/totalRawCounts*100); - trlmessage(MsgText); - - // put the CTE corrected data back into the SingleGroup structure - if ((status = alignAmp(&Image, ampID))) - { - freeOnExit(&ptrReg); - return (status); - } - - if ((status = insertAmp(chipImage, &Image, ampID, ctePars))) - { - freeOnExit(&ptrReg); - return (status); - } - - // free mem used by our amp arrays - freePtr(&ptrReg, &Image); - freePtr(&ptrReg, &columnMajorImage); - freePtr(&ptrReg, &smoothedImage); - }} - if (acs->printtime) - TimeStamp("CTE corrections complete",""); - - freeOnExit(&ptrReg); - return (status); -} - -static int extractAmp(SingleGroup * amp, const SingleGroup * image, const unsigned ampID, CTEParamsFast * ctePars) -{ - extern int status; - - if (!amp || !amp->sci.data.data || !image || !image->sci.data.data) - return (status = ALLOCATION_PROBLEM); - - //WARNING - assumes row major storage - assert(amp->sci.data.storageOrder == ROWMAJOR && image->sci.data.storageOrder == ROWMAJOR); - - if (ampID != AMP_A && ampID != AMP_B && ampID != AMP_C && ampID != AMP_D) - { - trlerror("Amp number not recognized, must be 0-3."); - return (status = ERROR_RETURN); - } - - const unsigned nRows = amp->sci.data.ny; - const unsigned nColumns = amp->sci.data.nx; - - unsigned rowSkipLength = image->sci.data.nx; - unsigned offset = ctePars->columnOffset; - - copyOffsetSingleGroup(amp, image, nRows, nColumns, 0, offset, nColumns, rowSkipLength); - return status; -} - -static int insertAmp(SingleGroup * image, const SingleGroup * amp, const unsigned ampID, CTEParamsFast * ctePars) -{ - extern int status; - - if (!amp || !amp->sci.data.data || !image || !image->sci.data.data) - return (status = ALLOCATION_PROBLEM); - - //WARNING - assumes row major storage - assert(amp->sci.data.storageOrder == ROWMAJOR && image->sci.data.storageOrder == ROWMAJOR); - - if (ampID != AMP_A && ampID != AMP_B && ampID != AMP_C && ampID != AMP_D) - { - trlerror("Amp number not recognized, must be 0-3."); - return (status = ERROR_RETURN); - } - - const unsigned nRows = amp->sci.data.ny; - const unsigned nColumns = amp->sci.data.nx; - - unsigned rowSkipLength = image->sci.data.nx; - unsigned offset = ctePars->columnOffset; - - copyOffsetSingleGroup(image, amp, nRows, nColumns, offset, 0, rowSkipLength, nColumns); - return status; -} - -static int alignAmp(SingleGroup * amp, const unsigned ampID) -{ - extern int status; - if (!amp) - return (status = ALLOCATION_PROBLEM); - - //sci data - if (amp->sci.data.data) - { - if ((status = alignAmpData(&->sci.data, ampID))) - return status; - } - //err data - if (amp->err.data.data) - { - if ((status = alignAmpData(&->err.data, ampID))) - return status; - } - //dq data - if (amp->dq.data.data) - assert(0);//unimplemented - - return status; -} - -static int alignAmpData(FloatTwoDArray * amp, const unsigned ampID) -{ - //NOTE: There is a similar version of this in wfc3 - code changes should be reflected in both. - - //Align image quadrants such that the amps are at the bottom left, i.e. aligned with amp C. - extern int status; - - if (!amp || !amp->data) - return (status = ALLOCATION_PROBLEM); - - //Amp C is already correctly aligned - if (ampID == AMP_C) - return status; - - PtrRegister ptrReg; - initPtrRegister(&ptrReg); - - //WARNING - assumes row major storage - assert(amp->storageOrder == ROWMAJOR); - - const unsigned nColumns = amp->nx; - const unsigned nRows = amp->ny; - - //Flip about y axis, i.e. about central column - if (ampID == AMP_B || ampID == AMP_D) - { - //grab a row, flip it, put it back - const unsigned rowLength = nColumns; - {unsigned i; -#ifdef _OPENMP - #pragma omp parallel for shared(amp) private(i) schedule(static) -#endif - for (i = 0; i < nRows; ++i) - { - //find row - float * row = amp->data + i*nColumns; - //flip right to left - float tempPixel; - {unsigned j; - for (j = 0; j < rowLength/2; ++j) - { - tempPixel = row[j]; - row[j] = row[rowLength-1-j]; - row[rowLength-1-j] = tempPixel; - }} - }} - } - - //Only thing left is to flip AB chip upside down - //Flip about x axis, i.e. central row - if (ampID == AMP_A || ampID == AMP_B) - { - //either physically align all or propagate throughout a mechanism to work on the array upside down (flip b quad though) - //we'll just flip all for now. See if there's info in the header specifying amp location rel to pix in file, - //i.e. a way to know whether the chip is 'upside down'. Could then reverse cte corr trail walk direction - Bool allocationFail = False; -#ifdef _OPENMP - #pragma omp parallel shared(amp, allocationFail) -#endif - { - float * tempRow = NULL; - size_t rowSize = nColumns*sizeof(*tempRow); - tempRow = malloc(rowSize); - if (!tempRow) - { -#ifdef _OPENMP - #pragma omp critical(allocationCheck) // This really isn't needed -#endif - { - allocationFail = True; - } - } - -#ifdef _OPENMP - #pragma omp barrier -#endif - if (!allocationFail) - { - {unsigned i; -#ifdef _OPENMP - #pragma omp for schedule(static) -#endif - for (i = 0; i < nRows/2; ++i) - { - float * topRow = amp->data + i*nColumns; - float * bottomRow = amp->data + (nRows-1-i)*nColumns; - memcpy(tempRow, topRow, rowSize); - memcpy(topRow, bottomRow, rowSize); - memcpy(bottomRow, tempRow, rowSize); - }} - } - if (tempRow) - free(tempRow); - }//close parallel block - if (allocationFail) - { - sprintf(MsgText, "Out of memory for 'tempRow' in 'alignAmpData'"); - trlerror(MsgText); - return (status = OUT_OF_MEMORY); - } - } - - return status; -} diff --git a/pkg/acs/lib/acscte/dopcte-gen3.c b/pkg/acs/lib/acscte/dopcte-gen3.c new file mode 100644 index 000000000..3cb7d3c76 --- /dev/null +++ b/pkg/acs/lib/acscte/dopcte-gen3.c @@ -0,0 +1,649 @@ +#include +#include +#include +#include +#include +#include + +#include "hstcal_memory.h" +#include "hstcal.h" +#include "hstio.h" + +#include "acs.h" +#include "acsinfo.h" +#include "hstcalerr.h" +#include "trlbuf.h" + +#include "pcte.h" + +# ifdef _OPENMP +# include +# endif +# include "../../../../ctegen2/ctegen2.h" +#include + +int get_amp_array_size_acs_cte(const ACSInfo *acs, SingleGroup *amp, + char *amploc, char *ccdamp, + int *xsize, int *ysize, int *xbeg, + int *xend, int *ybeg, int *yend); + +static int extractAmp(SingleGroup * amp, const SingleGroup * image, const unsigned ampID, CTEParamsFast * ctePars); +static int insertAmp(SingleGroup * amp, const SingleGroup * image, const unsigned ampID, CTEParamsFast * ctePars); +static int alignAmpData(FloatTwoDArray * amp, const unsigned ampID); +static int alignAmp(SingleGroup * amp, const unsigned ampID); +static int rotateAmp(SingleGroup * amp, const unsigned ampID, bool derotate, char ccdamp); +static int rotateAmpData(FloatTwoDArray * amp, const unsigned ampID); +static int derotateAmpData(FloatTwoDArray * amp, const unsigned ampID); + +static void transpose(FloatTwoDArray *amp); +static void side2sideFlip(FloatTwoDArray *amp); +static void top2bottomFlip(FloatTwoDArray *amp); + +int doPCTEGen3 (ACSInfo *acs, CTEParamsFast * ctePars, SingleGroup * chipImage, const bool forwardModelOnly, char * corrType, char *ccdamp, int nthAmp, char *amploc, int ampID) + +{ + + /* arguments: + acs i: calibration switches, etc + x io: image to be calibrated; written to in-place + */ + + extern int status; + + if (!acs || !ctePars || !chipImage) + return (status = ALLOCATION_PROBLEM); + + PtrRegister ptrReg; + initPtrRegister(&ptrReg); +#ifdef _OPENMP + if (acs->nThreads > 1) + trlmessage("(pctecorr) Using multiprocessing provided by OpenMP inside CTE routine."); +#endif + + int nRows, nColumns; /* int for C array dimension sizes */ + int amp_xsize, amp_ysize; /* int for amp array size (x/y in CCD coords) */ + int amp_xbeg, amp_xend; /* int for beg and end of amp arrays on chip */ + int amp_ybeg, amp_yend; + + /* functions from calacs/lib */ + void TimeStamp (char *message, char *rootname); + int PutKeyDbl (Hdr *hd, char *keyword, double value, char *comment); + + if (acs->printtime) { + TimeStamp("Starting CTE correction...",""); + } + + sprintf(MsgText, "(pctecorr) Performing CTE correction type %s for amp %c.", corrType, ccdamp[nthAmp]); + trlmessage(MsgText); + + ctePars->rn_amp = acs->readnoise[ampID]; + sprintf(MsgText, "(pctecorr) Read noise level from CCDTAB: %f.", ctePars->rn_amp); + trlmessage(MsgText); + + /* get amp array size */ + if (get_amp_array_size_acs_cte(acs, chipImage, amploc, ccdamp, + &_xsize, &_ysize, &_xbeg, + &_xend, &_ybeg, &_yend)) + { + freeOnExit(&ptrReg); + return (status); + } + + nRows = amp_ysize; + nColumns = amp_xsize; + ctePars->nRows = nRows; + ctePars->nColumns = nColumns; + ctePars->columnOffset = amp_xbeg; + ctePars->rowOffset = amp_ybeg; + // razColumnOffset is used to align the image with the column scaling (SCLBYCOL) in the PCTETAB. + // The SCLBYCOL array is 8192 cols wide which is 2048*4 and therefore does NOT include the physical + // overscan columns (24 for ACS WFC) like this array does for calwf3. + // The raz format is all 4 amps side by side in CDAB order. + ctePars->razColumnOffset = ctePars->columnOffset; + if (ampID == AMP_A || ampID == AMP_B) + ctePars->razColumnOffset += ACS_WFC_N_COLUMNS_PER_CHIP_EXCL_OVERSCAN; + + //This is used for the final output + SingleGroup ampImage; + initSingleGroup(&Image); + addPtr(&ptrReg, &Image, &freeSingleGroup); + if (allocSingleGroupExts(&Image, nColumns, nRows, SCIEXT | ERREXT, False) != 0) + { + freeOnExit(&ptrReg); + return (status = OUT_OF_MEMORY); + } + + // read data from the SingleGroup into an array containing data from + // just one amp + if ((status = extractAmp(&Image, chipImage, ampID, ctePars))) + { + freeOnExit(&ptrReg); + return (status); + } + + // rotate amp in preparation for the serial CTE correction + // This step is done ONLY to accommodate the serial CTE correction + // which preceeds the standard parallel CTE correction. + if (strcmp(corrType, "serial") == 0) + { + sprintf(MsgText, "\n(pctecorr) Invoking rotation for CTE Correction Type: %s", corrType); + trlmessage(MsgText); + if ((status = rotateAmp(&Image, ampID, False, ccdamp[nthAmp]))) + { + freeOnExit(&ptrReg); + return (status); + } + } + + if ((status = alignAmp(&Image, ampID))) + { + freeOnExit(&ptrReg); + return (status); + } + + //copy to column major storage + SingleGroup columnMajorImage; + initSingleGroup(&columnMajorImage); + addPtr(&ptrReg, &columnMajorImage, &freeSingleGroup); + if (allocSingleGroupExts(&columnMajorImage, nColumns, nRows, SCIEXT, False) != 0) + { + freeOnExit(&ptrReg); + return (status = OUT_OF_MEMORY); + } + if ((status = copySingleGroup(&columnMajorImage, &Image, COLUMNMAJOR))) + { + freeOnExit(&ptrReg); + return (status = ALLOCATION_PROBLEM); + } + + if (forwardModelOnly) + { + // Not needed as we won't be smoothing. + } + else + { + //CALCULATE THE SMOOTH READNOISE IMAGE + trlmessage("(pctecorr) Calculating smooth readnoise image..."); + if (ctePars->noise_mit != 0) + { + trlerror("(pctecorr) Only noise model 0 implemented!"); + freeOnExit(&ptrReg); + return (status = ERROR_RETURN); + } + } + SingleGroup smoothedImage; + initSingleGroup(&smoothedImage); + addPtr(&ptrReg, &smoothedImage, &freeSingleGroup); + if (allocSingleGroupExts(&smoothedImage, nColumns, nRows, SCIEXT, False) != 0) + { + freeOnExit(&ptrReg); + return (status = OUT_OF_MEMORY); + } + setStorageOrder(&smoothedImage, COLUMNMAJOR); + + if (forwardModelOnly) + { + // Don't actually smooth but copy to smoothedImage + if ((status = copySingleGroup(&smoothedImage, &columnMajorImage, COLUMNMAJOR))) + { + freeOnExit(&ptrReg); + return (status = ALLOCATION_PROBLEM); + } + } + else + { + // do some smoothing on the data so we don't amplify the read noise. + if ((status = cteSmoothImage(&columnMajorImage, &smoothedImage, ctePars, ctePars->rn_amp))) + { + freeOnExit(&ptrReg); + return (status); + } + } + trlmessage("(pctecorr) ...complete."); + + trlmessage("(pctecorr) Creating charge trap image..."); + SingleGroup trapPixelMap; + initSingleGroup(&trapPixelMap); + addPtr(&ptrReg, &trapPixelMap, &freeSingleGroup); + if (allocSingleGroupExts(&trapPixelMap, nColumns, nRows, SCIEXT, False) != 0) + { + freeOnExit(&ptrReg); + return (status = OUT_OF_MEMORY); + } + setStorageOrder(&trapPixelMap, COLUMNMAJOR); + if ((status = populateTrapPixelMap(&trapPixelMap, ctePars))) + { + freeOnExit(&ptrReg); + return status; + } + trlmessage("(pctecorr) ...complete."); + + SingleGroup * cteCorrectedImage = &columnMajorImage; + if (forwardModelOnly) + { + trlmessage("(pctecorr) Running forward model simulation..."); + //perform CTE correction + if ((status = forwardModel(&smoothedImage, cteCorrectedImage, &trapPixelMap, ctePars))) + { + freeOnExit(&ptrReg); + return status; + } + } + else + { + trlmessage("(pctecorr) Running correction algorithm..."); + //perform CTE correction + if ((status = inverseCTEBlur(&smoothedImage, cteCorrectedImage, &trapPixelMap, ctePars))) + { + freeOnExit(&ptrReg); + return status; + } + } + freePtr(&ptrReg, &trapPixelMap); + trlmessage("(pctecorr) ...complete."); + + // add 10% correction to error in quadrature. + float totalCounts = 0; + float totalRawCounts = 0; +#ifdef _OPENMP + #pragma omp parallel shared(nRows, nColumns, cteCorrectedImage, smoothedImage, ampImage) +#endif + { + double temp_err; + float threadCounts = 0; + float threadRawCounts = 0; + float delta; + { unsigned k; +#ifdef _OPENMP + #pragma omp for schedule(static) +#endif + //This loop order is row major as more row major storage is accessed than column. + //MORE: look into worth splitting out ops, prob needs a order swap (copy) so perhaps not worth it. + for (k = 0; k < nRows; ++k) + { + { unsigned m; + for (m = 0; m < nColumns; ++m) + { + delta = (PixColumnMajor(cteCorrectedImage->sci.data,k,m) - PixColumnMajor(smoothedImage.sci.data,k,m)); + threadCounts += delta; + threadRawCounts += Pix(ampImage.sci.data, m, k); + temp_err = 0.1 * fabs(delta); + Pix(ampImage.sci.data, m, k) += delta; + + float err2 = Pix(ampImage.err.data, m, k); + err2 *= err2; + Pix(ampImage.err.data, m, k) = sqrt(err2 + temp_err*temp_err); + } + } + } + } + +#ifdef _OPENMP + #pragma omp critical(deltaAggregate) +#endif + { + totalCounts += threadCounts; + totalRawCounts += threadRawCounts; + } + }//close parallel block + sprintf(MsgText, "(pctecorr) Total count difference (corrected-raw) incurred from correction: %f (%f%%)", totalCounts, totalCounts/totalRawCounts*100); + trlmessage(MsgText); + + // put the CTE corrected data back into the SingleGroup structure + if ((status = alignAmp(&Image, ampID))) + { + freeOnExit(&ptrReg); + return (status); + } + + // Derotate amp so the parallel CTE correction can proceed + if (strcmp(corrType, "serial") == 0) + { + sprintf(MsgText, "\n(pctecorr) Invoking de-rotation for CTE Correction Type: %s", corrType); + trlmessage(MsgText); + if ((status = rotateAmp(&Image, ampID, True, ccdamp[nthAmp]))) + { + freeOnExit(&ptrReg); + return (status); + } + } + + if ((status = insertAmp(chipImage, &Image, ampID, ctePars))) + { + freeOnExit(&ptrReg); + return (status); + } + + // free mem used by our amp arrays + freePtr(&ptrReg, &Image); + freePtr(&ptrReg, &columnMajorImage); + freePtr(&ptrReg, &smoothedImage); + if (acs->printtime) + TimeStamp("CTE corrections complete",""); + + freeOnExit(&ptrReg); + return (status); +} + +static int extractAmp(SingleGroup * amp, const SingleGroup * image, const unsigned ampID, CTEParamsFast * ctePars) +{ + extern int status; + + if (!amp || !amp->sci.data.data || !image || !image->sci.data.data) + return (status = ALLOCATION_PROBLEM); + + //WARNING - assumes row major storage + assert(amp->sci.data.storageOrder == ROWMAJOR && image->sci.data.storageOrder == ROWMAJOR); + + if (ampID != AMP_A && ampID != AMP_B && ampID != AMP_C && ampID != AMP_D) + { + trlerror("Amp number not recognized, must be 0-3."); + return (status = ERROR_RETURN); + } + + const unsigned nRows = amp->sci.data.ny; + const unsigned nColumns = amp->sci.data.nx; + + unsigned rowSkipLength = image->sci.data.nx; + unsigned offset = ctePars->columnOffset; + + copyOffsetSingleGroup(amp, image, nRows, nColumns, 0, offset, nColumns, rowSkipLength); + return status; +} + +/* + The rotation routines support the serial CTE correction such that the data + is configured to mimic the parallel CTE data. The idea is to + rotate the amp to put the serial trails in the same orientation as + those of the parallel trails and then the existing CTE functions are applied. + Essentially, the CTE correction routines can then be applied in the identical + manner for BOTH the serial and the parallel CTE correction cases. + + Once the serial CTE correction has been performed, the amps then need to be + "de-rotated" so the parallel CTE correction can then be applied. +*/ +static int rotateAmp(SingleGroup * amp, const unsigned ampID, bool derotate, char ccdamp) +{ + extern int status; + + if (!amp) + return (status = ALLOCATION_PROBLEM); + + // The case of rotating the amp into position to mimic parallel processing + if (!derotate) + { + //sci data + if (amp->sci.data.data) + { + sprintf(MsgText, "(pctecorr) Rotation for Amp: %c\n", ccdamp); + trlmessage(MsgText); + if ((status = rotateAmpData(&->sci.data, ampID))) + return status; + } + + //err data + if (amp->err.data.data) + { + if ((status = rotateAmpData(&->err.data, ampID))) + return status; + } + } + // The case of returning the amp to its original orientation + else + { + //sci data + if (amp->sci.data.data) + { + sprintf(MsgText, "(pctecorr) Derotation for Amp: %c\n", ccdamp); + trlmessage(MsgText); + if ((status = derotateAmpData(&->sci.data, ampID))) + return status; + } + + //err data + if (amp->err.data.data) + { + if ((status = derotateAmpData(&->err.data, ampID))) + return status; + } + } + + return status; +} + +static int rotateAmpData(FloatTwoDArray * amp, const unsigned ampID) +{ + extern int status; + if (!amp || !amp->data) + return (status = ALLOCATION_PROBLEM); + + /* + Rotate the amp to put the serial trails in the same orientation + as the parallel trails would be. A rotation requires a transpose + and then a flip. Always transpose the data first. + + */ + transpose(amp); + + /* + To complete the correct rotation, the flip either has to be + from side-to-side about the central column or top-to-bottom + about the central row. + */ + if (ampID == AMP_B || ampID == AMP_C) { + side2sideFlip(amp); + } else { + top2bottomFlip(amp); + } + + return status; +} + +static int derotateAmpData(FloatTwoDArray * amp, const unsigned ampID) +{ + extern int status; + if (!amp || !amp->data) + return (status = ALLOCATION_PROBLEM); + + /* + To derotate the amp, you are reversing the direction of the initial + rotation. Always transpose the data first, and then apply the + side-to-side or top-to-bottom flip to put the amp back into its + original orientation so the parallel CTE correction can proceed. + */ + transpose(amp); + + /* + The rotation here is in the opposite direction from the initial + rotation for the amp in question, so if a side-to-side flip were + done initially, now do a top-to-bottom flip (for example). + */ + if (ampID == AMP_B || ampID == AMP_C) { + top2bottomFlip(amp); + } else { + side2sideFlip(amp); + } + + return status; +} + +/* + Transpose the amp data +*/ +static void transpose(FloatTwoDArray * amp) +{ + const unsigned nRows = amp->ny; + const unsigned nColumns = amp->nx; + float temp; + + for (unsigned i = 0; i < nRows; i++) { + for (unsigned j = i; j < nColumns; j++) { + temp = PPix(amp, j, i); + PPix(amp, j, i) = PPix(amp, i, j); + PPix(amp, i, j) = temp; + } + } +} + +/* + Flip the amp about the Y-axis central column (i.e., flip from left to right) +*/ +static void side2sideFlip(FloatTwoDArray * amp) +{ + const unsigned nRows = amp->ny; + const unsigned nColumns = amp->nx; + float temp; + + for (unsigned i = 0; i < nRows; i++) { + for (unsigned j = 0; j < nColumns/2; j++) { + temp = PPix(amp, j, i); + PPix(amp, j, i) = PPix(amp, nColumns - j - 1, i); + PPix(amp, nColumns - j - 1, i) = temp; + } + } +} + +static int insertAmp(SingleGroup * image, const SingleGroup * amp, const unsigned ampID, CTEParamsFast * ctePars) +{ + extern int status; + + if (!amp || !amp->sci.data.data || !image || !image->sci.data.data) + return (status = ALLOCATION_PROBLEM); + + //WARNING - assumes row major storage + assert(amp->sci.data.storageOrder == ROWMAJOR && image->sci.data.storageOrder == ROWMAJOR); + + if (ampID != AMP_A && ampID != AMP_B && ampID != AMP_C && ampID != AMP_D) + { + trlerror("Amp number not recognized, must be 0-3."); + return (status = ERROR_RETURN); + } + + const unsigned nRows = amp->sci.data.ny; + const unsigned nColumns = amp->sci.data.nx; + + unsigned rowSkipLength = image->sci.data.nx; + unsigned offset = ctePars->columnOffset; + + copyOffsetSingleGroup(image, amp, nRows, nColumns, offset, 0, rowSkipLength, nColumns); + return status; +} + +static int alignAmp(SingleGroup * amp, const unsigned ampID) +{ + extern int status; + if (!amp) + return (status = ALLOCATION_PROBLEM); + + //sci data + if (amp->sci.data.data) + { + if ((status = alignAmpData(&->sci.data, ampID))) + return status; + } + //err data + if (amp->err.data.data) + { + if ((status = alignAmpData(&->err.data, ampID))) + return status; + } + //dq data + if (amp->dq.data.data) + assert(0);//unimplemented + + return status; +} + +static int alignAmpData(FloatTwoDArray * amp, const unsigned ampID) +{ + //NOTE: There is a similar version of this in wfc3 - code changes should be reflected in both. + + //Align image quadrants such that the amps are at the bottom left, i.e. aligned with amp C. + extern int status; + + if (!amp || !amp->data) + return (status = ALLOCATION_PROBLEM); + + //Amp C is already correctly aligned + if (ampID == AMP_C) + return status; + + //WARNING - assumes row major storage + assert(amp->storageOrder == ROWMAJOR); + + //Flip about y axis, i.e. about central column + if (ampID == AMP_B || ampID == AMP_D) + { + side2sideFlip(amp); + } + + //Only thing left is to flip AB chip upside down + //Flip about x axis, i.e. central row + if (ampID == AMP_A || ampID == AMP_B) + { + top2bottomFlip(amp); + } + + return status; +} + +/* + Flip the amp about the X-axis central row (i.e., flip from top to bottom) + + Note: This routine was originally in alignAmpData flow of processing. It + is now encapulated here as it is used multiple times due to the rotations + needed to accommodate the serial CTE correction. The original OPENMP + specifications have been left intact. +*/ +static void top2bottomFlip(FloatTwoDArray * amp) +{ + + const unsigned nColumns = amp->nx; + const unsigned nRows = amp->ny; + + Bool allocationFail = False; +#ifdef _OPENMP + #pragma omp parallel shared(amp, allocationFail) +#endif + { + float * tempRow = NULL; + size_t rowSize = nColumns*sizeof(*tempRow); + tempRow = malloc(rowSize); + if (!tempRow) + { +#ifdef _OPENMP + #pragma omp critical(allocationCheck) // This really isn't needed +#endif + { + allocationFail = True; + } + } + +#ifdef _OPENMP + #pragma omp barrier +#endif + if (!allocationFail) + { + { unsigned i; +#ifdef _OPENMP + #pragma omp for schedule(static) +#endif + for (i = 0; i < nRows/2; ++i) + { + float * topRow = amp->data + i*nColumns; + float * bottomRow = amp->data + (nRows-1-i)*nColumns; + memcpy(tempRow, topRow, rowSize); + memcpy(topRow, bottomRow, rowSize); + memcpy(bottomRow, tempRow, rowSize); + } + } + } + if (tempRow) + free(tempRow); + }//close parallel block + if (allocationFail) + { + sprintf(MsgText, "Out of memory for 'tempRow' in 'alignAmpData'"); + trlerror(MsgText); + } +} diff --git a/pkg/acs/lib/acscte/dopcte.c b/pkg/acs/lib/acscte/dopcte.c index ef44641f3..dc2f5977d 100644 --- a/pkg/acs/lib/acscte/dopcte.c +++ b/pkg/acs/lib/acscte/dopcte.c @@ -15,7 +15,7 @@ int get_amp_array_size_acs_cte(const ACSInfo *acs, SingleGroup *x, - const int amp, char *amploc, char *ccdamp, + char *amploc, char *ccdamp, int *xsize, int *ysize, int *xbeg, int *xend, int *ybeg, int *yend); static int make_amp_array(const ACSInfo *acs, const SingleGroup *im, @@ -32,232 +32,6 @@ static int unmake_amp_array(const ACSInfo *acs, const SingleGroup *im, double amp_err_array[arr1*arr2]); -/* Perform a pixel based CTE correction on the SCI data extension of ACS CCD - data. Parameters of the CTE characterization are read from the PCTE reference - file. - - Originally adapted from code written by Jay Anderson and rewritten by - Pey Lian Lim as a standalone application that operated on FLT files. - - MRD 10 Mar 2011 -*/ -int doPCTEGen1 (ACSInfo *acs, SingleGroup *x) { - - /* arguments: - acs i: calibration switches, etc - x io: image to be calibrated; written to in-place - */ - - extern int status; - - /* interpolated cte profile shape parameters */ - double chg_leak[MAX_TAIL_LEN*NUM_LOGQ]; - double chg_open[MAX_TAIL_LEN*NUM_LOGQ]; - - /* interpolated profile charge parameters */ - double dtde_q[MAX_PHI]; - - /* arrays interpolated at each parameterized charge level */ - double chg_leak_lt[MAX_TAIL_LEN*NUM_LEV]; - double chg_open_lt[MAX_TAIL_LEN*NUM_LEV]; - double dpde_l[NUM_LEV]; - - /* structure to hold CTE parameters from file */ - ACSCTEParams pars; - - /* temporary variable used during final error calculation */ - double temp_err; - - /* iteration variable */ - int i, k, m; - - char ccdamp[strlen(AMPSTR1)+1]; /* string to hold amps on current chip */ - int numamps; /* number of amps on chip */ - int amp; /* index amp A:0, B:1, etc. */ - char * amploc; /* pointer to amp character in AMPSORDER */ - int amp_arr1, amp_arr2; /* int for C array dimension sizes */ - int amp_xsize, amp_ysize; /* int for amp array size (x/y in CCD coords) */ - int amp_xbeg, amp_xend; /* int for beg and end of amp arrays on chip */ - int amp_ybeg, amp_yend; - - /* make arrays to hold data amp by amp. - these can be large arrays so it's best to declare them as pointers and - get the space for the ararys using malloc */ - double * amp_sci_arr; /* original sci data */ - double * amp_err_arr; /* original err data */ - double * amp_sig_arr; /* decomposed signal */ - double * amp_nse_arr; /* decomposed readout error */ - double * amp_cor_arr; /* cte corrected data */ - - /* in this algorithm each pixel has it's own CTE scaling, - so we need an array for that. */ - double * cte_frac_arr; - - /* functions from calacs/lib */ - void parseWFCamps (char *acsamps, int chip, char *ccdamp); - void TimeStamp (char *message, char *rootname); - int PutKeyDbl (Hdr *hd, char *keyword, double value, char *comment); - - /* test whether this we've been given ACS WFC data, since the CTE algorithm - is currently only valid for WFC data. */ - if (acs->detector != WFC_CCD_DETECTOR) { - trlerror("(pctecorr) only valid for WFC CCD data, PCTECORR should be OMIT for all others."); - return (status = ERROR_RETURN); - } - - if (acs->printtime) { - TimeStamp("Starting CTE correction...",""); - } - - /**************** read and calculate parameters of CTE model ************/ - if (PixCteParams(acs->pcte.name, acs->expstart, &pars)) { - return (status); - } - - if (CompareCteParams(x, &pars)) { - return (status); - } - - sprintf(MsgText, "(pctecorr) Read noise level PCTERNCL: %f", pars.rn_clip); - trlmessage(MsgText); - sprintf(MsgText, "(pctecorr) Readout simulation iterations PCTESMIT: %i", - pars.sim_nit); - trlmessage(MsgText); - sprintf(MsgText, "(pctecorr) Number of readout shifts PCTESHFT: %i", - pars.shft_nit); - trlmessage(MsgText); - sprintf(MsgText, "(pctecorr) CTE_FRAC: %f", pars.cte_frac); - trlmessage(MsgText); - - /* also add cte_frac as header keyword */ - if ((acs->chip == 2) || (acs->subarray == YES)) { - if (PutKeyDbl(x->globalhdr, "PCTEFRAC", pars.cte_frac, - "CTE scaling factor")) { - trlerror("(pctecorr) Error writing PCTEFRAC to image header"); - return (status = HEADER_PROBLEM); - } - } - - if (InterpolatePsi(pars.chg_leak, pars.psi_node, chg_leak, chg_open)) { - return (status); - } - if (InterpolatePhi(pars.dtde_l, pars.q_dtde, pars.shft_nit, dtde_q)) { - return (status); - } - if (FillLevelArrays(chg_leak, chg_open, dtde_q, pars.levels, chg_leak_lt, - chg_open_lt, dpde_l)) { - return (status); - } - /********* done reading and calculating CTE model parameters ************/ - - /* need to figure out which amps are on this chip */ - ccdamp[0] = '\0'; /* "reset" the string for reuse */ - parseWFCamps(acs->ccdamp, acs->chip, ccdamp); - - /* loop over amps on this chip and do CTE correction */ - numamps = strlen(ccdamp); - for (i = 0; i < numamps; i++) { - sprintf(MsgText, "(pctecorr) Performing CTE correction for amp %c", - ccdamp[i]); - trlmessage(MsgText); - - /* get the amp letter and number where A:0, B:1, etc. */ - amploc = strchr(AMPSORDER, ccdamp[i]); - amp = *amploc - AMPSORDER[0]; - - /* get amp array size */ - if (get_amp_array_size_acs_cte(acs, x, amp, amploc, ccdamp, - &_xsize, &_ysize, &_xbeg, - &_xend, &_ybeg, &_yend)) { - return (status); - } - - amp_arr1 = amp_ysize; - amp_arr2 = amp_xsize; - - /* allocate space to hold this amp's data in its various forms */ - amp_sci_arr = (double *) malloc(amp_arr1 * amp_arr2 * sizeof(double)); - amp_err_arr = (double *) malloc(amp_arr1 * amp_arr2 * sizeof(double)); - amp_sig_arr = (double *) malloc(amp_arr1 * amp_arr2 * sizeof(double)); - amp_nse_arr = (double *) malloc(amp_arr1 * amp_arr2 * sizeof(double)); - amp_cor_arr = (double *) malloc(amp_arr1 * amp_arr2 * sizeof(double)); - cte_frac_arr = (double *) malloc(amp_arr1 * amp_arr2 * sizeof(double)); - - /* read data from the SingleGroup into an array containing data from - just one amp */ - if (make_amp_array(acs, x, amp, amp_arr1, amp_arr2, amp_xbeg, amp_ybeg, - amp_sci_arr, amp_err_arr)) { - return (status); - } - - /* fill cte_frac_arr - - To match cte_frac_arr: v8.1.3 -> v8.2 - k*amp_arr2+m -> (k+20)*(amp_arr2+24) + (m+24) - - offsety (LTV2) always 0 for fullframe and -1 for 2K subarray. - */ - for (k = 0; k < amp_arr1; k++) { - for (m = 0; m < amp_arr2; m++) { - cte_frac_arr[k*amp_arr2 + m] = pars.cte_frac * - pars.col_scale[m*NAMPS + amp] * - ((double) k - acs->offsety) / CTE_REF_ROW; - } - } - - /* do some smoothing on the data so we don't amplify the read noise. - data should be in electrons. */ - if (DecomposeRN(amp_arr1, amp_arr2, amp_sci_arr, - pars.rn_clip, pars.noise_model, - amp_sig_arr, amp_nse_arr)) { - return (status); - } - - /* perform CTE correction */ - int onecpu = acs->nThreads <= 1 ? YES : NO; - if (FixYCte(amp_arr1, amp_arr2, amp_sig_arr, amp_cor_arr, pars.sim_nit, - pars.shft_nit, pars.sub_thresh, cte_frac_arr, pars.levels, - dpde_l, chg_leak_lt, chg_open_lt, onecpu)) { - return (status); - } - - /* add readout noise back and convert corrected data back to DN. - add 10% correction to error in quadrature. */ - for (k = 0; k < amp_arr1; k++) { - for (m = 0; m < amp_arr2; m++) { - amp_cor_arr[k*amp_arr2 + m] = amp_cor_arr[k*amp_arr2 + m] + - amp_nse_arr[k*amp_arr2 + m]; - - temp_err = 0.1 * fabs(amp_cor_arr[k*amp_arr2 + m] - - amp_sci_arr[k*amp_arr2 + m]); - amp_err_arr[k*amp_arr2 + m] = sqrt( - pow(amp_err_arr[k*amp_arr2 + m],2) + pow(temp_err,2)); - } - } - - /* put the CTE corrected data back into the SingleGroup structure */ - if (unmake_amp_array(acs, x, amp, amp_arr1, amp_arr2, amp_xbeg, amp_ybeg, - amp_cor_arr, amp_err_arr)) { - return (status); - } - - /* free space used by our amp arrays */ - free(amp_sci_arr); - free(amp_err_arr); - free(amp_sig_arr); - free(amp_nse_arr); - free(amp_cor_arr); - free(cte_frac_arr); - } - - if (acs->printtime) { - TimeStamp("CTE corrections complete...",""); - } - - return (status); -} - - /* Returns the x/y dimensions for an array that holds data readout through a single amp. currently only works for ACS WFC data. @@ -269,7 +43,7 @@ int doPCTEGen1 (ACSInfo *acs, SingleGroup *x) { - MRD 14 Mar 2011 */ int get_amp_array_size_acs_cte(const ACSInfo *acs, SingleGroup *x, - const int amp, char *amploc, char *ccdamp, + char *amploc, char *ccdamp, int *xsize, int *ysize, int *xbeg, int *xend, int *ybeg, int *yend) { extern int status;