IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 31362


Ignore:
Timestamp:
Apr 24, 2011, 10:26:57 PM (15 years ago)
Author:
eugene
Message:

change growth curve calculation to use sources not psf model; calculate and use minimum kron radius (that appropriate to stars); errMag renamed psfMagErr

Location:
branches/eam_branches/ipp-20110404/psphot/src
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110404/psphot/src/psphot.h

    r31313 r31362  
    313313
    314314bool psphotMakeFluxScale (psImage *image, psMetadata *recipe, pmPSF *psf);
    315 bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf);
     315bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources);
    316316bool psphotDumpPSFStars (pmReadout *readout, pmPSFtry *try, float radius, psImageMaskType maskVal, psImageMaskType markVal);
    317317

  • branches/eam_branches/ipp-20110404/psphot/src/psphotApResid.c

    r31154 r31362  
    206206
    207207        // XXX make this user-configurable?
    208         if (source->errMag > 0.01) continue;
     208        if (source->psfMagErr > 0.01) continue;
    209209
    210210        // aperture residual for this source
    … …  
    222222                 source->peak->xf, source->peak->yf,
    223223                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS],
    224                  source->psfMag, source->apMag, source->errMag,
     224                 source->psfMag, source->apMag, source->psfMagErr,
    225225                 source->modelPSF->params->data.F32[PM_PAR_I0],
    226226                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY],
    … …  
    228228# endif
    229229        if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
    230         if (!isfinite(source->errMag)) psAbort ("nan in errMag");
     230        if (!isfinite(source->psfMagErr)) psAbort ("nan in psfMagErr");
    231231        if (!isfinite(source->apMag)) psAbort ("nan in apMag");
    232232        if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
    … …  
    234234
    235235        psVectorAppend (mag, source->psfMag);
    236         psVectorAppend (dMag,source->errMag);
     236        psVectorAppend (dMag,source->psfMagErr);
    237237        psVectorAppend (apResid, dap);
    238238        psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]);

  • branches/eam_branches/ipp-20110404/psphot/src/psphotChoosePSF.c

    r31154 r31362  
    350350
    351351    // build curve-of-growth vector for this psf
    352     if (!psphotMakeGrowthCurve (readout, recipe, try->psf)) {
     352    if (!psphotMakeGrowthCurve (readout, recipe, try->psf, try->sources)) {
    353353        psError(PSPHOT_ERR_PSF, false, "Unable to construct flux scale for PSF");
    354354        psFree (models);

  • branches/eam_branches/ipp-20110404/psphot/src/psphotEfficiency.c

    r30560 r31362  
    468468                    source->modelPSF->params->data.F32[PM_PAR_XPOS],
    469469                    source->modelPSF->params->data.F32[PM_PAR_YPOS],
    470                     magRef, source->psfMag, source->errMag);
     470                    magRef, source->psfMag, source->psfMagErr);
    471471#endif
    472472            magDiff->data.F32[j] = source->psfMag - magRef;
    473             magErr->data.F32[j] = source->errMag;
     473            magErr->data.F32[j] = source->psfMagErr;
    474474        }
    475475        magDiff->n = sources->n;

  • branches/eam_branches/ipp-20110404/psphot/src/psphotExtendedSourceFits.c

    r31154 r31362  
    333333        // this uses the footprint to judge both radius and aperture?
    334334        // XXX save the psf-based moments for output
    335         if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) {
     335        if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) {
    336336            fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
    337337            // subtract the best fit from the object, leave local sky

  • branches/eam_branches/ipp-20110404/psphot/src/psphotKronMasked.c

    r31337 r31362  
    131131        }
    132132    }           
    133     psphotSaveImage (NULL, kronMask, "kronmask.fits");
     133    // psphotSaveImage (NULL, kronMask, "kronmask.fits");
    134134    psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
    135135

  • branches/eam_branches/ipp-20110404/psphot/src/psphotLoadSRCTEXT.c

    r31337 r31362  
    6767
    6868            source->psfMag    = 0.0;
    69             source->errMag    = 0.0;
     69            source->psfMagErr    = 0.0;
    7070            source->apMag     = 0.0;
    7171

  • branches/eam_branches/ipp-20110404/psphot/src/psphotMakeGrowthCurve.c

    r31154 r31362  
    11# include "psphotInternal.h"
    22
    3 bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf) {
     3bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources) {
    44
    55    bool status;
    … …  
    2222
    2323    // measure the aperture loss as a function of radius for PSF
    24     bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
    2524    float REF_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_REF_RADIUS");
    2625    float PSF_FIT_PAD   = psMetadataLookupF32 (&status, recipe, "PSF_FIT_PADDING");
     26   
     27    float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
     28    if (!status) {
     29        gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
     30    }
     31    float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
     32    float PSF_APERTURE = (int)(apScale*gaussSigma);
    2733
    2834    psf->growth = pmGrowthCurveAlloc (PSF_FIT_PAD, 100.0, REF_RADIUS);
    2935
     36# if (0)
     37    bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
    3038    if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH, maskVal, markVal)) {
    3139        psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
    … …  
    3341        return false;
    3442    }
     43# else
     44    bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
     45    if (!pmGrowthCurveGenerateFromSources (readout, psf, sources, INTERPOLATE_AP, maskVal, markVal)) {
     46        psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
     47        psFree(psf->growth); psf->growth = NULL;
     48        return false;
     49    }
     50# endif
     51
     52    float offset = pmGrowthCurveCorrect (psf->growth, PSF_APERTURE);
     53    fprintf (stderr, "correction from %f to %f: %f mags\n", PSF_APERTURE, REF_RADIUS, offset);
    3554
    3655    psLogMsg ("psphot", PS_LOG_MINUTIA, "built growth curve: %f sec\n", psTimerMark ("psphot.growth"));

  • branches/eam_branches/ipp-20110404/psphot/src/psphotOutput.c

    r31154 r31362  
    7575                 source->peak->xf, source->peak->yf,
    7676                 model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS],
    77                  source->psfMag, source->apMag, source->errMag,
     77                 source->psfMag, source->apMag, source->psfMagErr,
    7878                 model->params->data.F32[PM_PAR_I0],
    7979                 model->params->data.F32[PM_PAR_SXX], model->params->data.F32[PM_PAR_SXY], model->params->data.F32[PM_PAR_SYY],

  • branches/eam_branches/ipp-20110404/psphot/src/psphotSetThreads.c

    r31154 r31362  
    2525    psFree(task);
    2626
    27     task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 10);
     27    task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 11);
    2828    task->function = &psphotSourceStats_Threaded;
    2929    psThreadTaskAdd(task);

  • branches/eam_branches/ipp-20110404/psphot/src/psphotSourceFits.c

    r31313 r31362  
    235235    // XXX save the psf-based moments for output
    236236    psfMoments = *source->moments;
    237     if (!pmSourceMoments (source, radius, 0.0, 0.5, maskVal)) {
     237    if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
    238238      *source->moments = psfMoments;
    239239      return false;

  • branches/eam_branches/ipp-20110404/psphot/src/psphotSourceSize.c

    r31154 r31362  
    213213
    214214        psVectorAppend (ApOff, dMag);
    215         psVectorAppend (ApErr, source->errMag);
     215        psVectorAppend (ApErr, source->psfMagErr);
    216216    }
    217217    if (num == 0) {
    … …  
    520520        // set nSigmaMAG to include both systematic and poisson error terms.  we include a hard
    521521        // floor on the Ap Sys Err (to be a bit generous).  XXX put the floor in the recipe...
    522         float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.02));
     522        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->psfMagErr, hypot(options->ApSysErr, 0.02));
    523523        source->extNsigma = nSigmaMAG;
    524524

  • branches/eam_branches/ipp-20110404/psphot/src/psphotSourceStats.c

    r31328 r31362  
    165165        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    166166    }
     167    float MIN_KRON_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_MIN_KRON");
     168    if (!status) {
     169        MIN_KRON_RADIUS = 0.75*SIGMA;
     170    }
    167171
    168172    // threaded measurement of the source magnitudes
    … …  
    192196            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,  PS_TYPE_F32);
    193197            PS_ARRAY_ADD_SCALAR(job->args, SIGMA,   PS_TYPE_F32);
     198            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32);
    194199
    195200            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
    … …  
    221226            } else {
    222227                psScalar *scalar = NULL;
    223                 scalar = job->args->data[7];
     228                scalar = job->args->data[8];
    224229                Nmoments += scalar->data.S32;
    225                 scalar = job->args->data[8];
     230                scalar = job->args->data[9];
    226231                Nfail += scalar->data.S32;
    227                 scalar = job->args->data[9];
     232                scalar = job->args->data[10];
    228233                Nfaint += scalar->data.S32;
    229234            }
    … …  
    382387    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[3],F32);
    383388    float SIGMA                     = PS_SCALAR_VALUE(job->args->data[4],F32);
    384 
    385     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    386     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
     389    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[5],F32);
     390
     391    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
     392    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
    387393
    388394    // if no valid pixels, massive swing or very large Mrf, likely saturated source,
    … …  
    440446        if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) {
    441447            // measure basic source moments (no S/N clipping on input pixels)
    442             status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal);
     448            status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, MIN_KRON_RADIUS, maskVal);
    443449        } else {
    444450            // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but
    … …  
    456462
    457463            psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    458             status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, maskVal);
     464            status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, MIN_KRON_RADIUS, maskVal);
    459465            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
    460466        }
    … …  
    507513    float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0};
    508514    float Sout[NSIGMA];
     515    float Rmin[NSIGMA];
    509516    int   Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass
    510517
    … …  
    528535
    529536            // measure basic source moments (no S/N clipping on input pixels)
    530             status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, maskVal);
     537            status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, 0.0, maskVal);
    531538        }
    532539
    … …  
    539546        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);
    540547        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
     548
     549        Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);
    541550
    542551#if 0
    … …  
    564573    // we are looking for sigma for which Sout = 0.65 (or some other value)
    565574    int Nstars = 0;
     575    float minKronRadius = NAN;
    566576    float Sigma = NAN;
    567577    float minS = Sout[0];
    … …  
    580590        Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
    581591        Nstars = 0.5*(Nout[i] + Nout[i+1]);
     592        minKronRadius = Rmin[i] + (0.65 - Sout[i])*(Rmin[i+1] - Rmin[i])/(Sout[i+1] - Sout[i]);
    582593    }
    583594    psAssert (isfinite(Sigma), "did we miss a case?");
    … …  
    590601    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
    591602    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_NSTARS", PS_META_REPLACE, "number of stars in clump", Nstars);
    592 
    593     psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
     603    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_MIN_KRON", PS_META_REPLACE, "minimum Kron Radius", minKronRadius);
     604
     605    psLogMsg ("psphot", 3, "using window function with sigma = %f (Min Kron Radius = %f, Nstars = %d)\n", Sigma, minKronRadius, Nstars);
    594606    return true;
    595607}

  • branches/eam_branches/ipp-20110404/psphot/src/psphotVisual.c

    r31154 r31362  
    18671867            continue;
    18681868        }
    1869         if (source->errMag > 0.1) {
     1869        if (source->psfMagErr > 0.1) {
    18701870            xLOW->data.F32[nLOW] = source->moments->Mxx;
    18711871            yLOW->data.F32[nLOW] = source->moments->Myy;
    … …  
    25912591        x->data.F32[n] = source->psfMag;
    25922592        y->data.F32[n] = dMag;
    2593         dy->data.F32[n] = source->errMag;
     2593        dy->data.F32[n] = source->psfMagErr;
    25942594        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
    25952595        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
Note: See TracChangeset for help on using the changeset viewer.