Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/include/relphot.h (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/include/relphot.h (revision 36976) @@ -48,5 +48,5 @@ float dMsys; unsigned short nFitPhotom; - short Xm; + float Mchisq; float secz; float ubercalDist; @@ -119,5 +119,5 @@ float M; float dM; - float Xm; + float Mchisq; int Nsec; unsigned int objID; Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/BrightCatalog.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/BrightCatalog.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/BrightCatalog.c (revision 36976) @@ -182,12 +182,12 @@ // need to create and assign to flat-field correction - GET_COLUMN(M, "MAG", float); - GET_COLUMN(dM, "MAG_ERR", float); - GET_COLUMN(Xm, "MAG_CHI", float); - GET_COLUMN(flags, "FLAGS", int); - GET_COLUMN(Ncode, "NCODE", short); - GET_COLUMN(Nused, "NUSED", short); - GET_COLUMN(Mmin, "MAG_MIN", float); - GET_COLUMN(Mmax, "MAG_MAX", float); + GET_COLUMN(M, "MAG", float); + GET_COLUMN(dM, "MAG_ERR", float); + GET_COLUMN(Mchisq, "MAG_CHI", float); + GET_COLUMN(flags, "FLAGS", int); + GET_COLUMN(Ncode, "NCODE", short); + GET_COLUMN(Nused, "NUSED", short); + GET_COLUMN(Mmin, "MAG_MIN", float); + GET_COLUMN(Mmax, "MAG_MAX", float); gfits_free_header (&theader); gfits_free_table (&ftable); @@ -196,23 +196,23 @@ ALLOCATE (secfilt, SecFilt, Nrow); for (i = 0; i < Nrow; i++) { - secfilt[i].M = M[i]; - secfilt[i].dM = dM[i]; - secfilt[i].Xm = Xm[i]; - secfilt[i].flags = flags[i]; - secfilt[i].Ncode = Ncode[i]; - secfilt[i].Nused = Nused[i]; - secfilt[i].Mmin = Mmin[i]; - secfilt[i].Mmax = Mmax[i]; + secfilt[i].M = M[i]; + secfilt[i].dM = dM[i]; + secfilt[i].Mchisq = Mchisq[i]; + secfilt[i].flags = flags[i]; + secfilt[i].Ncode = Ncode[i]; + secfilt[i].Nused = Nused[i]; + secfilt[i].Mmin = Mmin[i]; + secfilt[i].Mmax = Mmax[i]; } fprintf (stderr, "loaded data for %lld secfilt\n", (long long) Nrow); - free (M ); - free (dM ); - free (Xm ); - free (flags); - free (Ncode); - free (Nused); - free (Mmin ); - free (Mmax ); + free (M ); + free (dM ); + free (Mchisq); + free (flags ); + free (Ncode ); + free (Nused ); + free (Mmin ); + free (Mmax ); catalog->secfilt = secfilt; // assert Nsecfilt * Naverage = Nrow? @@ -447,5 +447,5 @@ float *M ; ALLOCATE (M , float, Nsec); float *dM ; ALLOCATE (dM , float, Nsec); - float *Xm ; ALLOCATE (Xm , float, Nsec); + float *Mchisq ; ALLOCATE (Mchisq , float, Nsec); int *flags ; ALLOCATE (flags , int, Nsec); short *Ncode ; ALLOCATE (Ncode , short, Nsec); @@ -457,27 +457,27 @@ SecFilt *secfilt = catalog->secfilt; for (i = 0; i < Nsec; i++) { - M [i] = secfilt[i]. M ; - dM [i] = secfilt[i]. dM ; - Xm [i] = secfilt[i]. Xm ; - flags [i] = secfilt[i]. flags ; - Ncode [i] = secfilt[i]. Ncode ; - Nused [i] = secfilt[i]. Nused ; - Mmin [i] = secfilt[i]. Mmin ; - Mmax [i] = secfilt[i]. Mmax ; + M [i] = secfilt[i].M ; + dM [i] = secfilt[i].dM ; + Mchisq[i] = secfilt[i].Mchisq ; + flags [i] = secfilt[i].flags ; + Ncode [i] = secfilt[i].Ncode ; + Nused [i] = secfilt[i].Nused ; + Mmin [i] = secfilt[i].Mmin ; + Mmax [i] = secfilt[i].Mmax ; } // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "MAG", M , Nsec); - gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM , Nsec); - gfits_set_bintable_column (&theader, &ftable, "MAG_CHI", Xm , Nsec); - gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, Nsec); - gfits_set_bintable_column (&theader, &ftable, "NCODE", Ncode, Nsec); - gfits_set_bintable_column (&theader, &ftable, "NUSED", Nused, Nsec); - gfits_set_bintable_column (&theader, &ftable, "MAG_MIN", Mmin , Nsec); - gfits_set_bintable_column (&theader, &ftable, "MAG_MAX", Mmax , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG", M , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_CHI", Mchisq, Nsec); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nsec); + gfits_set_bintable_column (&theader, &ftable, "NCODE", Ncode , Nsec); + gfits_set_bintable_column (&theader, &ftable, "NUSED", Nused , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_MIN", Mmin , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_MAX", Mmax , Nsec); free (M ); free (dM ); - free (Xm ); + free (Mchisq ); free (flags ); free (Ncode ); Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/ImageOps.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/ImageOps.c (revision 36976) @@ -701,5 +701,5 @@ double MaxOffset, MaxScatter, MedOffset; - if (FREEZE_IMAGES) return; + // if (FREEZE_IMAGES) return; if (VERBOSE) fprintf (stderr, "marking poor images\n"); @@ -712,4 +712,8 @@ for (i = N = 0; i < Nimage; i++) { if (image[i].flags & IMAGE_BAD) continue; + + int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500)); + if (FREEZE_IMAGES && isGPC1cam) continue; + mlist[N] = image[i].Mcal; slist[N] = image[i].dMcal; @@ -803,5 +807,5 @@ Graphdata graphdata; - if (FREEZE_IMAGES) return; + // if (FREEZE_IMAGES) return; ALLOCATE (xlist, double, Nimage); @@ -816,8 +820,15 @@ float mindMcal = +100.0; float maxdMcal = -100.0; - for (i = 0; i < Nimage; i++) { - Mlist[i] = image[i].Mcal; - dlist[i] = image[i].dMcal; - xlist[i] = image[i].secz; + + int Nplot = 0; + + for (i = 0; i < Nimage; i++) { + + int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500)); + if (FREEZE_IMAGES && isGPC1cam) continue; + + Mlist[Nplot] = image[i].Mcal; + dlist[Nplot] = image[i].dMcal; + xlist[Nplot] = image[i].secz; minAirmass = MIN (image[i].secz, minAirmass); maxAirmass = MAX (image[i].secz, maxAirmass); @@ -826,4 +837,6 @@ mindMcal = MIN (image[i].dMcal, mindMcal); maxdMcal = MAX (image[i].dMcal, maxdMcal); + + Nplot ++; } @@ -858,5 +871,8 @@ for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i; bzero (Mlist, NBIN*sizeof(double)); - for (i = 0; i < Nimage; i++) { + + for (i = 0; i < Nimage; i++) { + int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500)); + if (FREEZE_IMAGES && isGPC1cam) continue; bin = image[i].dMcal / 0.00025; bin = MAX (0, MIN (NBIN - 1, bin)); @@ -883,5 +899,7 @@ StatType stats; bzero (&stats, sizeof (StatType)); - if (FREEZE_IMAGES) return (stats); + + // we no longer blindly apply FREEZE_IMAGES to all images, only to mosaics + // if (FREEZE_IMAGES) return (stats); ALLOCATE (list, double, Nimage); @@ -891,4 +909,7 @@ for (i = 0; i < Nimage; i++) { if (image[i].flags & IMAGE_BAD) continue; + + int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500)); + if (FREEZE_IMAGES && isGPC1cam) continue; N = 0; @@ -926,5 +947,5 @@ bzero (&stats, sizeof (StatType)); - if (FREEZE_IMAGES) return (stats); + // if (FREEZE_IMAGES) return (stats); ALLOCATE (list, double, Nimage); @@ -935,4 +956,7 @@ if (image[i].flags & IMAGE_BAD) continue; + + int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500)); + if (FREEZE_IMAGES && isGPC1cam) continue; list[n] = pow (10.0, 0.01*image[i].Xm); @@ -956,5 +980,5 @@ bzero (&stats, sizeof (StatType)); - if (FREEZE_IMAGES) return (stats); + // if (FREEZE_IMAGES) return (stats); ALLOCATE (list, double, Nimage); @@ -965,4 +989,7 @@ if (image[i].flags & IMAGE_BAD) continue; + + int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500)); + if (FREEZE_IMAGES && isGPC1cam) continue; list[n] = image[i].Mcal; @@ -986,5 +1013,5 @@ bzero (&stats, sizeof (StatType)); - if (FREEZE_IMAGES) return (stats); + // if (FREEZE_IMAGES) return (stats); ALLOCATE (list, double, Nimage); @@ -995,4 +1022,7 @@ if (image[i].flags & IMAGE_BAD) continue; + + int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500)); + if (FREEZE_IMAGES && isGPC1cam) continue; list[n] = image[i].dMcal; Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/MeanMagIO.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/MeanMagIO.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/MeanMagIO.c (revision 36976) @@ -60,5 +60,5 @@ GET_COLUMN (M, "MAG", float); GET_COLUMN (dM, "MAG_ERR", float); - GET_COLUMN (Xm, "MAG_CHISQ", float); + GET_COLUMN (Mchisq, "MAG_CHISQ", float); GET_COLUMN (Nsec, "NSEC", int); GET_COLUMN (objID, "OBJ_ID", int); @@ -73,19 +73,19 @@ ALLOCATE (meanmags, MeanMag, Nrow); for (i = 0; i < Nrow; i++) { - meanmags[i].M = M [i]; - meanmags[i].dM = dM [i]; - meanmags[i].Xm = Xm [i]; - meanmags[i].Nsec = Nsec [i]; - meanmags[i].objID = objID[i]; - meanmags[i].catID = catID[i]; + meanmags[i].M = M [i]; + meanmags[i].dM = dM [i]; + meanmags[i].Mchisq = Mchisq [i]; + meanmags[i].Nsec = Nsec [i]; + meanmags[i].objID = objID [i]; + meanmags[i].catID = catID [i]; } fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow); - free (M ); - free (dM ); - free (Xm ); - free (Nsec ); - free (objID); - free (catID); + free (M ); + free (dM ); + free (Mchisq); + free (Nsec ); + free (objID ); + free (catID ); *nmeanmags = Nrow; @@ -125,5 +125,5 @@ gfits_create_table (&theader, &ftable); - float *M, *dM, *Xm; + float *M, *dM, *Mchisq; int *Nsec; unsigned int *objID, *catID; @@ -132,5 +132,5 @@ ALLOCATE (M, float, Nmeanmags); ALLOCATE (dM, float, Nmeanmags); - ALLOCATE (Xm, float, Nmeanmags); + ALLOCATE (Mchisq, float, Nmeanmags); ALLOCATE (Nsec, int, Nmeanmags); ALLOCATE (objID, unsigned int, Nmeanmags); @@ -139,10 +139,10 @@ // assign the storage arrays for (i = 0; i < Nmeanmags; i++) { - M [i] = meanmags[i].M ; - dM [i] = meanmags[i].dM ; - Xm [i] = meanmags[i].Xm ; - Nsec [i] = meanmags[i].Nsec ; - objID[i] = meanmags[i].objID; - catID[i] = meanmags[i].catID; + M [i] = meanmags[i].M ; + dM [i] = meanmags[i].dM ; + Mchisq[i] = meanmags[i].Mchisq; + Nsec [i] = meanmags[i].Nsec ; + objID [i] = meanmags[i].objID ; + catID [i] = meanmags[i].catID ; } @@ -150,15 +150,15 @@ gfits_set_bintable_column (&theader, &ftable, "MAG", M, Nmeanmags); gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, Nmeanmags); - gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Xm, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Mchisq, Nmeanmags); gfits_set_bintable_column (&theader, &ftable, "NSEC", Nsec, Nmeanmags); gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, Nmeanmags); gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, Nmeanmags); - free (M ); - free (dM ); - free (Xm ); - free (Nsec ); - free (objID); - free (catID); + free (M ); + free (dM ); + free (Mchisq); + free (Nsec ); + free (objID ); + free (catID ); FILE *f = fopen (filename, "w"); Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/MosaicOps.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/MosaicOps.c (revision 36976) @@ -96,12 +96,12 @@ /* a new mosaic, define ranges -- preserve the original values incase this image is not used */ - mosaic[Nmosaic].start = start; - mosaic[Nmosaic].stop = stop; - mosaic[Nmosaic].Mcal = 0.0; // note : mosaic stores only offsets relative to the original image values - mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images - mosaic[Nmosaic].dMsys = 0.0; - mosaic[Nmosaic].Xm = 0.0; + mosaic[Nmosaic].start = start; + mosaic[Nmosaic].stop = stop; + mosaic[Nmosaic].Mcal = 0.0; // note : mosaic stores only offsets relative to the original image values + mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images + mosaic[Nmosaic].dMsys = 0.0; + mosaic[Nmosaic].Mchisq = 0.0; mosaic[Nmosaic].flags = image[i].flags; - mosaic[Nmosaic].secz = image[i].secz; + mosaic[Nmosaic].secz = image[i].secz; mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); @@ -172,8 +172,14 @@ // generate a list of all subset image start times ALLOCATE (startTimes, unsigned int, Nsubset); + int Nmoschip = 0; for (i = 0; i < Nsubset; i++) { - startTimes[i] = subset[i].tzero; - } - sort_times (startTimes, Nsubset); + int isGPC1cam = ((subset[i].photcode >= 10000) && (subset[i].photcode <= 10500)); + if (!isGPC1cam) { + continue; + } + startTimes[Nmoschip] = subset[i].tzero; + Nmoschip ++; + } + sort_times (startTimes, Nmoschip); MARKTIME("create array of all image obstimes: %f sec\n", dtime); @@ -184,5 +190,5 @@ // generate a list of the unique start times (these define the mosaics) - for (i = 0; i < Nsubset; i++) { + for (i = 0; i < Nmoschip; i++) { if (startTimes[i] < startTimesMosaic[Nmosaic]) { fprintf (stderr, "error?\n"); @@ -210,14 +216,14 @@ for (i = 0; i < Nmosaic; i++) { /* a new mosaic, define ranges */ - mosaic[i].start = startTimesMosaic[i]; - mosaic[i].stop = 0; - mosaic[i].Mcal = 0.0; - mosaic[i].dMcal = 0.0; - mosaic[i].dMsys = 0.0; - mosaic[i].Xm = 0.0; - mosaic[i].flags = 0; - mosaic[i].secz = NAN; + mosaic[i].start = startTimesMosaic[i]; + mosaic[i].stop = 0; + mosaic[i].Mcal = 0.0; + mosaic[i].dMcal = 0.0; + mosaic[i].dMsys = 0.0; + mosaic[i].Mchisq = 0.0; + mosaic[i].flags = 0; + mosaic[i].secz = NAN; mosaic[i].photcode = 0; - mosaic[i].skipCal = FALSE; + mosaic[i].skipCal = FALSE; memset (&mosaic[i].coords, 0, sizeof(Coords)); @@ -238,9 +244,18 @@ } - /* select valid mosaic images by photcode */ - pname = GetPhotcodeNamebyCode (image[i].photcode); - if (!pname) continue; - status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); - if (status) continue; + // this is gpc1-specific. needs a more general solution + int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500)); + if (!isGPC1cam) { + continue; + } + + if (0) { + // XXX this version fails if we have stack and fforce data: GPC1.g.SkyChip + /* select valid mosaic images by photcode */ + pname = GetPhotcodeNamebyCode (image[i].photcode); + if (!pname) continue; + status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); + if (status) continue; + } /* set image time range */ @@ -262,11 +277,22 @@ // assign each image to a mosaic + int Nsimple = 0; for (i = 0; i < Nsubset; i++) { ImageToMosaic[i] = -1; - /* select valid mosaic images by photcode */ - pname = GetPhotcodeNamebyCode (subset[i].photcode); - status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); - if (status) continue; + // this is gpc1-specific. needs a more general solution + int isGPC1cam = ((subset[i].photcode >= 10000) && (subset[i].photcode <= 10500)); + if (!isGPC1cam) { + Nsimple ++; + continue; + } + + if (0) { + // XXX this version fails if we have stack and fforce data: GPC1.g.SkyChip + /* select valid mosaic images by photcode */ + pname = GetPhotcodeNamebyCode (subset[i].photcode); + status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); + if (status) continue; + } start = subset[i].tzero; @@ -302,5 +328,5 @@ mosaic[j].Mcal = 0.0; mosaic[j].dMcal = 0.0; - mosaic[j].Xm = 0.0; + mosaic[j].Mchisq = 0.0; mosaic[j].dMsys = subset[i].flags; mosaic[j].flags = subset[i].flags; @@ -316,5 +342,5 @@ initMosaicGrid (subset, Nsubset); - fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nsubset, (int) Nmosaic); + fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple); return; } @@ -384,5 +410,5 @@ mosaic[i].dMcal = 0.0; mosaic[i].dMsys = 0.0; - mosaic[i].Xm = 0.0; + mosaic[i].Mchisq = 0.0; mosaic[i].flags = 0; mosaic[i].secz = NAN; @@ -443,5 +469,5 @@ mosaic[j].Mcal = 0.0; mosaic[j].dMcal = 0.0; - mosaic[j].Xm = 0.0; + mosaic[j].Mchisq = 0.0; mosaic[j].dMsys = image[i].flags; mosaic[j].flags = image[i].flags; @@ -566,5 +592,5 @@ Mcal += image[m].Mcal; dMcal += image[m].dMcal; - Xm += image[m].Xm; + Xm += image[m].Mchisq; // for ubercal images, we (elsewhere) keep Mcal frozen @@ -573,5 +599,5 @@ image[m].Mcal = 0.0; image[m].dMcal = NAN; - image[m].Xm = NAN_S_SHORT; + image[m].Mchisq = NAN_S_SHORT; } @@ -617,7 +643,7 @@ mosaic[i].coords.pc2_1 = 0.0; - mosaic[i].Mcal = Mcal / MosaicN_Image[i]; - mosaic[i].dMcal = dMcal / MosaicN_Image[i]; - mosaic[i].Xm = Xm / MosaicN_Image[i]; + mosaic[i].Mcal = Mcal / MosaicN_Image[i]; + mosaic[i].dMcal = dMcal / MosaicN_Image[i]; + mosaic[i].Mchisq = Xm / MosaicN_Image[i]; } return; @@ -669,5 +695,5 @@ Mcal += image[m].Mcal; dMcal += image[m].dMcal; - Xm += image[m].Xm; + Xm += image[m].Mchisq; // XXX: how does this work with UBERCAL? We want to keep the Mcal values supplied by ubercal, but @@ -679,7 +705,7 @@ /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ - image[m].Mcal = 0.0; - image[m].dMcal = NAN; - image[m].Xm = NAN_S_SHORT; + image[m].Mcal = 0.0; + image[m].dMcal = NAN; + image[m].Mchiqs = NAN; } dS /= MosaicN_Image[i]; @@ -697,7 +723,7 @@ RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords); - mosaic[i].Mcal = Mcal / MosaicN_Image[i]; - mosaic[i].dMcal = dMcal / MosaicN_Image[i]; - mosaic[i].Xm = Xm / MosaicN_Image[i]; + mosaic[i].Mcal = Mcal / MosaicN_Image[i]; + mosaic[i].dMcal = dMcal / MosaicN_Image[i]; + mosaic[i].Mchisq = Xm / MosaicN_Image[i]; } if (!USE_GRID) return; @@ -729,5 +755,5 @@ image[im].Mcal += mosaic[i].Mcal; image[im].dMcal = mosaic[i].dMcal; - image[im].Xm = mosaic[i].Xm; + image[im].Mchisq = mosaic[i].Mchisq; image[im].dMagSys = mosaic[i].dMsys; image[im].nFitPhotom = mosaic[i].nFitPhotom; @@ -840,6 +866,7 @@ mosID = ImageToMosaic[idx]; if (mosID < 0) { - Image *image = getimage(idx); - fprintf (stderr, "unmatched image %s\n", image[0].name); + // Image *image = getimage(idx); + // fprintf (stderr, "unmatched image %s\n", image[0].name); + // skip measurements from simple chips (not matched to a mosaic by definition) return; } @@ -882,6 +909,7 @@ if (!MOSAIC_ZEROPT) return (0); + // unassigned measurements belong to simple chips i = MeasureToMosaic[cat][meas]; - if (i == -1) return (NAN); + if (i == -1) return (0.0); if (mosaic[i].flags & IMAGE_BAD) return (NAN); @@ -1219,5 +1247,5 @@ myMosaic[0].dMcal = stats.error; myMosaic[0].nFitPhotom = N; - myMosaic[0].Xm = 100.0*log10(stats.chisq); + myMosaic[0].Mchisq = stats.chisq; if (testImage) { @@ -1239,5 +1267,5 @@ if (testImage) { - fprintf (stderr, "%f %f : %d %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, myMosaic[0].Xm, pow(10.0, 0.01*myMosaic[0].Xm)); + fprintf (stderr, "%f %f : %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, myMosaic[0].Mchisq); } @@ -1699,5 +1727,5 @@ if (mosaic[i].flags & IMAGE_BAD) continue; if (mosaic[i].skipCal) continue; - list[n] = pow(10.0, 0.01*mosaic[i].Xm); + list[n] = mosaic[i].Mchisq; dlist[n] = 1; n++; @@ -1784,5 +1812,4 @@ off_t i, j, m, c, N, Nimage; double *xlist, *ylist; - // double Xmin, Xmax, Ymin, Ymax; char string[64]; Graphdata graphdata; @@ -1803,6 +1830,4 @@ for (i = 0; i < Nmosaic; i++) { N = 0; - // Xmin = Ymin = +360.0; - // Xmax = Ymax = -360.0; for (j = 0; j < N_onMosaic[i]; j++) { Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/StarOps.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/StarOps.c (revision 36976) @@ -153,6 +153,6 @@ int i; - fprintf (stderr, "this version of setMrel is invalid for now\n"); - exit (1); + // fprintf (stderr, "this version of setMrel is invalid for now\n"); + // exit (1); int Nsecfilt = GetPhotcodeNsecfilt (); @@ -673,7 +673,7 @@ liststats (list, dlist, wlist, N, &stats); - catalog[Nc].secfilt[Nsecfilt*j+Nsec].M = stats.mean; - catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error; - catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].M = stats.mean; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mchisq = (stats.Nmeas > 1) ? stats.chisq : NAN; // when running -averages, we have no information about the images, so we cannot set this @@ -953,5 +953,5 @@ int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan; - float dM, Xm; + float dM; double Chisq, MaxScatter, MaxChisq; double *xlist, *slist, *dlist; @@ -981,8 +981,6 @@ for (j = 0; j < catalog[i].Naverage; j++) { if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue; - Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm; - if (Xm == -1) continue; - Chisq = pow (10.0, 0.01*Xm); - xlist[Ntot] = Chisq; + Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; + xlist[Ntot] = Mchisq; slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; dlist[Ntot] = 1; @@ -1002,13 +1000,12 @@ for (j = 0; j < catalog[i].Naverage; j++) { dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; - Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm; - Chisq = pow (10.0, 0.01*Xm); - mark = (dM > MaxScatter) || (Xm == NAN_S_SHORT) || (Chisq > MaxChisq); + Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; + mark = (dM > MaxScatter) || (isnan(Mchisq) || (Mchisq > MaxChisq); if (mark) { catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR; Ndel ++; - if (dM > MaxScatter) { Nscat ++; } - if (Xm == NAN_S_SHORT) { Nnan ++; } - if (Chisq > MaxChisq) { Nchi ++; } + if (dM > MaxScatter) { Nscat ++; } + if (isnan(Mchisq)) { Nnan ++; } + if (Mchisq > MaxChisq) { Nchi ++; } } else { catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR; @@ -1236,5 +1233,5 @@ off_t j, Ntot; - int i, n, Xm; + int i, n; double *list, *dlist; StatType stats; @@ -1257,7 +1254,7 @@ if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm; - if (Xm == NAN_S_SHORT) continue; - list[n] = pow (10.0, 0.01*Xm); + Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; + if (isnan(Mchisq)) continue; + list[n] = Mchisq; dlist[n] = 1; n++; @@ -1380,7 +1377,7 @@ if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].M; - value = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm; - if (value == NAN_S_SHORT) continue; - ylist[N] = 0.01*value; + value = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; + if (isnan(value)) continue; + ylist[N] = value; N++; } Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/args.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/args.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/args.c (revision 36976) @@ -227,8 +227,8 @@ } - SHOW_PARAMS = FALSE; + SHOW_PARAMS = TRUE; if ((N = get_argument (argc, argv, "-params"))) { remove_argument (N, &argc, argv); - SHOW_PARAMS = TRUE; + SHOW_PARAMS = FALSE; } Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/assign_images.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/assign_images.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/assign_images.c (revision 36976) @@ -27,5 +27,5 @@ makeMosaics (image, Nimage); - // center coords and Mcal, dMcal, Xm for the mosaics + // center coords and Mcal, dMcal, Mchisq for the mosaics setMosaicCenters (image, Nimage); MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime); Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/bcatalog.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/bcatalog.c (revision 36976) @@ -22,4 +22,7 @@ ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); Nmeasure = Naverage = 0; + + int NmMax = 0; + int NmMin = catalog[0].Nmeasure; Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0; @@ -82,5 +85,5 @@ // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN; // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Xm = NAN; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mchisq = NAN; // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20 = NAN_S_SHORT; // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80 = NAN_S_SHORT; @@ -187,5 +190,5 @@ } myAssert (Nmeasure < NMEASURE, "realloc failure"); - } + } // end of catalog.average.Nmeasure loop // skip object if it is likely to be a galaxy @@ -212,4 +215,6 @@ continue; } + NmMin = MIN (Nm, NmMin); + NmMax = MAX (Nm, NmMax); // for w-band photometry (& other cases?) convert gri photometry @@ -242,4 +247,11 @@ subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt; + if (VERBOSE) { + fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" ("OFF_T_FMT" measures) for catalog %s, %d < Nm < %d\n", + subcatalog[0].Naverage, subcatalog[0].Nmeasure, catalog[0].Naverage, catalog[0].Nmeasure, catalog[0].filename, NmMin, NmMax); + fprintf (stderr, "rejections: %d stars have too few measures:\n %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", + Nfew, Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor); + } + // limit the total number of stars in the catalog if (MaxDensityUse) { @@ -247,10 +259,4 @@ } - if (VERBOSE) { - fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", - subcatalog[0].Naverage, subcatalog[0].Nmeasure, i, catalog[0].filename); - fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", - Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor); - } return (TRUE); } Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/initialize.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/initialize.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/initialize.c (revision 36976) @@ -36,5 +36,14 @@ if (SHOW_PARAMS) { - fprintf (stderr, "current parameter settings:\n"); + int Ns; + fprintf (stderr, "subset selection criteria:\n"); + fprintf (stderr, " photcodes "); + for (Ns = 0; Ns < Nphotcodes; Ns++) { + if (Ns == Nphotcodes - 1) { + fprintf (stderr, "%s\n", photcodes[Ns][0].name); + } else { + fprintf (stderr, "%s, ", photcodes[Ns][0].name); + } + } if (TimeSelect) { fprintf (stderr, "TimeSelect: TRUE (%s - %s)\n", ohana_sec_to_date (TSTART), ohana_sec_to_date (TSTOP)); @@ -42,4 +51,18 @@ fprintf (stderr, "TimeSelect: FALSE\n"); } + if (DophotSelect) { + fprintf (stderr, "DophotSelect: TRUE (%d)\n", DophotValue); + } else { + fprintf (stderr, "DophotSelect: FALSE\n"); + } + fprintf (stderr, "PSF_QF limit: 0.85 (hardwired)\n"); + + // fprintf (stderr, "Photom Bad Mask: 0x%08x, Photom Poor Mask: 0x%08x\n"); + + fprintf (stderr, "MAG_LIM: %f, SIGMA_LIM: %f\n", MAG_LIM, SIGMA_LIM); + fprintf (stderr, "INST_MAG_MIN: %f, INST_MAG_MAX: %f\n", ImagMin, ImagMax); + + fprintf (stderr, "STAR_TOOFEW: %d\n", STAR_TOOFEW); + fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF); fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n", @@ -49,5 +72,4 @@ (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING)); - fprintf (stderr, "MAG_LIM %lf\n", MAG_LIM); fprintf (stderr, "STAR_SCATTER %lf\n", STAR_SCATTER); fprintf (stderr, "IMAGE_SCATTER %lf\n", IMAGE_SCATTER); @@ -56,5 +78,4 @@ fprintf (stderr, "GSCFILE %s\n", GSCFILE); fprintf (stderr, "CATDIR %s\n", CATDIR); - exit (0); } Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelCatalog.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelCatalog.c (revision 36976) @@ -354,7 +354,7 @@ liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats); - secfilt[Nsec].M = psfstats->mean; - secfilt[Nsec].dM = psfstats->error; - secfilt[Nsec].Xm = (psfstats->Nmeas > 1) ? 100.0*log10(psfstats->chisq + 1e-4) : NAN_S_SHORT; + secfilt[Nsec].M = psfstats->mean; + secfilt[Nsec].dM = psfstats->error; + secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; // when running -averages, we have no information about the images, so we cannot set this @@ -613,7 +613,7 @@ secfilt[Nsec].dFapStk = zpFactor * measure[k].dFluxAp; - secfilt[Nsec].MpsfStk = (measure[k].FluxPSF > 0.0) ? zp -2.5*log10(measure[k].FluxPSF) : NAN; - secfilt[Nsec].MkronStk = (measure[k].FluxKron > 0.0) ? zp -2.5*log10(measure[k].FluxKron) : NAN; - secfilt[Nsec].MapStk = (measure[k].FluxAp > 0.0) ? zp -2.5*log10(measure[k].FluxAp) : NAN; + secfilt[Nsec].MpsfStk = (measure[k].FluxPSF > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN; + secfilt[Nsec].MkronStk = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN; + secfilt[Nsec].MapStk = (measure[k].FluxAp > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN; secfilt[Nsec].stackDetectID = ID; @@ -716,5 +716,5 @@ if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); } - float Fcal = MagToFlux(Mcal); + float Fcal = MagToFlux(-Mcal); // in the calculations below, @@ -767,5 +767,5 @@ secfilt[Nsec].sFpsfWrp = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision) secfilt[Nsec].NusedWrp = psfstats->Nmeas; - secfilt[Nsec].MpsfWrp = isnan(secfilt[Nsec].FpsfWrp) ? NAN : -2.5*log(secfilt[Nsec].FpsfWrp); // XXX need to add zero-point + secfilt[Nsec].MpsfWrp = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy // NOTE : use the modified weight for apmags as well as psf mags @@ -775,5 +775,5 @@ secfilt[Nsec].sFapWrp = Nap > 0 ? apstats->sigma : NAN; secfilt[Nsec].NusedApWrp = Nap; - secfilt[Nsec].MapWrp = isnan(secfilt[Nsec].FapWrp) ? NAN : -2.5*log(secfilt[Nsec].FapWrp); // XXX need to add zero-point + secfilt[Nsec].MapWrp = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats); @@ -782,5 +782,5 @@ secfilt[Nsec].sFkronWrp = Nkron > 0 ? kronstats->sigma : NAN; secfilt[Nsec].NusedKronWrp = Nkron; - secfilt[Nsec].MkronWrp = isnan(secfilt[Nsec].FkronWrp) ? NAN : -2.5*log(secfilt[Nsec].FkronWrp); // XXX need to add zero-point + secfilt[Nsec].MkronWrp = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy } return (TRUE); Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/share_mean_mags.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/share_mean_mags.c (revision 36975) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/share_mean_mags.c (revision 36976) @@ -123,5 +123,5 @@ meanmags->M = secfilt->M; meanmags->dM = secfilt->dM; - meanmags->Xm = secfilt->Xm; + meanmags->Mchisq = secfilt->Mchisq; meanmags->Nsec = Nsec; // key to secfilt entry