Index: /branches/eam_branches/ipp-20140610/Ohana/src/libautocode/def/secfilt-ps1-v5.d =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/libautocode/def/secfilt-ps1-v5.d (revision 36916) +++ /branches/eam_branches/ipp-20140610/Ohana/src/libautocode/def/secfilt-ps1-v5.d (revision 36917) @@ -2,5 +2,5 @@ EXTNAME DVO_SECFILT_PS1_V5 TYPE BINTABLE -SIZE 152 +SIZE 160 DESCRIPTION DVO SecFilt : Secondary Filter Data @@ -10,6 +10,8 @@ FIELD Map, MAG_AP, float, ave aperture mag in this band, mags FIELD dMap, MAG_AP_ERR, float, ave aperture mag in this band, mags +FIELD sMap, MAG_AP_STDEV, float, standard deviation of ap mags, mags FIELD Mkron, MAG_KRON, float, ave kron mag in this band, mags FIELD dMkron, MAG_KRON_ERR, float, formal error on average kron mag, mags +FIELD sMkron, MAG_KRON_STDEV, float, standard deviation of kron mags, mags # these statistics are PSF-specific Index: /branches/eam_branches/ipp-20140610/Ohana/src/libautocode/def/secfilt.d =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/libautocode/def/secfilt.d (revision 36916) +++ /branches/eam_branches/ipp-20140610/Ohana/src/libautocode/def/secfilt.d (revision 36917) @@ -2,5 +2,5 @@ EXTNAME DVO_SECFILT TYPE BINTABLE -SIZE 152 +SIZE 160 DESCRIPTION DVO SecFilt : Secondary Filter Data @@ -10,6 +10,8 @@ FIELD Map, MAG_AP, float, average mag in this band, mags FIELD dMap, MAG_AP_ERR, float, formal error on average mag, mags +FIELD sMap, MAG_AP_STDEV, float, standard deviation of ap mags, mags FIELD Mkron, MAG_KRON, float, ave kron mag in this band, mags FIELD dMkron, MAG_KRON_ERR, float, formal error on average kron mag, mags +FIELD sMkron, MAG_KRON_STDEV, float, standard deviation of kron mags, mags # these statistics are PSF-specific Index: /branches/eam_branches/ipp-20140610/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/libdvo/include/dvo.h (revision 36916) +++ /branches/eam_branches/ipp-20140610/Ohana/src/libdvo/include/dvo.h (revision 36917) @@ -599,4 +599,5 @@ /* pointers for data manipulation */ off_t *found; + off_t *foundWarp; off_t *image; off_t *mosaic; Index: /branches/eam_branches/ipp-20140610/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/libdvo/src/dvo_catalog.c (revision 36916) +++ /branches/eam_branches/ipp-20140610/Ohana/src/libdvo/src/dvo_catalog.c (revision 36917) @@ -166,6 +166,8 @@ secfilt->Map = NAN; secfilt->dMap = NAN; + secfilt->sMap = NAN; secfilt->Mkron = NAN; secfilt->dMkron = NAN; + secfilt->sMkron = NAN; secfilt->Mstdev = NAN; Index: /branches/eam_branches/ipp-20140610/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 36916) +++ /branches/eam_branches/ipp-20140610/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 36917) @@ -284,5 +284,6 @@ out[i].dM = in[i].dM; out[i].Map = in[i].Map; - out[i].dMap = in[i].dM; + out[i].dMap = in[i].dMap; + out[i].sMap = in[i].sMap; out[i].Mkron = in[i].Mkron; out[i].dMkron = in[i].dMkron; @@ -329,7 +330,7 @@ out[i].sFapWrp = in[i].sFapWrp; - out[i].NusedWrp = in[i].Nused; - out[i].NusedKronWrp = in[i].NusedKron; - out[i].NusedApWrp = in[i].NusedAp; + out[i].NusedWrp = in[i].NusedWrp; + out[i].NusedKronWrp = in[i].NusedKronWrp; + out[i].NusedApWrp = in[i].NusedApWrp; out[i].ubercalDist = in[i].ubercalDist; @@ -350,7 +351,9 @@ out[i].dM = in[i].dM; out[i].Map = in[i].Map; - out[i].dMap = in[i].dM; + out[i].dMap = in[i].dMap; + out[i].sMap = in[i].sMap; out[i].Mkron = in[i].Mkron; out[i].dMkron = in[i].dMkron; + out[i].sMkron = in[i].sMkron; out[i].Mstdev = in[i].Mstdev; @@ -395,7 +398,7 @@ out[i].sFapWrp = in[i].sFapWrp; - out[i].NusedWrp = in[i].Nused; - out[i].NusedKronWrp = in[i].NusedKron; - out[i].NusedApWrp = in[i].NusedAp; + out[i].NusedWrp = in[i].NusedWrp; + out[i].NusedKronWrp = in[i].NusedKronWrp; + out[i].NusedApWrp = in[i].NusedApWrp; out[i].ubercalDist = in[i].ubercalDist; Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/StarOps.c (revision 36916) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/StarOps.c (revision 36917) @@ -153,4 +153,7 @@ int i; + fprintf (stderr, "this version of setMrel is invalid for now\n"); + exit (1); + int Nsecfilt = GetPhotcodeNsecfilt (); Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelCatalog.c (revision 36916) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelCatalog.c (revision 36917) @@ -1,3 +1,4 @@ # include "relphot.h" +# define UBERCAL_WEIGHT 100.0 # define SKIP_THIS_MEAS(REASON) { \ @@ -13,4 +14,9 @@ continue; } +static float MagToFlux (float Mag) { + float Flux = pow(10.0, -0.4*(Mag)); + return (Flux); +} + int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) { @@ -29,4 +35,6 @@ return (TRUE); } + +int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found); int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { @@ -59,4 +67,10 @@ setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]); } + + // only measure force-warp mean values if issetMrelFinal (make it optional?) + if (isSetMrelFinal) { + setMrelAverageForcedWarp (m, Nc, pass, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], &catalog[Nc].foundWarp[Nsecfilt*j]); + setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]); + } } if (primaryCell) free (primaryCell); @@ -64,4 +78,5 @@ } +// set mean of chip measurements (selected by photcode range for now): int setMrelAverageExposure (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) { @@ -199,6 +214,11 @@ } - // gpc1 stack data - if ((measure[k].photcode >= 11000) && (measure[k].photcode <= 11400)) { + // SKIP gpc1 stack data + if ((measure[k].photcode >= 11000) && (measure[k].photcode <= 11500)) { + continue; + } + + // SKIP gpc1 forced-warp data + if ((measure[k].photcode >= 12000) && (measure[k].photcode <= 12500)) { continue; } @@ -258,6 +278,6 @@ if (!isnan(Map)) { Maplist[Nap] = Map - Mcal - Mmos - Mgrid; - daplist[Nap] = dMpsf; - waplist[Nap] = isUbercal ? 10.0 : 1.0; + daplist[Nap] = dMpsf;// XXX check on this... + waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0; Nap ++; } @@ -267,5 +287,5 @@ Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid; dkronlist[Nkron] = measure[k].dMkron; - wkronlist[Nkron] = isUbercal ? 10.0 : 1.0; + wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; Nkron ++; } @@ -280,5 +300,5 @@ Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid; dpsflist[Npsf] = dMpsf; - wpsflist[Npsf] = isUbercal ? 10.0 : 1.0; + wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0; // NOTE: @@ -293,5 +313,5 @@ if (refPhotcode) { if (code->code == refPhotcode->code) { - wpsflist[Npsf] = 100.0; + wpsflist[Npsf] = UBERCAL_WEIGHT; } } @@ -304,5 +324,5 @@ if (isSetMrelFinal && (pass == 0)) { if ((thisCode < 6) || (thisCode == 9)) { - secfilt[Nsec].Ncode = NexpPS1; + secfilt[Nsec].Ncode = NexpPS1; } else { secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists @@ -346,18 +366,23 @@ found[Nsec] = TRUE; - secfilt[Nsec].Mstdev = 1000.0*psfstats->sigma; // Mstdev is in millimags (not enough space for more precision) + secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision) // secfilt[Nsec].Ncode = Ncode; secfilt[Nsec].Nused = psfstats->Nmeas; - secfilt[Nsec].Mmax = 1000 * psfstats->max; - secfilt[Nsec].Mmin = 1000 * psfstats->min; + secfilt[Nsec].Mmax = psfstats->max; + secfilt[Nsec].Mmin = psfstats->min; // NOTE : use the modified weight for apmags as well as psf mags liststats (Maplist, daplist, waplist, Nap, apstats); secfilt[Nsec].Map = Nap > 0 ? apstats->mean : NAN; + secfilt[Nsec].dMap = Nap > 0 ? apstats->error : NAN; + secfilt[Nsec].sMap = Nap > 0 ? apstats->sigma : NAN; + secfilt[Nsec].NusedAp = Nap; liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats); secfilt[Nsec].Mkron = Nkron > 0 ? kronstats->mean : NAN; secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN; + secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN; + secfilt[Nsec].NusedKron = Nkron; // NOTE: for 2MASS measurements, Next should be 1, as should N @@ -585,4 +610,10 @@ secfilt[Nsec].FkronStk = zpFactor * measure[k].FluxKron; secfilt[Nsec].dFkronStk = zpFactor * measure[k].dFluxKron; + secfilt[Nsec].FapStk = zpFactor * measure[k].FluxAp; + 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].stackDetectID = ID; @@ -607,4 +638,153 @@ } +# undef SKIP_THIS_MEAS +# define SKIP_THIS_MEAS(REASON) { \ + measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + results->REASON ++; \ + continue; } + +// set mean of forced-warp measurements (selected by photcode range for now): +// somewhat simplified relative to chip-photometry: +// * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr +int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found) { + + off_t k; + float Fsys = 0, Mcal= 0; + + // we are measuring means for 3 types of FLUXes: psf, ap, kron. I am using the psf mag + // error for the ap mags, but krons have their own errors. it is an open question if I + // should be doing weighted or unweighted fits (this is a user option) + double *Fpsflist = results->Mpsflist; + double *dpsflist = results->dpsflist; + double *wpsflist = results->wpsflist; + + double *Faplist = results->Maplist; + double *daplist = results->daplist; + double *waplist = results->waplist; + + double *Fkronlist = results->Mkronlist; + double *dkronlist = results->dkronlist; + double *wkronlist = results->wkronlist; + + StatType *psfstats = &results->psfstats; + StatType *apstats = &results->apstats; + StatType *kronstats = &results->kronstats; + + // option for a test print + if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) { + fprintf (stderr, "test obj\n"); + print_measure_set_alt (average, secfilt, measure); + } + + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + /* calculate the average mag in this SEC photcode for a single star */ + + /* star/photcodes already calibrated */ + if (found[Nsec]) continue; + + off_t meas = measureOffset; + + int Nap = 0; + int Npsf = 0; + int Nkron = 0; + for (k = 0; k < average[0].Nmeasure; k++, meas++) { + + // skip measurements which are not ForcedWarp values + if ((measure[k].photcode < 12000) || (measure[k].photcode > 12500)) continue; + + // skip measurements that do not match the current photcode + PhotCode *code = GetPhotcodebyCode (measure[k].photcode); + if (!code) continue; + if (code->equiv != thisCode) { continue; } + + if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); + + if (getImageEntry (meas, cat) < 0) { + // measurements without an image are either external reference photometry or + // data for which the associated image has not been loaded (probably because of + // overlaps). Msys + measure.Mcal is our best guess of the true magnitude + Mcal = measure[k].Mcal; // check that this is zero for loaded REF value + } else { + // use getMcal not getMcal_alt? + Mcal = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd); + // Mcal = getMcal (meas, cat); + if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); + } + float Fcal = MagToFlux(Mcal); + + // in the calculations below, + // ...list gives the error per measurement, wlist gives the weight + // we can modify the error and weight in a few ways: + // 1) MIN_ERROR guarantees a floor + // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode + + // skip some absurd values NAN + // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms + // are small (and not measured) + float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF); + if (isnan(Fsys)) SKIP_THIS_MEAS(Nsys); + // if (Msys < 0.0) SKIP_THIS_MEAS(Nsys); + // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys); + + float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF); + dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error + Fpsflist[Npsf] = Fpsf * Fcal; + dpsflist[Npsf] = dFpsf * Fcal; + wpsflist[Npsf] = 1.0; + Npsf ++; + + float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER); + float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER); + if (!isnan(Fap)) { + Faplist[Nap] = Fap * Fcal; + daplist[Nap] = dFap * Fcal; + waplist[Nap] = 1.0; // drop weight lists? + Nap ++; + } + + float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON); + float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON); + if (!isnan(Fkron)) { + Fkronlist[Nkron] = Fkron * Fcal; + dkronlist[Nkron] = dFkron * Fcal; + wkronlist[Nkron] = 1.0; + Nkron ++; + } + } + if (Npsf < 1) continue; + + found[Nsec] = TRUE; + + liststats (Fpsflist, dpsflist, wpsflist, Npsf, psfstats); + + secfilt[Nsec].FpsfWrp = psfstats->mean; + secfilt[Nsec].dFpsfWrp = psfstats->error; + 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 + + // NOTE : use the modified weight for apmags as well as psf mags + liststats (Faplist, daplist, waplist, Nap, apstats); + secfilt[Nsec].FapWrp = Nap > 0 ? apstats->mean : NAN; + secfilt[Nsec].dFapWrp = Nap > 0 ? apstats->error : NAN; + 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 + + liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats); + secfilt[Nsec].FkronWrp = Nkron > 0 ? kronstats->mean : NAN; + secfilt[Nsec].dFkronWrp = Nkron > 0 ? kronstats->error : NAN; + 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 + } + return (TRUE); +} + int setGlobalObjStats (Average *average, Measure *measure) { Index: /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelFinal.c (revision 36916) +++ /branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelFinal.c (revision 36917) @@ -86,6 +86,8 @@ /* set catalog[0].found[i] = FALSE */ ALLOCATE (catalog[0].found, off_t, MAX (1, Nsecfilt*catalog[0].Naverage)); + ALLOCATE (catalog[0].foundWarp, off_t, MAX (1, Nsecfilt*catalog[0].Naverage)); for (i = 0; i < Nsecfilt*catalog[0].Naverage; i++) { catalog[0].found[i] = FALSE; + catalog[0].foundWarp[i] = FALSE; }