Index: /branches/eam_branches/relphot.20200519/Makefile =================================================================== --- /branches/eam_branches/relphot.20200519/Makefile (revision 41361) +++ /branches/eam_branches/relphot.20200519/Makefile (revision 41362) @@ -32,4 +32,6 @@ $(SRC)/BrightCatalog.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/TGroupOps.$(ARCH).o \ +$(SRC)/ReadTGroupFile.$(ARCH).o \ $(SRC)/SetSignals.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o \ @@ -89,4 +91,6 @@ $(SRC)/BrightCatalog.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/TGroupOps.$(ARCH).o \ +$(SRC)/ReadTGroupFile.$(ARCH).o \ $(SRC)/SetSignals.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o \ Index: /branches/eam_branches/relphot.20200519/include/relphot.h =================================================================== --- /branches/eam_branches/relphot.20200519/include/relphot.h (revision 41361) +++ /branches/eam_branches/relphot.20200519/include/relphot.h (revision 41362) @@ -90,4 +90,29 @@ } Mosaic; +typedef struct { + unsigned int start; + unsigned int stop; + short photcode; + float McalPSF; + float McalAPER; + float dMcal; + float dMsys; + float McalChiSq; + unsigned short nFitPhotom; + unsigned int flags; + + off_t Nimage; + off_t NIMAGE; + off_t *image; + + off_t Nmeasure; + off_t NMEASURE; + off_t *measure; + + off_t Ncatalog; + off_t NCATALOG; + off_t *catalog; +} TGroup; + typedef enum { STATS_NONE, @@ -145,4 +170,5 @@ int Ncal; int Nmos; + int Ngrp; int Ngrid; @@ -305,4 +331,5 @@ int VERBOSE_IMAGE; int MOSAIC_ZEROPT; +int TGROUP_ZEROPT; int TEST_IMAGE1; @@ -403,4 +430,5 @@ void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr)); void clean_mosaics PROTO((void)); +void clean_tgroups PROTO((void)); void clean_stars PROTO((Catalog *catalog, int Ncatalog)); int corner_check PROTO((double *x1, double *y1, double *x2, double *y2)); @@ -410,4 +438,5 @@ void findImages PROTO((Catalog *catalog, int Ncatalog, int doImageList)); int findMosaics PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); +int findTGroups PROTO((Catalog *catalog, int Ncatalog)); void clearImages (void); @@ -428,4 +457,6 @@ void freeImageBins PROTO((int Ncatalog, int doImageList)); void freeMosaicBins PROTO((int Ncatalog, int doMosaicList)); +void freeTGroupBins PROTO((int Ncatalog, int doTGroupList)); + void free_catalogs PROTO((Catalog *catalog, int Ncatalog)); int gcatalog PROTO((Catalog *catalog, int FINAL)); @@ -437,4 +468,5 @@ float getMgrid PROTO((off_t meas, int cat)); float getMmos PROTO((off_t meas, int cat)); +float getMgrp PROTO((off_t meas, int cat)); float getMrel PROTO((Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source)); short getUbercalDist PROTO((off_t meas, int cat)); @@ -452,4 +484,6 @@ void initMosaicGrid PROTO((Image *image, off_t Nimage)); void initMosaics PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage)); +void initTGroups PROTO((Image *subset, off_t *LineNumber, off_t Nsubset)); +void initTGroupBins PROTO((Catalog *catalog, int Ncatalog)); void initMrel PROTO((Catalog *catalog, int Ncatalog)); RelphotMode initialize PROTO((int argc, char **argv)); @@ -460,4 +494,10 @@ int liststats_irls PROTO((StatDataSet *dataset, int Npoints, StatType *stats)); +unsigned int *ReadTGroupFile PROTO((FILE *f, int *nelem)); +void loadTGroups PROTO((char *filename)); +void initTGroupsMcal PROTO((void)); +TGroup *getTGroupForImage PROTO((off_t im)); +off_t findTGroup PROTO((unsigned int start, int photcode)); + Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile)); Catalog *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog, char *syncfile)); @@ -470,4 +510,6 @@ void matchImage PROTO((Catalog *catalog, off_t meas, int cat, int doImageList)); void matchMosaics PROTO((Catalog *catalog, off_t meas, int cat, int doMosaicList)); +void matchTGroups PROTO((Catalog *catalog, off_t meas, int cat)); + double opening_angle PROTO((double x1, double y1, double x2, double y2, double x3, double y3)); void plot_chisq PROTO((Catalog *catalog, int Ncatalog)); @@ -491,8 +533,10 @@ int setExclusions PROTO((Catalog *catalog, int Ncatalog, int verbose)); void setMcal PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); -void setMcalFinal PROTO((void)); +void setMcalFromMosaics PROTO((void)); +void setMcalFromTGroups PROTO((void)); int setMcalOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); void setMgrid PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); int setMmos PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); +int setMgrp PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage)); Index: /branches/eam_branches/relphot.20200519/src/MosaicOps.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/MosaicOps.c (revision 41361) +++ /branches/eam_branches/relphot.20200519/src/MosaicOps.c (revision 41362) @@ -720,5 +720,5 @@ } -void setMcalFinal () { +void setMcalFromMosaics () { off_t i, j, im, Nimage; Index: /branches/eam_branches/relphot.20200519/src/ReadTGroupFile.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/ReadTGroupFile.c (revision 41362) +++ /branches/eam_branches/relphot.20200519/src/ReadTGroupFile.c (revision 41362) @@ -0,0 +1,100 @@ +# include "relphot.h" + +# define NBUFFER 0x10000 +# define IsCSV FALSE +# define COLUMN_NUMBER 1 + +unsigned int *ReadTGroupFile (FILE *f, int *nelem) { + + // we allocate one extra byte into which we never read so there will always be a NULL terminating the string + ALLOCATE_PTR (buffer, char, NBUFFER + 1); + bzero (buffer, NBUFFER + 1); + + int Nstart = 0; // location of the last valid byte in the buffer (start filling here) + int EndOfFile = FALSE; + + // int TimeFormat = TIME_MJD; + // time_t TimeReference = 0; + + int Nelem = 0; + int NELEM = 1000; + ALLOCATE_PTR (elem, unsigned int, NELEM); + + while (!EndOfFile) { + + int Nbytes = NBUFFER - Nstart; + bzero (&buffer[Nstart], Nbytes + 1); + // ^-- we have allocated one extra byte into which we never read so there will always be a NULL terminating the string + + int nread = fread (&buffer[Nstart], 1, Nbytes, f); + if (ferror (f)) { + perror ("error reading data file"); + return NULL; + } + + // we still need to parse the rest of the buffer, but there might not be an EOL on the last line + if (nread == 0) { + EndOfFile = TRUE; + } + + int bufferStatus = TRUE; + + char *c0 = buffer; // c0 always marks the start of a line + while (bufferStatus) { + char *c1 = strchr (c0, '\n'); // find the end of this current line (also valid for a Mac: \r\n) + if (!c1) { + c1 = strchr (c0, '\r'); // try \r for Windows files + } + if (!c1) { + Nstart = strlen (c0); + if (EndOfFile) { + // if we have reached EOF, we need to do one last pass in case there is a line without a return + c1 = c0 + Nstart; + bufferStatus = FALSE; + if (Nstart == 0) continue; // if we have reached EOF and c0 points at the last valid character, we are done + } else { + // if we have not reached EOF, we need to shift the buffer to the start of this line and read more data + memmove (buffer, c0, Nstart); + bufferStatus = FALSE; + continue; + } + } + *c1 = 0; // mark the end of the line + + // skip to the next line (but if EOF, do not overrun buffer) + if (*c0 == '#') { if (!EndOfFile) { c0 = c1 + 1; } continue; } + if (*c0 == '!') { if (!EndOfFile) { c0 = c1 + 1; } continue; } + + // the times are written as MJD in the file: + double dvalue; + int readStatus = IsCSV ? dparse_csv (&dvalue, COLUMN_NUMBER, c0) : dparse (&dvalue, COLUMN_NUMBER, c0); + if (!readStatus) continue; // skip invalid entries + + // convert the (double) mjd value to (unsigned int) time + elem[Nelem] = ohana_mjd_to_sec (dvalue); + + Nelem ++; + if (Nelem == NELEM) { + NELEM += 1000; + REALLOCATE (elem, unsigned int, NELEM); + } + + if (!EndOfFile) { + c0 = c1 + 1; + } + } + } + FREE (buffer); + + *nelem = Nelem; + return elem; +} + +/* + time_t tvalue; + int readStatus = IsCSV ? tparse_csv (&tvalue, COLUMN_NUMBER, c0) : tparse (&tvalue, COLUMN_NUMBER, c0); + if (!readStatus) continue; // skip invalid entries + + dvalue = TimeValue (tvalue, TimeReference, TimeFormat); + elem[Nelem] = (int) dvalue; +*/ Index: /branches/eam_branches/relphot.20200519/src/StarOps.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/StarOps.c (revision 41361) +++ /branches/eam_branches/relphot.20200519/src/StarOps.c (revision 41362) @@ -105,4 +105,5 @@ results->Ncal = 0; results->Nmos = 0; + results->Ngrp = 0; results->Ngrid = 0; } @@ -264,4 +265,5 @@ summary->Ncal += results->Ncal ; summary->Nmos += results->Nmos ; + summary->Ngrp += results->Ngrp ; summary->Ngrid += results->Ngrid; } @@ -313,5 +315,5 @@ SetMrelInfoAccum (&summary, &results); } - if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrp, summary.Ngrid, summary.Nsys); SetMrelInfoFree (&results); @@ -336,5 +338,5 @@ SetMrelInfoAccum (&summary, &results); } - if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrp, summary.Ngrid, summary.Nsys); SetMrelInfoFree (&results); @@ -393,5 +395,5 @@ // report stats & summary from the threads for (i = 0; i < NTHREADS; i++) { - if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", + if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", i, threadinfo[i].summary.Ncode, @@ -400,4 +402,5 @@ threadinfo[i].summary.Ncal, threadinfo[i].summary.Nmos, + threadinfo[i].summary.Ngrp, threadinfo[i].summary.Ngrid, threadinfo[i].summary.Nsys); @@ -405,5 +408,5 @@ SetMrelInfoFree (&threadinfo[i].summary); } - if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrp, summary.Ngrid, summary.Nsys); free (threadinfo); @@ -488,4 +491,6 @@ float Mmos = getMmos (m, i); if (isnan(Mmos)) continue; + float Mgrp = getMgrp (m, i); + if (isnan(Mgrp)) continue; float Mgrid = getMgrid (m, i); if (isnan(Mgrid)) continue; @@ -502,6 +507,7 @@ // set the output calibration - catalog[i].measure[m].McalPSF = McalPSF + Mmos + Mgrid; - catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Mmos + Mgrid : McalPSF + Mmos + Mgrid; + float Moff = Mmos + Mgrp + Mgrid; + catalog[i].measure[m].McalPSF = McalPSF + Moff; + catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Moff : McalPSF + Moff; if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { @@ -662,6 +668,6 @@ off_t *ilist; double *tlist, *list, *dlist; - float Msys, Mcal, Mmos, Mgrid; - int Ncal, Nmos, Ngrid, Nfew; + float Msys, Mcal, Mmos, Mgrp, Mgrid; + int Ncal, Nmos, Ngrp, Ngrid, Nfew; int Nsecfilt = GetPhotcodeNsecfilt (); @@ -689,5 +695,5 @@ Ndel = Nave = 0; - Ncal = Nmos = Ngrid = Nfew = 0; + Ncal = Nmos = Ngrp = Ngrid = Nfew = 0; for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { @@ -724,4 +730,6 @@ Mmos = getMmos (m, i); if (isnan(Mmos)) { Nmos ++; continue; } + Mgrp = getMgrp (m, i); + if (isnan(Mgrp)) { Ngrp ++; continue; } Mgrid = getMgrid (m, i); if (isnan(Mgrid)) { Ngrid ++; continue; } Index: /branches/eam_branches/relphot.20200519/src/TGroupOps.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/TGroupOps.c (revision 41362) +++ /branches/eam_branches/relphot.20200519/src/TGroupOps.c (revision 41362) @@ -0,0 +1,1186 @@ +# include "relphot.h" + +// XXX is this needed? +void plot_setMcal (double *list, int Npts, StatType *stats, float clouds); + +// see discussion in ImagesOps.c re: IDX_T + +// array of tgroup definition structures +static off_t Ntgroup; +static TGroup *tgroup; + +// tgroup index for given image : ImageToTGroup[ImageIndex] = TGroupIndex (ImageIndex : 0 < Nimage) +static off_t *ImageToTGroup; + +// elsewhere, we have loaded a set of catalogs with measures (catalog[cat].measure[meas]) + +// each tgroup has N_onTGroup[TGroupIndex] measurements +// static off_t *N_onTGroup; // actual number of measurements on tgroup +// static off_t *N_ONTGROUP; // allocated number of measurements on tgroup + +// relationships between the measure,catalog set and the tgroups: +static off_t **MeasureToTGroup; // TGroup index from measure,catalog : MeasureToTGroup[cat][meas] = TGroupIndex +// static off_t **TGroupToCatalog; // catalog for given measure on tgroup : TGroupCatalog[TGroupIndex][i] = cat (i : 0 < NonTGroup[TGroupIndex]) +// static off_t **TGroupToMeasure; // measure for given measure on tgroup : TGroupMeasure[TGroupIndex][i] = cat (i : 0 < NonTGroup[TGroupIndex]) + +// tgroupTimes carries the times of the photometric nights +static int NtgroupTimes; +static unsigned int *tgroupTimes; + +void sort_times (unsigned int *T, int N); + +/* find tgroups (unique time periods) (NOTE : we do NOT require matching photcodes...) + we need to have a tgroup for each photcode available for a given photometric night + */ + +// load the photometric nights from a file : these will be used to define the tgroups +// XXX call this in args.c where TGROUP_ZEROPT is set +void loadTGroups (char *filename) { + + if (!TGROUP_ZEROPT) return; + + // read start times from a file (currently a list of MJD values) + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "tgroup file %s not found\n", filename); + exit (4); + } + tgroupTimes = ReadTGroupFile (f, &NtgroupTimes); + if (!tgroupTimes) { + fprintf (stderr, "error reading tgroup file %s\n", filename); + exit (4); + } + + // sort the times for quick bisection + sort_times (tgroupTimes, NtgroupTimes); + + // we define a tgroup for each photcode being considered + // below, if there are no images in this tgroup with a given + // photcode, it is left empty (no images) + Ntgroup = NtgroupTimes*Nphotcodes; + + ALLOCATE (tgroup, TGroup, Ntgroup); + + // the start of an MJD day is 14:00 HST, so the end time can be the start time + 23h 59m 59s (86399s) + // XXX NOTE that this is HST-centric + + for (int i = 0; i < NtgroupTimes; i++) { + for (int j = 0; j < Nphotcodes; j++) { + /* a new tgroup, define ranges */ + int Nt = i*Nphotcodes + j; + tgroup[Nt].start = tgroupTimes[i]; + tgroup[Nt].stop = tgroupTimes[i] + 86399; + tgroup[Nt].McalPSF = 0.0; // note : at the end, tgroup.Mcal is added back to the input images + tgroup[Nt].McalAPER = 0.0; // note : tgroup stores only offsets relative to the original image values + tgroup[Nt].dMcal = 0.0; // note : at the end, tgroup.Mcal is added back to the input images + tgroup[Nt].dMsys = 0.0; + tgroup[Nt].McalChiSq = 0.0; // NAN or 0.0? + tgroup[Nt].flags = 0; + tgroup[Nt].photcode = photcodes[j][0].code; + + tgroup[Nt].NIMAGE = 1000; + tgroup[Nt].Nimage = 0; + ALLOCATE (tgroup[Nt].image, off_t, tgroup[Nt].NIMAGE); + tgroup[Nt].image[0] = -1; + + tgroup[Nt].Nmeasure = 0; + tgroup[Nt].NMEASURE = 1000; + ALLOCATE (tgroup[Nt].measure, off_t, tgroup[Nt].NMEASURE); + ALLOCATE (tgroup[Nt].catalog, off_t, tgroup[Nt].NMEASURE); + } + } + return; +} + +// assign the images to the tgroups +void initTGroups (Image *subset, off_t *LineNumber, off_t Nsubset) { + + if (!TGROUP_ZEROPT) return; + + /* a 'tgroup' in relphot is a virtual concept: there is no + * entry in the image table that represents this tgroup. Instead, it is an + * internal construct that defines a group of related images + */ + + ALLOCATE (ImageToTGroup, off_t, Nsubset); // tgroup to which image belongs + + // assign each image to a tgroup + int Nsimple = 0; + for (off_t i = 0; i < Nsubset; i++) { + ImageToTGroup[i] = -1; + + // ignore non-GPC1 images + if (!isGPC1chip(subset[i].photcode)) { + Nsimple ++; + continue; + } + + off_t j = findTGroup(subset[i].tzero, subset[i].photcode); + if (j == -1) { + // images which do not belong to a tgroup should be ignored in this analysis + // XXX mark the images and the measurements somehow + continue; + } + + // add reference from image to tgroup + ImageToTGroup[i] = j; + + // add image to tgroup image list + tgroup[j].image[tgroup[j].Nimage] = i; // should this be 'i' or LineNumber[i] ?? + tgroup[j].Nimage ++; + if (tgroup[j].Nimage == tgroup[j].NIMAGE) { + tgroup[j].NIMAGE += 1000; + REALLOCATE (tgroup[j].image, off_t, tgroup[j].NIMAGE); + } + } + + initTGroupsMcal (); + + fprintf (stderr, "matched %d images to %d tgroups, %d simple chips not matched to tgroups\n", (int) (Nsubset - Nsimple), (int) Ntgroup, (int) Nsimple); + return; +} + +void initTGroupsMcal (void) { + + off_t Nimages; + off_t *LineNumber; + Image *images = getimages (&Nimages, &LineNumber); + + if (!TGROUP_ZEROPT) return; + + StatType stats; + liststats_setmode (&stats, "MEDIAN"); + + // init the tgroup array values + for (int i = 0; i < Ntgroup; i++) { + tgroup[i].McalPSF = 0.0; // note : at the end, tgroup.Mcal is added back to the input images + tgroup[i].McalAPER = 0.0; // note : tgroup stores only offsets relative to the original image values + tgroup[i].dMcal = 0.0; // note : at the end, tgroup.Mcal is added back to the input images + tgroup[i].dMsys = 0.0; + tgroup[i].McalChiSq = 0.0;// NAN or 0.0? + + if (tgroup[i].Nimage == 0) continue; // no images, ignore the tgroup + + // calculate the median of the Mcal values for the images in this tgroup: + + ALLOCATE_PTR (McalPSF, double, tgroup[i].Nimage); + ALLOCATE_PTR (McalAPER, double, tgroup[i].Nimage); + for (int j = 0; j < tgroup[i].Nimage; j++) { + int Nsub = LineNumber[tgroup[i].image[j]]; + McalPSF[j] = images[Nsub].McalPSF; + McalAPER[j] = images[Nsub].McalAPER; + + // set these back to 0.0 to have no future effect + images[Nsub].McalPSF = 0.0; + images[Nsub].McalAPER = 0.0; + } + + liststats (McalPSF, NULL, NULL, tgroup[i].Nimage, &stats); + tgroup[i].McalPSF = stats.median; + FREE (McalPSF); + + liststats (McalAPER, NULL, NULL, tgroup[i].Nimage, &stats); + tgroup[i].McalAPER = stats.median; + FREE (McalAPER); + } + return; +} + +TGroup *getTGroupForImage (off_t im) { + + if (im < 0) return NULL; + if (!ImageToTGroup) return NULL; + + off_t m = ImageToTGroup[im]; + if (m < 0) return NULL; + if (m >= Ntgroup) return NULL; + + return (&tgroup[m]); +} + +// use bisection to find the overlapping tgroup (returns exact match) +// tgroupTimes is a sorted, unique list of times +// assume a 1day range for now +off_t findTGroup (unsigned int start, int photcode) { + + // first find the matching photcode + int Ns = -1; + for (int i = 0; (Ns < 0) && (i < Nphotcodes); i++) { + int ecode = GetPhotcodeEquivCodebyCode (photcode); + if (photcodes[i][0].code != ecode) continue; + Ns = i; + } + if (Ns < 0) return -1; + + off_t Nlo, Nhi, N; + + // find the last tgroup before start + Nlo = 0; // first valid tgroupTimes value + Nhi = Ntgroup - 1; // last valid tgroupTimes value + + // if start is not in this range, return -1 + if (start < tgroupTimes[Nlo]) return (-1); + if (start > tgroupTimes[Nhi]) return (-1); + + while (Nhi - Nlo > 4) { + N = 0.5*(Nlo + Nhi); + if (tgroupTimes[N] < start) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N, Ntgroup - 1); + } + } + // we now have : tgroupTimes[Nlo] < start <= tgroupTimes[Nhi] + + // find a matched tgroup starting from Nlo, or return -1 + // tgroupTimes is a lower bound, the upper bound is tgroupTimes + 86399 (1d) + for (N = Nlo; N <= Nhi; N++) { + if (start < tgroupTimes[N]) continue; + if (start > tgroupTimes[N] + 86399) continue; // XXX since they are sorted, probably do not need this check + return (N*Nphotcodes + Ns); // we have the tgroupTime and photcodes, so we can construct the tgroup index + } + return (-1); +} + +void setMcalFromTGroups () { + + off_t i, j, im, Nimage; + Image *image; + + if (!TGROUP_ZEROPT) return; + + image = getimages (&Nimage, NULL); + + fprintf (stderr, "*** return Mcal from tgroup.Mcal to image.Mcal ***\n"); + + // copy the tgroup results to the images. set the tgroup Mcal to 0.0 since we have moved its + // impact to the images + for (i = 0; i < Ntgroup; i++) { + for (j = 0; j < tgroup[i].Nimage; j++) { + im = tgroup[i].image[j]; + image[im].McalPSF += tgroup[i].McalPSF; + image[im].McalAPER += tgroup[i].McalAPER; + image[im].dMcal = tgroup[i].dMcal; + image[im].McalChiSq = tgroup[i].McalChiSq; + image[im].dMagSys = tgroup[i].dMsys; + image[im].nFitPhotom = tgroup[i].nFitPhotom; + image[im].flags |= (tgroup[i].flags & ID_IMAGE_PHOTOM_FEW); + image[im].flags |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR); + } + tgroup[i].McalPSF = 0.0; + tgroup[i].McalAPER = 0.0; + } +} + +void initTGroupBins (Catalog *catalog, int Ncatalog) { + + off_t i, j; + + /* measure -> tgroup */ + if (!TGROUP_ZEROPT) return; + + ALLOCATE (MeasureToTGroup, off_t *, Ncatalog); + for (i = 0; i < Ncatalog; i++) { + ALLOCATE (MeasureToTGroup[i], off_t, MAX (catalog[i].Nmeasure, 1)); + for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToTGroup[i][j] = -1; + } +} + +void freeTGroupBins (int Ncatalog, int doTGroupList) { + + off_t i; + + /* measure -> tgroup */ + if (!TGROUP_ZEROPT) return; + + for (i = 0; i < Ncatalog; i++) { + free (MeasureToTGroup[i]); + } + free (MeasureToTGroup); +} + +int findTGroups (Catalog *catalog, int Ncatalog) { + + int i, ecode, found, Ns; + off_t j; + + if (!TGROUP_ZEROPT) return (FALSE); + // if we are calibrating by tgroup, redefine myDet == on one of my tgroups + + int Nmatch = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Nmeasure; j++) { + catalog[i].measureT[j].myDet = FALSE; // a detetion is not mine until proven otherwise + + if (TimeSelect) { + if (catalog[i].measureT[j].t < TSTART) continue; + if (catalog[i].measureT[j].t > TSTOP) continue; + } + + ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[j].photcode); + found = FALSE; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (ecode == photcodes[Ns][0].code) found = TRUE; + } + if (!found) continue; + + matchTGroups (catalog, j, i); + Nmatch ++; + } + } + // fprintf (stderr, "matched %d detections to tgroups\n", Nmatch); + return (TRUE); +} + +void matchTGroups (Catalog *catalog, off_t meas, int cat) { + + MeasureTiny *measure = &catalog[cat].measureT[meas]; + + off_t ID = measure[0].imageID; // ID is a unique ID in DVO + off_t idx = getImageByID (ID); // idx is the sequence number in the subset list + if (idx == -1) { + if (VERBOSE2) fprintf (stderr, "missed measurement "OFF_T_FMT", %d\n", meas, cat); + return; + } + + off_t grpID = ImageToTGroup[idx]; + if (grpID < 0) { + // Image *image = getimage(idx); + // fprintf (stderr, "unmatched image %s\n", image[0].name); + // skip measurements from simple chips (not matched to a tgroup by definition) + return; + } + + // test to check we got the right match: + { + Image *image = getimage(idx); + // XXX we are now matching with just tzero. be careful for cameras with drift + // unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1); + unsigned int imageStart = image[0].tzero; + if ((imageStart < tgroup[grpID].start) || (imageStart > tgroup[grpID].stop)) { + fprintf (stderr, "error in image to tgroup match\n"); + abort(); + } + } + + // this measurement is on one of my tgroups, mark it as mine. + catalog[cat].measureT[meas].myDet = TRUE; + MeasureToTGroup[cat][meas] = grpID; + + tgroup[grpID].catalog[tgroup[grpID].Nmeasure] = cat; + tgroup[grpID].measure[tgroup[grpID].Nmeasure] = meas; + tgroup[grpID].Nmeasure ++; + + if (tgroup[grpID].Nmeasure == tgroup[grpID].NMEASURE) { + tgroup[grpID].NMEASURE += 1000; + REALLOCATE (tgroup[grpID].catalog, off_t, tgroup[grpID].NMEASURE); + REALLOCATE (tgroup[grpID].measure, off_t, tgroup[grpID].NMEASURE); + } + return; +} + +float getMgrp (off_t meas, int cat) { + + off_t i; + float value; + + if (!TGROUP_ZEROPT) return (0); + + // unassigned measurements belong to simple chips + i = MeasureToTGroup[cat][meas]; + if (i == -1) return (0.0); + + // if (tgroup[i].flags & IMAGE_BAD) return (NAN); + value = tgroup[i].McalPSF; + return (value); +} + +int getTGroupFlags (off_t meas, int cat) { + + if (!TGROUP_ZEROPT) return (0); + + // unassigned measurements belong to simple chips + int i = MeasureToTGroup[cat][meas]; + if (i == -1) return (0); + + return (tgroup[i].flags); +} + +typedef struct { + int Nfew; + int Nbad; + int Ncal; + int Nmos; + int Ngrid; + int Nrel; + int Nsys; + int Nskip; + off_t Nmax; + int PoorImages; + double *psfMagList; + double *psfErrList; + double *brightMagList; + double *brightErrList; + double *kronMagList; + double *kronErrList; +} SetMgrpInfo; + +enum {THREAD_RUN, THREAD_DONE}; + +typedef struct { + int entry; + int state; + Catalog *catalog; + Image *image; + FlatCorrectionTable *flatcorr; + SetMgrpInfo info; +} ThreadInfo; + +int setMgrp_tgroup (TGroup *tgroup, off_t Nmos, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr); +void *setMgrp_worker (void *data); +int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr); + +void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists, int PoorImages) { + info->Nfew = 0; + info->Nbad = 0; + info->Ncal = 0; + info->Nmos = 0; + info->Nrel = 0; + info->Ngrid = 0; + info->Nskip = 0; + info->Nsys = 0; + + info->Nmax = Nmax; + info->PoorImages = PoorImages; + + if (allocLists) { + ALLOCATE (info->psfMagList, double, Nmax); + ALLOCATE (info->psfErrList, double, Nmax); + ALLOCATE (info->kronMagList, double, Nmax); + ALLOCATE (info->kronErrList, double, Nmax); + ALLOCATE (info->brightMagList, double, Nmax); + ALLOCATE (info->brightErrList, double, Nmax); + } +} + +void SetMgrpInfoFree (SetMgrpInfo *info) { + free (info->psfMagList ); + free (info->psfErrList ); + free (info->kronMagList ); + free (info->kronErrList ); + free (info->brightMagList); + free (info->brightErrList); +} + +void SetMgrpInfoAccum (SetMgrpInfo *summary, SetMgrpInfo *results) { + summary->Nfew += results->Nfew ; + summary->Nsys += results->Nsys ; + summary->Nbad += results->Nbad ; + summary->Ncal += results->Ncal ; + summary->Nmos += results->Nmos ; + summary->Nrel += results->Nrel ; + summary->Ngrid += results->Ngrid; + summary->Nskip += results->Nskip; +} + +static int npass_output = 0; + +// mutex to lock setMgrp_worker operations +static pthread_mutex_t setMgrp_mutex = PTHREAD_MUTEX_INITIALIZER; +static int nextTGroup = 0; + +// we have an array of tgroups (tgroup, Ntgroup). we need to hand out tgroups one at a time to +// the worker threads as they need +off_t getNextTGroupForThread () { + + pthread_mutex_lock (&setMgrp_mutex); + if (nextTGroup >= Ntgroup) { + pthread_mutex_unlock (&setMgrp_mutex); + return (-1); + } + int thisTGroup = nextTGroup; + nextTGroup ++; + + pthread_mutex_unlock (&setMgrp_mutex); + return (thisTGroup); +} + +int setMgrp (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { + + off_t i, N, Nmax; + Image *image; + + if (!TGROUP_ZEROPT) return (FALSE); + + // plots cannot be done in a threaded context (the plot commands collide) + // so do not run setMgrp in threaded mode if PLOTSTUFF is set + if (NTHREADS && !PLOTSTUFF) { + int status = setMgrp_threaded (catalog, PoorImages, flatcorr); + return status; + } + + image = getimages (&N, NULL); + + fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); + + if (PoorImages) { + // XXX use bad stars and measurements for PoorImages? or not? + // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; + IMAGE_BAD = 0; + } + + Nmax = 0; + for (i = 0; i < Ntgroup; i++) { + Nmax = MAX (Nmax, tgroup[i].Nmeasure); + } + + SetMgrpInfo info; + SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages); + + // int savePlotDelay = PLOTDELAY; + // if (PLOTSTUFF) PLOTDELAY = 0.0; + + for (i = 0; i < Ntgroup; i++) { + setMgrp_tgroup (&tgroup[i], i, image, catalog, &info, flatcorr); + } + SetMgrpInfoFree (&info); + + if (PLOTSTUFF) { + fprintf (stdout, "press return\n"); + if (fscanf (stdin, "%*c") != 1) fprintf (stderr, "\n"); + } + + npass_output ++; + + fprintf (stderr, "%d tgroups marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", info.Nfew, info.Nbad, info.Ncal, info.Ngrid, info.Nrel, info.Nsys); + + if (PoorImages) { + IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; + STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; + MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; + } + return (TRUE); +} + +# define NIGHT_TOOFEW 100 +# define NIGHT_GOOD_FRACTION 0.1 + +// 'tgroup' is a pointer to the current tgroup of interest (Nmos) +int setMgrp_tgroup (TGroup *myTGroup, off_t Ngrp, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr) { + + off_t j, NimageReal; + + StatType stats; + liststats_setmode (&stats, "INNER_WTMEAN"); + + double *psfMagList = info->psfMagList; + double *psfErrList = info->psfErrList; + double *kronMagList = info->kronMagList; + double *kronErrList = info->kronErrList; + double *brightMagList = info->brightMagList; + double *brightErrList = info->brightErrList; + + // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used + getimages (&NimageReal, NULL); + + /* on PoorImages run, skip good images */ + if (info->PoorImages) { + int bad = myTGroup[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP); + if (!bad) return TRUE; + } + + int Nsecfilt = GetPhotcodeNsecfilt (); + + int testImage = FALSE; + // testImage |= (abs(myTGroup[0].start - 1324104046) < 10); + // testImage |= (abs(myTGroup[0].start - 1324103823) < 10); + + FILE *fout = NULL; + if (testImage) { + char filename[64]; + snprintf (filename, 64, "test.%05d.%02d.dat", (int) Ngrp, npass_output); + fout = fopen (filename, "w"); + } + + // number of stars to measure the bright-end scatter + int Nbright = 0; + + int N = 0; + for (j = 0; j < myTGroup->Nmeasure; j++) { + + off_t m = myTGroup->measure[j]; + off_t c = myTGroup->catalog[j]; + + // NOTE : we are only using Mcal == McalPSF in this function; we set McalAPER to McalPSF + if (fout) { + + float Mcal = getMcal (m, c, MAG_CLASS_PSF); // image zero point + float Mmos = getMmos (m, c); // mosaic zero point + // XXX these two should both be 0.0 if we are fitting tgroups : skip those calls above? + + float Mgrid = getMgrid (m, c); // camera offset (deprecated?) + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude + float Mflat = getMflat (m, c, flatcorr, catalog); // flat-field correction + + off_t n = catalog[c].measureT[m].averef; + + // magnitude for this measurement + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + + // for historical reasons, Mflat is defined with the wrong sign + float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid + Mflat; + + int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD); + + fprintf (fout, "%f %f : %f %f %f %f %f : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad); + } + + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { + info->Nbad ++; + continue; + } + float Mcal = getMcal (m, c, MAG_CLASS_PSF); + if (isnan(Mcal)) { + info->Ncal++; + continue; + } + float Mmos = getMmos (m, c); + if (isnan(Mmos)) { + info->Nmos++; + continue; + } + float Mgrid = getMgrid (m, c); + if (isnan(Mgrid)) { + info->Ngrid ++; + continue; + } + + // Mrel* is the average magnitude for this star. For PS1 stacks, we have too much + // PSF variability. We need to calibrate the PSF magnitudes separately from the + // Aperture-like magnitues. (We have an option to use the kron magnitudes or the + // other apertures here). I basically need to do this analysis separately for each + // magnitude type + + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(MrelPSF)) { + info->Nrel ++; + continue; + } + float MrelKron = getMrel (catalog, m, c, MAG_CLASS_KRON, MAG_SRC_CHP); + + // image.Mcal is not supposed to include the flat-field correction, so we need to + // apply that offset as well here for this image (in other words, each detection is + // being compared to the model, excluding the zero point, Mcal. The model includes + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) + + float Mflat = getMflat (m, c, flatcorr, catalog); + + off_t n = catalog[c].measureT[m].averef; + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + if (isnan(MsysPSF)) { + info->Nsys++; + continue; + } + float MsysKron = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_KRON); + + PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); + if (!code) goto skip; + if (code->equiv < 1) goto skip; + int Nsec = GetPhotcodeNsec (code->equiv); + if (Nsec == -1) goto skip; + + skip: + assert (N < info->Nmax); + assert (N >= 0); + assert (Nbright < info->Nmax); + assert (Nbright >= 0); + + psfMagList[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat; + psfErrList[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + kronMagList[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat; + kronErrList[N] = psfErrList[N]; + if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { + brightMagList[Nbright] = psfMagList[N]; + brightErrList[Nbright] = psfErrList[N]; + Nbright ++; + } + N++; + } + if (fout) { fclose (fout); } + + /* too few good measurements or too many bad measurements (skip in PoorImages run) */ + if (!info->PoorImages) { + int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure); + if (mark) { + // XXX is this right (effect of subset?) : image[tgroup[Ngrp].image[0]].name + if (VERBOSE2) { fprintf (stderr, "marked tgroup %s ("OFF_T_FMT"), (%d < %d) || (%d < %f*"OFF_T_FMT")\n", image[tgroup[Ngrp].image[0]].name, Ngrp, N, NIGHT_TOOFEW, N, NIGHT_GOOD_FRACTION, myTGroup->Nmeasure); } + myTGroup[0].flags |= ID_IMAGE_PHOTOM_FEW; + info->Nfew ++; + if (testImage) { + fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) NIGHT_TOOFEW, (int) (NIGHT_GOOD_FRACTION*myTGroup->Nmeasure)); + } + } else { + myTGroup[0].flags &= ~ID_IMAGE_PHOTOM_FEW; + } + } + + liststats (psfMagList, psfErrList, NULL, N, &stats); + if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mgrp: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N); + + // for now, I have no reason to measure these separately for camera-level images + myTGroup[0].McalPSF = stats.mean; + myTGroup[0].dMcal = stats.error; + myTGroup[0].McalChiSq = stats.chisq; + myTGroup[0].nFitPhotom = N; + + // drop this for tgroups? + liststats (kronMagList, kronErrList, NULL, N, &stats); + myTGroup[0].McalAPER = stats.mean; + + if (testImage) { + fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Ngrp, myTGroup[0].start, stats.mean, stats.error, myTGroup[0].nFitPhotom); + } + + if (PLOTSTUFF) { + fprintf (stderr, "Mgrp: %6.3f %6.3f +/- %6.3f %5d %5d | %s\n", stats.mean, stats.median, stats.sigma, stats.Nmeas, N, image[tgroup[Ngrp].image[0]].name); + plot_setMcal (psfMagList, N, &stats, CLOUD_TOLERANCE); + } + + // bright end scatter + liststats (brightMagList, brightErrList, NULL, Nbright, &stats); + myTGroup[0].dMsys = stats.sigma; + + // keep this?? + if (myTGroup[0].McalPSF < -CLOUD_TOLERANCE) { + myTGroup[0].McalPSF = 0.0; + } + + if (testImage) { + fprintf (stderr, "%f %f : %f\n", myTGroup[0].McalPSF, myTGroup[0].dMsys, myTGroup[0].McalChiSq); + } + + return TRUE; +} + +int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { + + int i; + off_t N; + + Image *image = getimages (&N, NULL); + + fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); + + if (PoorImages) { + // XXX use bad stars and measurements for PoorImages? or not? + // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; + IMAGE_BAD = 0; + } + + off_t Nmax = 0; + for (i = 0; i < Ntgroup; i++) { + Nmax = MAX (Nmax, tgroup[i].Nmeasure); + } + + SetMgrpInfo summary; + SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages); + + pthread_attr_t attr; + pthread_attr_init (&attr); + pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED); + + pthread_t *threads; + ALLOCATE (threads, pthread_t, NTHREADS); + + ThreadInfo *threadinfo; + ALLOCATE (threadinfo, ThreadInfo, NTHREADS); + + // each time this function is called, we cycle through the available tgroups + // make sure we start at 0 + nextTGroup = 0;; + + // launch N worker threads + for (i = 0; i < NTHREADS; i++) { + threadinfo[i].entry = i; + threadinfo[i].state = THREAD_RUN; + threadinfo[i].catalog = catalog; + threadinfo[i].image = image; + threadinfo[i].flatcorr = flatcorr; + SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages); + pthread_create (&threads[i], NULL, setMgrp_worker, &threadinfo[i]); + } + pthread_attr_destroy (&attr); + + // wait until all threads have finished + while (1) { + int allDone = TRUE; + for (i = 0; i < NTHREADS; i++) { + if (threadinfo[i].state == THREAD_RUN) allDone = FALSE; + } + if (allDone) { + break; + } + usleep (500000); + } + + // all threads are done, free the threads array and grab the info + free (threads); + + // report stats & summary from the threads + for (i = 0; i < NTHREADS; i++) { + fprintf (stderr, "setMgrp thread %d : %d tgroups marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d), %d partials skipped\n", + i, + threadinfo[i].info.Nfew, + threadinfo[i].info.Nbad, + threadinfo[i].info.Ncal, + threadinfo[i].info.Nmos, + threadinfo[i].info.Ngrid, + threadinfo[i].info.Nrel, + threadinfo[i].info.Nsys, + threadinfo[i].info.Nskip); + SetMgrpInfoAccum (&summary, &threadinfo[i].info); + } + fprintf (stderr, "total : %d tgroups marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d), %d partials skipped\n", + summary.Nfew, + summary.Nbad, + summary.Ncal, + summary.Nmos, + summary.Ngrid, + summary.Nrel, + summary.Nsys, + summary.Nskip); + free (threadinfo); + + npass_output ++; + + if (PoorImages) { + IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; + STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; + MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; + } + return TRUE; +} + +void *setMgrp_worker (void *data) { + + ThreadInfo *threadinfo = data; + + SetMgrpInfo results; + SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here + + while (1) { + + off_t i = getNextTGroupForThread(); + if (i == -1) { + threadinfo->state = THREAD_DONE; + return NULL; + } + + Catalog *catalog = threadinfo->catalog; + FlatCorrectionTable *flatcorr = threadinfo->flatcorr; + Image *image = threadinfo->image; + + setMgrp_tgroup (&tgroup[i], i, image, catalog, &results, flatcorr); + SetMgrpInfoAccum (&threadinfo->info, &results); + } + + SetMgrpInfoFree (&results); + return NULL; +} + +StatType statsTGroupM (Catalog *catalog) { + OHANA_UNUSED_PARAM(catalog); + + off_t i; + int n; + double *list, *dlist; + StatType stats; + + bzero (&stats, sizeof (StatType)); + if (!TGROUP_ZEROPT) return (stats); + + ALLOCATE (list, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + n = 0; + for (i = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + list[n] = tgroup[i].McalPSF; + dlist[n] = 1; + n++; + } + + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); + free (list); + free (dlist); + return (stats); +} + +StatType statsTGroupdM (Catalog *catalog) { + OHANA_UNUSED_PARAM(catalog); + + off_t i, n; + double *list, *dlist; + StatType stats; + + bzero (&stats, sizeof (StatType)); + if (!TGROUP_ZEROPT) return (stats); + + ALLOCATE (list, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + n = 0; + for (i = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + list[n] = tgroup[i].dMcal; + dlist[n] = 1; + n++; + } + + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); + free (list); + free (dlist); + return (stats); +} + +StatType statsTGroupN (Catalog *catalog, FlatCorrectionTable *flatcorr) { + + off_t i, j, m, c, n, N; + double *list, *dlist; + float Mcal, Mgrid, Mrel, Mmos; + StatType stats; + + bzero (&stats, sizeof (StatType)); + if (!TGROUP_ZEROPT) return (stats); + + ALLOCATE (list, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + n = 0; + for (i = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + + N = 0; + for (j = 0; j < tgroup[i].Nmeasure; j++) { + + m = tgroup[i].measure[j]; + c = tgroup[i].catalog[j]; + + Mcal = getMcal (m, c, MAG_CLASS_PSF); + if (isnan(Mcal)) continue; + Mmos = getMmos (m, c); + if (isnan(Mmos)) continue; + Mgrid = getMgrid (m, c); + if (isnan(Mgrid)) continue; + Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(Mrel)) continue; + N++; + } + list[n] = N; + dlist[n] = 1; + n++; + } + // fprintf (stderr, "Ntgroup: "OFF_T_FMT", n: "OFF_T_FMT"\n", Ntgroup, n); + + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); + free (list); + free (dlist); + return (stats); +} + +StatType statsTGroupX (Catalog *catalog) { + OHANA_UNUSED_PARAM(catalog); + + off_t i, n; + double *list, *dlist; + StatType stats; + + bzero (&stats, sizeof (StatType)); + if (!TGROUP_ZEROPT) return (stats); + + ALLOCATE (list, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + n = 0; + for (i = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + + list[n] = tgroup[i].McalChiSq; + dlist[n] = 1; + n++; + } + + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); + free (list); + free (dlist); + return (stats); +} + +/* mark tgroup if: abs(Mcal - ) too large, dMcal too large */ +void clean_tgroups () { + + off_t i, N, mark, Nmark, Nscatter, Noffset; + double *mlist, *slist, *dlist; + double MaxOffset, MedOffset, MaxScatter; + + if (!TGROUP_ZEROPT) return; + + if (VERBOSE) fprintf (stderr, "marking poor tgroups\n"); + + ALLOCATE (mlist, double, Ntgroup); + ALLOCATE (slist, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + for (i = N = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + mlist[N] = tgroup[i].McalPSF; + slist[N] = tgroup[i].dMcal; + dlist[N] = 1; + N++; + } + + StatType stats; + liststats_setmode (&stats, "MEAN"); + + liststats (mlist, dlist, NULL, N, &stats); + MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma); + MedOffset = stats.median; + + liststats (slist, dlist, NULL, N, &stats); + MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median); + fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset); + + Nmark = Nscatter = Noffset = 0; + for (i = 0; i < Ntgroup; i++) { + // if we are keeping ubercal sacrosanct, then we should not be allowed to break them... + if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + + if (tgroup[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue; + + mark = FALSE; + if (tgroup[i].dMcal > MaxScatter) { + mark = TRUE; + Nscatter ++; + } + if (fabs(tgroup[i].McalPSF - MedOffset) > MaxOffset) { + mark = TRUE; + Noffset ++; + } + if (mark) { + Nmark ++; + tgroup[i].flags |= ID_IMAGE_PHOTOM_POOR; + } else { + tgroup[i].flags &= ~ID_IMAGE_PHOTOM_POOR; + } + } + + fprintf (stderr, OFF_T_FMT" tgroups marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n", Nmark, Nscatter, Noffset); + + free (mlist); + free (slist); + free (dlist); +} + +void plot_tgroup_fields (Catalog *catalog) { + + off_t i, j, m, c, N, Nimage; + double *xlist, *ylist; + char string[64]; + Graphdata graphdata; + + if (!TGROUP_ZEROPT) return; + + // Image *image = getimages (&Nimage, NULL); returned value ignored + getimages (&Nimage, NULL); + + N = 0; + for (i = 0; i < Ntgroup; i++) { + N = MAX (N, tgroup[i].Nmeasure); + } + + ALLOCATE (xlist, double, N); + ALLOCATE (ylist, double, N); + + for (i = 0; i < Ntgroup; i++) { + N = 0; + for (j = 0; j < tgroup[i].Nmeasure; j++) { + + m = tgroup[i].measure[j]; + c = tgroup[i].catalog[j]; + + if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue; + + // ave = catalog[c].measureT[m].averef; + xlist[N] = catalog[c].measureT[m].R; + ylist[N] = catalog[c].measureT[m].D; + N++; + } + + sprintf (string, "TGroup "OFF_T_FMT, i); + plot_defaults (&graphdata); + plot_list (&graphdata, xlist, ylist, N, string, NULL); + } + + free (ylist); + free (xlist); +} + +void plot_tgroups () { + + off_t i, bin; + double *xlist, *Mlist, *dlist; + Graphdata graphdata; + + if (!TGROUP_ZEROPT) return; + + ALLOCATE (xlist, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + ALLOCATE (Mlist, double, Ntgroup); + + for (i = 0; i < Ntgroup; i++) { + Mlist[i] = tgroup[i].McalPSF; + dlist[i] = tgroup[i].dMcal; + xlist[i] = tgroup[i].start; + } + + plot_defaults (&graphdata); + graphdata.xmin = 0.95; + graphdata.xmax = 2.50; + graphdata.ymin = PlotdMmin; + graphdata.ymax = PlotdMmax; + plot_list (&graphdata, xlist, Mlist, Ntgroup, "airmass vs Mcal", "%s.airmass.png", OUTROOT); + plot_defaults (&graphdata); + graphdata.size = 1.5; + graphdata.ptype = 7; + plot_list (&graphdata, Mlist, dlist, Ntgroup, "Mcal vs dMcal", "%s.MdM.png", OUTROOT); + +# define NBIN 200 + REALLOCATE (xlist, double, NBIN); + REALLOCATE (Mlist, double, NBIN); + + /**** dMcal histgram ****/ + for (i = 0; i < NBIN; i++) xlist[i] = 0.00005*i; + bzero (Mlist, NBIN*sizeof(double)); + for (i = 0; i < Ntgroup; i++) { + bin = tgroup[i].dMcal / 0.00005; + bin = MAX (0, MIN (NBIN - 1, bin)); + Mlist[bin] += 1.0; + } + plot_defaults (&graphdata); + graphdata.style = 1; + plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "%s.dMcalhist.png", OUTROOT); + + free (dlist); + free (xlist); + free (Mlist); +} + Index: /branches/eam_branches/relphot.20200519/src/args.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/args.c (revision 41361) +++ /branches/eam_branches/relphot.20200519/src/args.c (revision 41362) @@ -315,4 +315,13 @@ exit (2); } + } + + // calibrate the photometry by night / time period + TGROUP_ZEROPT = FALSE; + if ((N = get_argument (argc, argv, "-tgroup"))) { + TGROUP_ZEROPT = TRUE; + remove_argument (N, &argc, argv); + loadTGroups (argv[N]); + remove_argument (N, &argc, argv); } Index: /branches/eam_branches/relphot.20200519/src/load_images.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/load_images.c (revision 41361) +++ /branches/eam_branches/relphot.20200519/src/load_images.c (revision 41362) @@ -110,4 +110,6 @@ MARKTIME("init images: %f sec\n", dtime); + initTGroups (subset, LineNumber, Nsubset); + return TRUE; } Index: /branches/eam_branches/relphot.20200519/src/relphot_images.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/relphot_images.c (revision 41361) +++ /branches/eam_branches/relphot.20200519/src/relphot_images.c (revision 41362) @@ -70,4 +70,6 @@ initMosaicBins (catalog, Ncatalog, TRUE); + initTGroupBins (catalog, Ncatalog); + initGridBins (catalog, Ncatalog); initMrel (catalog, Ncatalog); @@ -77,4 +79,7 @@ findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ + MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); + + findTGroups (catalog, Ncatalog); /* also sets Grid values */ MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); @@ -129,4 +134,6 @@ setMmos (catalog, FALSE, flatcorr); + setMgrp (catalog, FALSE, flatcorr); + setMgrid (catalog, flatcorr); MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); @@ -172,7 +179,9 @@ setMcal (catalog, onlyPoorImages, flatcorr); setMmos (catalog, onlyPoorImages, flatcorr); + setMgrp (catalog, onlyPoorImages, flatcorr); MARKTIME("-- finalize Mcal values: %f sec\n", dtime); - setMcalFinal (); // copy per-mosaic calibrations to the images + setMcalFromMosaics (); // copy per-mosaic calibrations to the images + setMcalFromTGroups (); // copy per-tgroup calibrations to the images if (SAVE_IMAGE_UPDATES) { @@ -255,5 +264,6 @@ // If nloop == 0, the above pass is not performed, in which case // the Mcal values passed to the mosaics are not returned to the images... - setMcalFinal (); // copy per-mosaic calibrations to the images + setMcalFromMosaics (); // copy per-mosaic calibrations to the images + setMcalFromTGroups (); // copy per-tgroup calibrations to the images } Index: /branches/eam_branches/relphot.20200519/src/relphot_parallel_images.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/relphot_parallel_images.c (revision 41361) +++ /branches/eam_branches/relphot.20200519/src/relphot_parallel_images.c (revision 41362) @@ -116,4 +116,5 @@ setMcal (catalog, FALSE, flatcorr); setMmos (catalog, FALSE, flatcorr); + // setMgrp (catalog, FALSE, flatcorr); XXX think this through: this may not make sense for tgroups MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); @@ -152,7 +153,9 @@ setMcal (catalog, onlyPoorImages, flatcorr); setMmos (catalog, onlyPoorImages, flatcorr); + // setMgrp (catalog, onlyPoorImages, flatcorr); XXX see above MARKTIME("-- finalize Mcal values: %f sec\n", dtime); - setMcalFinal (); // copy per-mosaic calibrations to the images + setMcalFromMosaics (); // copy per-mosaic calibrations to the images + // XXX keep? setMcalFromTGroups (); // copy per-tgroup calibrations to the images share_image_mags (regionHosts, -1);