IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 14, 2014, 2:28:47 PM (12 years ago)
Author:
bills
Message:

Iniital implementation of the full force stages to the pipeline.
Maginitude and galactic coordinate limits on extended source fits

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psphot

  • trunk/psphot/src

  • trunk/psphot/src/psphotChooseAnalysisOptions.c

    r36375 r36441  
    6060}
    6161
    62 bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos);
     62static bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos);
    6363
    6464// this function use an internal flag to mark sources which have already been measured
    … …  
    121121    psSphereRot *toGal = NULL;
    122122    float GAL_LIMIT = 0.0;
     123    float GAL_LIMIT_BULGE = 0.0;
     124    float GAL_LIMIT_SIGMA2 = 0.0;
    123125    bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE");
    124126    assert (status);
    125127    if (useGAL_LIMIT) {
     128        toGal = psSphereRotICRSToGalactic();
    126129        GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT");
    127         toGal = psSphereRotICRSToGalactic();
     130        assert (status);
     131        GAL_LIMIT = DEG_TO_RAD(GAL_LIMIT);
     132        GAL_LIMIT_BULGE = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE");
     133        assert (status);
     134        GAL_LIMIT_BULGE = DEG_TO_RAD(GAL_LIMIT_BULGE);
     135        float GAL_LIMIT_SIGMA = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE.SIGMA");
     136        assert (status);
     137        assert(GAL_LIMIT_SIGMA > 0);
     138        GAL_LIMIT_SIGMA = DEG_TO_RAD(GAL_LIMIT_SIGMA);
     139        GAL_LIMIT_SIGMA2 = GAL_LIMIT_SIGMA * GAL_LIMIT_SIGMA;
    128140    }
    129141
    … …  
    135147    }
    136148
     149    float petroFluxLim = NAN;
    137150    float extFitFluxLim = NAN;
    138151
    … …  
    140153        float extFitMagLimDefault = NAN;
    141154        float extFitMagLim = NAN;
     155        float petroMagLimDefault = NAN;
     156        float petroMagLim = NAN;
    142157
    143158        // match to the given filter
    … …  
    156171            // find a matching filter or default to 'any'
    157172            if (!strcasecmp (thisFilter, "any")) {
    158                 psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
    159                 psAssert(extFitMagLimStr, "missing MAG.LIMIT");
     173                psString petroMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO");
     174                psAssert(petroMagLimStr, "missing MAG.LIMIT.PETRO");
     175                petroMagLimDefault = atof (petroMagLimStr);
     176
     177                psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT");
     178                psAssert(extFitMagLimStr, "missing MAG.LIMIT.EXTFIT");
    160179                extFitMagLimDefault = atof (extFitMagLimStr);
    161180            }
    162181
    163             // find a matching filter or default to 'any'
    164182            if (!strcasecmp (thisFilter, filterID)) {
    165                 psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
    166                 psAssert(extFitMagLimStr, "missing MAG.LIMIT");
     183                psString petroMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO");
     184                psAssert(petroMagLimStr, "missing MAG.LIMIT.PETRO");
     185                petroMagLim = atof (petroMagLimStr);
     186
     187                psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT");
     188                psAssert(extFitMagLimStr, "missing MAG.LIMIT.EXTFIT");
    167189                extFitMagLim = atof (extFitMagLimStr);
    168190                break;
    169191            }
    170192        }
     193        psFree(iter);
     194        if (!isfinite (petroMagLim)) petroMagLim = petroMagLimDefault;
    171195        if (!isfinite (extFitMagLim)) extFitMagLim = extFitMagLimDefault;
    172196
    … …  
    185209        psAssert (status, "missing FPA.ZP?");
    186210
     211        petroFluxLim = exptime * pow (10.0, 0.4*(zeropt - petroMagLim));
    187212        extFitFluxLim = exptime * pow (10.0, 0.4*(zeropt - extFitMagLim));
    188213    }
    … …  
    225250            psSphere ptGal, ptSky;
    226251            GetGalacticCoords (&ptGal, &ptSky, toGal, chip, source->peak->xf, source->peak->yf);
    227             if (fabs(ptGal.d) < GAL_LIMIT) continue;
     252            float b = ptGal.r;
     253            float limit = GAL_LIMIT + GAL_LIMIT_BULGE * exp(-0.5*(b*b/GAL_LIMIT_SIGMA2));
     254            if (fabs(ptGal.d) < limit) continue;
    228255            // include an exception for low density skycells below the limit?
    229256        }
    230257
    231258        // for petro and extFit, we will either use the mag limits or the S/N
     259        if (doPetrosian) {
     260            if (isfinite(petroFluxLim)) {
     261                if (source->moments->KronFlux > petroFluxLim) {
     262                    source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
     263                }
     264            } else if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) {
     265                source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
     266            }
     267        }
    232268        if (isfinite(extFitFluxLim)) {
    233269            if (source->moments->KronFlux > extFitFluxLim) {
    234270                source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
    235                 if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
    236271            }
    237272        } else {
    238273            if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) {
    239274                source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
    240                 if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
    241275            }
    242276        }
    … …  
    244278
    245279    psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.options"));
     280
     281    psFree(toGal);
    246282
    247283    return true;
    … …  
    287323    psSphereRot *toGal = NULL;
    288324    float GAL_LIMIT = 0.0;
     325    float GAL_LIMIT_BULGE = 0.0;
     326    float GAL_LIMIT_SIGMA2 = 0.0;
    289327    bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE");
    290328    assert (status);
    291329    if (useGAL_LIMIT) {
     330        toGal = psSphereRotICRSToGalactic();
    292331        GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT");
    293         toGal = psSphereRotICRSToGalactic();
     332        assert (status);
     333        GAL_LIMIT = DEG_TO_RAD(GAL_LIMIT);
     334        GAL_LIMIT_BULGE = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE");
     335        assert (status);
     336        GAL_LIMIT_BULGE = DEG_TO_RAD(GAL_LIMIT_BULGE);
     337        float GAL_LIMIT_SIGMA = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE.SIGMA");
     338        assert (status);
     339        assert(GAL_LIMIT_SIGMA > 0);
     340        GAL_LIMIT_SIGMA = DEG_TO_RAD(GAL_LIMIT_SIGMA);
     341        GAL_LIMIT_SIGMA2 = GAL_LIMIT_SIGMA * GAL_LIMIT_SIGMA;
    294342    }
    295343
    … …  
    365413    // find extFitFluxLim->data.F32[i] for i == image number
    366414    psVector *extFitFluxLim = NULL;
     415    psVector *petroFluxLim = NULL;
    367416    if (magLimits) {
    368417        extFitFluxLim = psVectorAlloc (inputFilters->n, PS_TYPE_F32);
     418        petroFluxLim = psVectorAlloc (inputFilters->n, PS_TYPE_F32);
    369419        psVectorInit (extFitFluxLim, NAN);
     420        psVectorInit (petroFluxLim, NAN);
    370421
    371422        float extFitMagLimDefault = NAN;
     423        float petroMagLimDefault = NAN;
    372424
    373425        // match mag limits (flux limits) to the filters
    … …  
    383435            // save the default value to assign to unset filters
    384436            if (!strcasecmp (thisFilter, "any")) {
    385                 psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
    386                 psAssert(magString, "missing MAG.LIMIT");
     437                psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT");
     438                psAssert(magString, "missing MAG.LIMIT.EXTFIT");
    387439                extFitMagLimDefault = atof (magString);
     440
     441                magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO");
     442                psAssert(magString, "missing MAG.LIMIT.PETRO");
     443                petroMagLimDefault = atof (magString);
    388444                continue;
    389445            }
    … …  
    393449                if (i == chisqNum) continue;
    394450                if (!strcasecmp (thisFilter, inputFilters->data[i])) {
    395                     psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
    396                     psAssert(magString, "missing MAG.LIMIT");
     451                    psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO");
     452                    psAssert(magString, "missing MAG.LIMIT.PETRO");
    397453                    float magvalue = atof (magString);
     454                    petroFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - magvalue));
     455
     456                    magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT");
     457                    psAssert(magString, "missing MAG.LIMIT.EXTFIT");
     458                    magvalue = atof (magString);
    398459                    extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - magvalue));
    399460                    break;
    … …  
    401462            }
    402463        }
     464        psFree(iter);
    403465
    404466        for (int i = 0; i < num; i++) {
    405467            if (i == chisqNum) continue;
    406             if (isfinite(extFitFluxLim->data.F32[i])) continue;
    407             extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - extFitMagLimDefault));
     468            if (!isfinite(petroFluxLim->data.F32[i])) {
     469                petroFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - petroMagLimDefault));
     470            }
     471            if (!isfinite(extFitFluxLim->data.F32[i])) {
     472                extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - extFitMagLimDefault));
     473            }
    408474        }
    409475    }
    … …  
    420486        bool doObjectRadial = false;
    421487        bool doObjectExtFit = false;
    422         for (int j = 0; !doObjectExtFit && !doObjectRadial && (j < object->sources->n); j++) {
     488        bool doObjectPetrosian = false;
     489        for (int j = 0; !doObjectExtFit && !doObjectRadial && !doObjectPetrosian && (j < object->sources->n); j++) {
    423490
    424491            pmSource *source = object->sources->data[j];
    … …  
    451518            if (extFitAll) {
    452519                doObjectExtFit = true;
     520                doObjectPetrosian = doPetrosian;
    453521                continue;
    454522            }
    … …  
    465533                psSphere ptGal, ptSky;
    466534                GetGalacticCoords (&ptGal, &ptSky, toGal, chips->data[imageID], source->peak->xf, source->peak->yf);
    467                 if (fabs(ptGal.d) < GAL_LIMIT) continue;
    468                 // include an exception for low density skycells below the limit?
    469             }
    470 
    471             float fluxLim = extFitFluxLim ? extFitFluxLim->data.F32[imageID] : NAN;
     535                float b = ptGal.r;
     536                float limit = GAL_LIMIT + GAL_LIMIT_BULGE * exp(-0.5*(b*b/GAL_LIMIT_SIGMA2));
     537                if (fabs(ptGal.d) < limit) continue;
     538            }
     539
     540            float fluxLim = NAN;
     541            // for petro and extFit, we will either use the mag limits or the S/N
     542            if (doPetrosian) {
     543                fluxLim = petroFluxLim ? petroFluxLim->data.F32[imageID] : NAN;
     544                if (isfinite(fluxLim)) {
     545                    if (source->moments->KronFlux > extFitFluxLim->data.F32[imageID]) {
     546                        doObjectPetrosian = true;
     547                    }
     548                } else {
     549                    if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) {
     550                        doObjectPetrosian = true;
     551                    }
     552                }
     553            }
     554            fluxLim = extFitFluxLim ? extFitFluxLim->data.F32[imageID] : NAN;
    472555            // for petro and extFit, we will either use the mag limits or the S/N
    473556            if (isfinite(fluxLim)) {
    … …  
    503586            if (source->modelPSF == NULL) continue;
    504587               
    505             // Do the fits if the recipe requests we do extended source fits to everything
    506588            if (doObjectExtFit) {
    507589                source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
    508                 if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
    509             }
     590            }
     591            if (doObjectPetrosian) {
     592                source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
     593            }
    510594
    511595            // Do the fits if the recipe requests we do extended source fits to everything
    … …  
    522606    psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld objects: %f sec\n", objects->n, psTimerMark ("psphot.options"));
    523607
     608    psFree(toGal);
     609
    524610    return true;
    525611}
    526612
    527 bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos) {
     613static bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos) {
    528614
    529615    pmFPA *fpa = chip->parent;
    … …  
    543629    psSphereRotApply (ptGal, toGal, ptSky);
    544630
     631    // psSphereRotApply insures that 0 < r < 2PI. We want -PI < b <= PI
     632    if (ptGal->r > M_PI) {
     633        ptGal->r -= 2.0 * M_PI;
     634    }
     635
    545636    return true;
    546637
Note: See TracChangeset for help on using the changeset viewer.