IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34086


Ignore:
Timestamp:
Jun 26, 2012, 11:33:10 AM (14 years ago)
Author:
eugene
Message:

re-enable MODEL_VAR option for linear photometry fit

Location:
trunk/psphot
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/configure.ac

    r23805 r34086  
    99AM_MAINTAINER_MODE
    1010
    11 IPP_STDCFLAGS
     11IPP_STDLDFLAGS
    1212
    1313AC_LANG(C)
    … …  
    199199dnl Set CFLAGS for build
    200200IPP_STDOPTS
    201 
    202 CFLAGS="${CFLAGS=} -Wall -Werror"
     201IPP_STDCFLAGS
     202
    203203echo "PSPHOT_CFLAGS: $PSPHOT_CFLAGS"
    204204echo "PSPHOT_LIBS: $PSPHOT_LIBS"
    … …  
    206206IPP_VERSION
    207207
     208dnl note: the PSPHOT_ entries below do not include PSLIB_ and PSMODULE_ entires
     209dnl those are passed down to lower Makefiles by the PKG_CHECK_MODULES call above
    208210AC_SUBST([PSPHOT_CFLAGS])
    209211AC_SUBST([PSPHOT_LIBS])

  • trunk/psphot/src

  • trunk/psphot/src/Makefile.am

    r33690 r34086  
    183183        psphotSourcePlots.c            \
    184184        psphotRadialPlot.c             \
    185         psphotKronMasked.c             \
    186185        psphotKronIterate.c            \
    187186        psphotRadialProfileWings.c     \
    … …  
    209208        psphotSetNFrames.c
    210209
    211 # re-instate these
     210# not currently used
    212211#       psphotIsophotal.c              \
    213212#       psphotAnnuli.c                 \
    214213#       psphotKron.c                   \
     214#       psphotKronMasked.c             \
    215215#
    216216

  • trunk/psphot/src/psphot.h

    r33994 r34086  
    6161bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    6262
     63bool            psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule);
     64bool            psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     65
    6366bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
    6467bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    … …  
    99102
    100103bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
    101 
    102 # if (HAVE_MODEL_VAR)
    103104bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
    104 # else
    105 bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
    106 # endif
    107105
    108106bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);
    … …  
    489487bool psphotSetNFramesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    490488
     489bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
     490pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
     491bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
     492
     493bool psphotModelBackgroundReadout(psImage *model,  // Model image
     494                                  psImage *modelStdev, // Model stdev image
     495                                  psMetadata *analysis, // Analysis metadata for outputs
     496                                  pmReadout *readout, // Readout for which to generate a background model
     497                                  psImageBinning *binning, // Binning parameters
     498                                  const pmConfig *config,// Configuration
     499                                  bool useVarianceImage
     500    );
     501
    491502#endif

  • trunk/psphot/src/psphotCleanup.c

    r29548 r34086  
    2929psExit psphotGetExitStatus (void) {
    3030
    31     psErrorCode err = psErrorCodeLast ();
     31    // gcc -Wswitch complains here if err is declared as type psErrorCode
     32    // the collection of ps*ErrorCode values are enums defined separately for
     33    // each module (psphot, pswarp, etc).  the lowest type, psErrorCode is only the base set and does
     34    // not include the possible psphot values
     35
     36    // for now, to get around this, we just use an int for the switch
     37
     38    // psErrorCode err = psErrorCodeLast ();
     39    int err = psErrorCodeLast ();
    3240    switch (err) {
    3341      case PS_ERR_NONE:

  • trunk/psphot/src/psphotEfficiency.c

    r33993 r34086  
    420420
    421421    // psphotFitSourcesLinearReadout subtracts the model fits
    422 # if (HAVE_MODEL_VAR)
    423422    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
    424 # else
    425     if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
    426 # endif
    427423        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
    428424        psFree(fakeSources);

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r33994 r34086  
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
    1313
    14 # if (HAVE_MODEL_VAR)
    15 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
    16 pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
    17 bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
    18 # endif
    19 
    2014// for now, let's store the detections on the readout->analysis for each readout
    2115bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final)
    … …  
    3024    assert (recipe);
    3125
    32 # if (HAVE_MODEL_VAR)
    3326    pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
    3427    if (!fitVarMode) {
    … …  
    4033    // do a single pass.
    4134    pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
    42 # endif
    4335
    4436    int num = psphotFileruleCount(config, filerule);
    … …  
    6860        psAssert (psf, "missing psf?");
    6961
    70 # if (HAVE_MODEL_VAR)
    71         if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1))
    72 # else
    73         if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final))
    74 # endif
    75         {
     62        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) {
    7663            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
    7764            return false;
    7865        }
    7966
    80 # if (HAVE_MODEL_VAR)
    8167        // the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance
    8268        // after we have determined the initial set of fits, then we can generate the variance image and
    … …  
    10490            }
    10591        }
    106 # endif
    10792
    10893        psphotVisualShowResidualImage (readout, (num > 0));
    … …  
    11398}
    11499
    115 # if (HAVE_MODEL_VAR)
    116100// look up the fit variance mode from the recipe; older recipes do not have the value
    117101// 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and
    … …  
    136120        return PM_SOURCE_PHOTFIT_IMAGE_VAR;
    137121    }
     122    if (!strcasecmp(fitVarModeString, "SKY")   || !strcasecmp(fitVarModeString, "MODEL_SKY")) {
     123        return PM_SOURCE_PHOTFIT_MODEL_SKY;
     124    }
    138125    if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
    139126        return PM_SOURCE_PHOTFIT_MODEL_VAR;
    … …  
    142129    return PM_SOURCE_PHOTFIT_NONE;
    143130}
    144 # endif
    145 
    146 # if (HAVE_MODEL_VAR)
    147 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode)
    148 # else
    149 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final)
    150 # endif
    151 {
     131
     132bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) {
    152133    bool status;
    153134    float x;
    … …  
    185166    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    186167
    187 # if (!HAVE_MODEL_VAR)
    188     bool CONSTANT_PHOTOMETRIC_WEIGHTS =
    189         psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
    190     if (!status) {
    191         psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
    192     }
    193 # endif
    194168    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
    195169    if (!status) {
    … …  
    299273            source->mode |= PM_SOURCE_MODE_PSFMODEL;
    300274        }           
    301        
     275
    302276        psArrayAdd (fitSources, 100, source);
    303277    }
    … …  
    319293    psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
    320294
    321 # if (HAVE_MODEL_VAR)
    322295    // if fitVarMode is MODEL_VAR, then we need to generate the model image variance
    323296    // XXX we have two possibilities here:
    … …  
    327300
    328301    // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?)
    329 # endif
    330302
    331303    // fill out the sparse matrix elements and border elements (B)
    … …  
    336308
    337309        // diagonal elements of the sparse matrix (auto-cross-product)
    338 # if (HAVE_MODEL_VAR)
    339         f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
    340 # else
    341         f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    342 # endif
    343         psSparseMatrixElement (sparse, i, i, f);
    344 
    345 # if (HAVE_MODEL_VAR)
    346         // if we have used CONSTANT errors, then we need to calculate the value of the parameter error
     310        float MM = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
     311        psSparseMatrixElement (sparse, i, i, MM);
     312
     313        // if we have used CONSTANT errors, then we need to re-calculate the value of the parameter error
    347314        if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
    348315            float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
    349316            errors->data.F32[i] = 1.0 / sqrt(var);
    350317        } else {
    351             errors->data.F32[i] = 1.0 / sqrt(f);
    352         }
    353 # else
    354         // the formal error depends on the weighting scheme
    355         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    356             float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
    357             errors->data.F32[i] = 1.0 / sqrt(var);
    358         } else {
    359             errors->data.F32[i] = 1.0 / sqrt(f);
    360         }
    361 # endif
     318            errors->data.F32[i] = 1.0 / sqrt(MM);
     319        }
    362320
    363321        // find the image x model value
    364 # if (HAVE_MODEL_VAR)
    365         f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
    366 # else
    367         f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    368 # endif
    369         psSparseVectorElement (sparse, i, f);
     322        float FM = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
     323        psSparseVectorElement (sparse, i, FM);
    370324
    371325        // add the per-source variances (border region)
    372326        switch (SKY_FIT_ORDER) {
    373327          case 1:
    374 # if (HAVE_MODEL_VAR)
    375328            f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal);
    376329            psSparseBorderElementB (border, i, 1, f);
    377330            f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal);
    378331            psSparseBorderElementB (border, i, 2, f);
    379 # else
    380             f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    381             psSparseBorderElementB (border, i, 1, f);
    382             f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    383             psSparseBorderElementB (border, i, 2, f);
    384 # endif
    385332
    386333          case 0:
    387 # if (HAVE_MODEL_VAR)
    388334            f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal);
    389 # else
    390             f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    391 # endif
    392335            psSparseBorderElementB (border, i, 0, f);
    393336            break;
    … …  
    409352
    410353            // got an overlap; calculate cross-product and add to output array
    411 # if (HAVE_MODEL_VAR)
    412354            f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal);
    413 # else
    414             f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    415 # endif
    416355            psSparseMatrixElement (sparse, j, i, f);
    417356        }
    … …  
    451390
    452391    // set the sky, sky_x, sky_y components of border matrix
    453 # if (HAVE_MODEL_VAR)
    454392    SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
    455 # else
    456     SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
    457 # endif
    458393
    459394    psSparseConstraint constraint;
    … …  
    500435    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
    501436
     437    // mean stats on fit windows
     438    float sumRadius = 0.0;
     439    float sumPixels = 0.0;
     440    float sumSource = 0.0;
     441
    502442    // measure chisq for each source
    503443    // for (int i = 0; final && (i < fitSources->n); i++) {
    … …  
    505445        pmSource *source = fitSources->data[i];
    506446        pmModel *model = pmSourceGetModel (NULL, source);
     447
     448        // accumulate fit windows statistics
     449        sumRadius += model->fitRadius;
     450        sumPixels += M_PI*PS_SQR(model->fitRadius);
     451        sumSource += 1.0;
     452
    507453        if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) {
    508454            model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value
    … …  
    511457    }
    512458    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
     459
     460    float meanRadius = sumRadius / sumSource;
     461    float meanPixels = sumPixels / sumSource;
    513462
    514463    // psFree (index);
    … …  
    520469    psFree (border);
    521470
    522     psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
     471    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs (mean radius = %f pixels, mean area = %f pixels: %f sec\n", meanRadius, meanPixels, psTimerMark ("psphot.linear"));
    523472
    524473    psphotVisualPlotChisq (sources);
    … …  
    531480    return true;
    532481}
    533 
    534 // XXX do we need this?
    535 // XXX disallow the simultaneous sky fit and remove this code...
    536482
    537483// Calculate the weight terms for the sky fit component of the matrix.  This function operates
    … …  
    614560}
    615561
    616 # if (HAVE_MODEL_VAR)
    617562bool psphotModelBackgroundReadout(psImage *model,  // Model image
    618563                                  psImage *modelStdev, // Model stdev image
    … …  
    624569    );
    625570
    626 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
     571 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
    627572
    628573    bool status = false;
    … …  
    650595    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
    651596
     597    // generate an image of the mean variance image in DN
    652598    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
    653599        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
    … …  
    664610        return false;
    665611    }
    666 
    667612    psFree (varModel);
    668613    psFree (varModelStdev);
    669614
    670     // XXX for a test:
    671     psphotSaveImage (NULL, modelVar, "model.bck.fits");
     615    float gain = 1.0;  // accept 1.0 as a default since it is not critical to the analysis
     616    pmCell *cell = readout->parent; // The parent cell
     617    if (cell) {
     618      gain = psMetadataLookupF32(&status, cell->concepts, "CELL.GAIN"); // Cell gain
     619      if (!status) {
     620        gain = 1.0;           // set note above
     621      }
     622    }
     623    if (gain > 2.0) { /* warn? */ }
     624    // XXX we are not actually using the gain, but need to test it to avoid gcc pedantic warnings
    672625
    673626    // insert all of the source models
    … …  
    690643
    691644        // add the source model to the model variance image
     645        // XXX note that this should be added with gain applied
     646        // var_DN = flux_DN / gain [e/DN]
     647        // to do this requires an API upgrade...
    692648        pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
    693649    }
    694 
    695     // XXX for a test:
    696     psphotSaveImage (NULL, modelVar, "model.var.fits");
    697     psphotSaveImage (NULL, readout->variance, "image.var.fits");
    698650
    699651    // we save the model variance for future reference
    … …  
    724676    return true;
    725677}
    726 # endif

  • trunk/psphot/src/psphotImageLoop.c

    r33691 r34086  
    164164    pmFPAfileActivate (config->files, true, "PSPHOT.EXPNUM");
    165165    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     166        if (! chip->process || ! chip->file_exists) { continue; }
    166167        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed attempting to load EXPNUM input for Chip in psphot.");
    167168        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {

  • trunk/psphot/src/psphotMakeResiduals.c

    r30624 r34086  
    177177
    178178                mflux = 0;
    179                 bool offImage = false;
    180179                if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
    181180                    // This pixel is off the image
    182                     offImage = true;
    183181                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
    184182                } else {

  • trunk/psphot/src/psphotMaskReadout.c

    r29936 r34086  
    9494
    9595    // test output of files at this stage
    96     if (psTraceGetLevel("psphot") >= 5) {
     96    if (psTraceGetLevel("psphot.imsave") >= 5) {
    9797        psphotSaveImage (NULL, readout->image,  "image.fits");
    9898        psphotSaveImage (NULL, readout->mask,   "mask.fits");
    … …  
    105105    return true;
    106106}
     107
     108// XXX this function and support below was created to test the theory that the faint-end
     109// bias results from the Poisson variation of the background pixels.  This is NOT the
     110// case.  Using the code below maintains the faint-end bias.
     111bool psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule) {
     112
     113    bool status = false;
     114
     115    // select the appropriate recipe information
     116    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     117    psAssert (recipe, "missing recipe?");
     118
     119    int num = psphotFileruleCount(config, filerule);
     120
     121    // loop over the available readouts
     122    for (int i = 0; i < num; i++) {
     123
     124        // Generate the mask and weight images, including the user-defined analysis region of interest
     125        if (!psphotUpdateVarianceReadout (config, view, filerule, i, recipe)) {
     126            psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for %s entry %d", filerule, i);
     127            return false;
     128        }
     129    }
     130    return true;
     131}
     132
     133// determine the mean variance image (equivalent to the background model, but for the variance image)
     134// set the variance image to the MAX(input, mean)
     135bool psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
     136
     137    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
     138    psAssert (file, "missing file?");
     139
     140    // find the currently selected readout
     141    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
     142    psAssert (readout, "missing readout?");
     143
     144    pmSourceFitVarMode varMode = psphotGetFitVarMode (recipe);
     145    if (varMode == PM_SOURCE_PHOTFIT_NONE) {
     146      psError (PSPHOT_ERR_CONFIG, false, "need valid LINEAR_FIT_VARIANCE_MODE");
     147      return false;
     148    }
     149
     150    // make this an option via the recipe
     151    if (varMode != PM_SOURCE_PHOTFIT_MODEL_SKY) return true;
     152
     153    // create a model variance image (full-scale image to take result of psImageUnbin below)
     154    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
     155
     156    // find the binning information
     157    psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config);
     158    assert (backBinning);
     159   
     160    psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     161    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     162
     163    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
     164        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
     165        psFree (varModel);
     166        psFree (varModelStdev);
     167        return false;
     168    }
     169
     170    // linear interpolation to full-scale
     171    if (!psImageUnbin (modelVar, varModel, backBinning)) {
     172        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
     173        psFree (varModel);
     174        psFree (varModelStdev);
     175        return false;
     176    }
     177
     178    // XXX save these?
     179    psFree (varModel);
     180    psFree (varModelStdev);
     181
     182    psImage *im = readout->image;
     183    psImage *wt = readout->variance;
     184    for (int j = 0; j < im->numRows; j++) {
     185      for (int i = 0; i < im->numCols; i++) {
     186        if (!isfinite(im->data.F32[j][i])) continue;
     187        if (!isfinite(wt->data.F32[j][i])) continue;
     188        // XXX for a test, make variance constant wt->data.F32[j][i] = PS_MAX(wt->data.F32[j][i], modelVar->data.F32[j][i]);
     189        wt->data.F32[j][i] = modelVar->data.F32[j][i];
     190      }
     191    }
     192
     193    // test output of files at this stage
     194    if (psTraceGetLevel("psphot.imsave") >= 5) {
     195        psphotSaveImage (NULL, readout->image,  "image.varsky.fits");
     196        psphotSaveImage (NULL, readout->mask,   "mask.varsky.fits");
     197        psphotSaveImage (NULL, readout->variance, "variance.varsky.fits");
     198    }
     199
     200    psFree (modelVar);
     201
     202    return true;
     203}

  • trunk/psphot/src/psphotPetrosianStats.c

    r32348 r34086  
    4646
    4747    bool anyPetro = false;
    48     bool manyPetro = false;
     48    // bool manyPetro = false;  XXX not used
    4949    bool above = true;
    5050    float Asum = 0.0;
    … …  
    122122            }
    123123            above = false;
    124             if (anyPetro) manyPetro = true;
     124            // if (anyPetro) manyPetro = true;
    125125            anyPetro = true;
    126126        }
    … …  
    212212    source->extpars->petrosianR50    = R50;
    213213    source->extpars->petrosianR90    = R90;
    214     source->extpars->petrosianFill      = petApix / petArea;
     214    source->extpars->petrosianFill   = petApix / petArea;
    215215   
    216216    // XXX add the errors

  • trunk/psphot/src/psphotReadout.c

    r33915 r34086  
    99}
    1010
     11# if (0)
     12// TEST CODE, can be removed
     13bool psphotDumpFlux (pmConfig *config, const pmFPAview *view, const char *filerule) {
     14
     15    bool status = false;
     16
     17    // find the currently selected readout
     18    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest
     19    psAssert (file, "missing file?");
     20
     21    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     22    psAssert (readout, "missing readout?");
     23
     24    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     25    psAssert (detections, "missing detections?");
     26
     27    psArray *sources = detections->allSources;
     28    psAssert (sources, "missing sources?");
     29
     30    static int npass;
     31    char filename[64];
     32    snprintf (filename, 64, "mags.%d.dat", npass);
     33    FILE *ftest = fopen (filename, "w");
     34    for (int j = 0; j < sources->n; j++) {
     35        pmSource *source = sources->data[j];
     36
     37        float psfMag;
     38        status = pmSourcePhotometryModel (&psfMag, NULL, source->modelPSF);
     39
     40        float psfMagNorm;
     41        float Io = source->modelPSF->params->data.F32[PM_PAR_I0];
     42        source->modelPSF->params->data.F32[PM_PAR_I0] = 1.0;
     43        status = pmSourcePhotometryModel (&psfMagNorm, NULL, source->modelPSF);
     44        source->modelPSF->params->data.F32[PM_PAR_I0] = Io;
     45
     46        // double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y);
     47        fprintf (ftest, "%d %d %d  %f %f %f %f  %f %f\n", j, source->peak->x, source->peak->y, source->modelPSF->params->data.F32[PM_PAR_I0], source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_SXY], psfMag, psfMagNorm);
     48    }
     49    fclose (ftest);
     50    npass++;
     51
     52    return true;
     53}
     54# endif
     55
    1156bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) {
    1257

  • trunk/psphot/src/psphotReadoutCleanup.c

    r29936 r34086  
    77
    88    // remove internal pmFPAfiles, if created
    9     if (psErrorCodeLast() == PSPHOT_ERR_DATA) {
     9    if (psErrorCodeLast() == (psErrorCode) PSPHOT_ERR_DATA) {
    1010        psErrorStackPrint(stderr, "Error in the psphot readout analysis");
    1111        psErrorClear();

  • trunk/psphot/src/psphotSignificanceImage.c

    r31154 r34086  
    88    float SIGMA_SMTH, NSIGMA_SMTH;
    99    bool status = false;
    10     bool guess = false;
    1110
    1211    // smooth the image and variance map
    … …  
    3534        SIGMA_SMTH  = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrt(2.0*log(2.0)));
    3635        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    37         guess = false;
    3836    } else {
    3937        // if we do not know the FWHM, use the guess smoothing kernel supplied.
    … …  
    4341        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    4442        PS_ASSERT (status, NULL);
    45         guess = true;
    4643    }
    4744    // record the actual smoothing sigma

  • trunk/psphot/src/psphotSourceFits.c

    r32744 r34086  
    224224    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    225225    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
    226 
    227 # define TEST_X -420.0
    228 # define TEST_Y 300.0
    229    
    230     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    231         fprintf (stderr, "test galaxy\n");
    232     }
    233 
    234 # undef TEST_X
    235 # undef TEST_Y
    236226
    237227    // set the radius based on the footprint (also sets the mask pixels)
    … …  
    513503    }
    514504
    515 # define TEST_X -540.0
    516 # define TEST_Y 540.0
    517    
    518     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    519         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
    520     }
    521 
    522505    float t1, t2, t4, t5;
    523506    if (TIMING) { psTimerStart ("psphotFitPCM"); }
    … …  
    569552    }
    570553
    571     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    572         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
    573     }
    574 
    575554    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    576555    psFree (pcm);
    … …  
    578557    return model;
    579558}
    580 
    581 # undef TEST_X
    582 # undef TEST_Y
    583559
    584560// note that these should be 1/2n of the standard sersic index
    … …  
    603579    float xMin = NAN;
    604580    float chiSquare[N_INDEX_GUESS];
    605 
    606 # define TEST_X -540.0
    607 # define TEST_Y 540.0
    608    
    609     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    610         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
    611     }
    612581
    613582    for (int i = 0; i < N_INDEX_GUESS; i++) {
    … …  
    635604    assert (iMin >= 0);
    636605
    637     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    638         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
    639     }
    640 
    641606    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    642607    model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
    … …  
    666631    float xMin = NAN;
    667632    float chiSquare[N_INDEX_GUESS];
    668 
    669     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    670         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
    671     }
    672633
    673634    for (int i = 0; i < N_INDEX_GUESS; i++) {
    … …  
    701662    }
    702663    assert (iMin >= 0);
    703    
    704     if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    705         psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
    706     }
    707664
    708665    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it

  • trunk/psphot/src/psphotStackImageLoop.c

  • trunk/psphot/test/tap_psphot_varmodel.pro

    r33963 r34086  
    22# -*-sh-*-
    33
    4 # $KAPA = kapa -noX
    5 
    6 # PSF.CONVOLVE : if true, we insert delta functions (and optionally
    7 #                galaxies) and smooth the image with the psf model
    8 #                (uses a GAUSS regardless of the model). Note that
    9 #                PSF.CONVOLVE = T is faster than F, but (a) only
    10 #                allows Gauss models and (b) only yields quantized
    11 #                locations
    12 
    13 # config for ppImage to generate chip, mask, weight
    14 $ppImageConfig = -recipe PPIMAGE PPIMAGE_N
    15 $ppImageConfig = $ppImageConfig -Db BACKGROUND T
    16 $ppImageConfig = $ppImageConfig -Db CHIP.FITS T
    17 $ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
    18 $ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
    19 $ppImageConfig = $ppImageConfig -Db BASE.FITS F
    20 $ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
    21 $ppImageConfig = $ppImageConfig -Db PHOTOM F
    22 
    23 # basic options for the these images (filter, location, obstype)
    24 $BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0
    25 $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
    26 
    27 # create an image with fake sources and insert the resulting cmf file into a dvodb
    28 $RefConfig = -camera SIMTEST
    29 $RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE
    30 $RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
    31 $RefConfig = $RefConfig -Db PSF.CONVOLVE F
    32 
    33 # options for the reference image
    34 $RefOptions = $BaseOptions -exptime 100.0
    35 $RefOptions = $RefOptions -seeing 1.0
    36 $RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS
    37 $RefOptions = $RefOptions -Df STARS.DENSITY 10.0
    38 $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
    39 
    40 # basic config for ppSim with randomly distributed stars and NO galaxies
    41 $RealConfig = -camera SIMTEST
    42 $RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN
    43 $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
    44 $RealConfig = $RealConfig -Db STARS.FAKE F
    45 $RealConfig = $RealConfig -Db STARS.REAL T
    46 $RealConfig = $RealConfig -Db MATCH.DENSITY F
    47 $RealConfig = $RealConfig -Db PSF.CONVOLVE F
    48 $RealConfig = $RealConfig -Df STARS.DENSITY 10.0
    49 $RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5
    50 $RealConfig = $RealConfig -Db GALAXY.FAKE F
    51 $RealConfig = $RealConfig -Db GALAXY.GRID F
    52 
    53 # options for the repeated images
    54 $RealOptions = $BaseOptions -exptime 30.0
    55  
    56 $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
    57 
    58 # sample alternate options:
    59 # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
    60 # $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
    61 # $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
    62 # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
    63 
    64 list fwhm
    65  1.0
    66  1.1
    67  1.2
    68  1.5
     4# This script includes a set of tests to demonstrate the dependence of the faint-end bias on the weighting scheme
     5# We have 3 weighting schemes :
     6# CONSTANT -- per-pixel weight is fixed (disadvantage: lower S/N, especially for higher sky?)
     7# IMAGE_VAR - per-pixel weight is Poisson from image (disadvantage: faint-end bias)
     8# MODEL_VAR - per-pixel weight is Poisson from model (disadvantage: 2 linear-fit passes)
     9
     10# Functions:
     11#
     12# init  : initialize variables
     13# mkref : generate a fake reference catalog in a DVO database
     14# mkexp : generate a fake exposure from fake catalog and detrend it (saves basename.in.cmf & basename.fits)
     15# runphot : run psphot on an exposure (saves basename.cmf)
     16# ckchip.mags : generate a set of summary plots for the psphot analysis
     17
     18# I would like to produce a grid of tests:
     19# sky (bright, middle, dark)
     20# fwhm (1.0, 1.3, 1.6)
     21# PSF_MODEL input (GAUSS PS1_V1)
     22# PSF_MODEL apply (GAUSS PS1_V1)
     23# variance mode: CONSTANT, IMAGE_VAR, MODEL_VAR
     24
     25macro ppsimtest
     26  data $1
     27  # sf1 : measured flux inserted in image
     28  # sm2 : theoretical magnitude for star
     29  read x 1 y 2 sf1 3 sm2 5 Io 12
     30  set sm1 = -2.5*log(sf1)
     31  set dsm = sm1 - sm2
     32  lim -n 2 sm1 -0.2 0.2; clear; box; plot sm1 dsm
     33end
     34
     35macro go.one
     36  mkdir test
     37
     38  local sky fwhm psf_in psf_out
     39
     40  # mkref creates refimage.* and catdir.ref
     41  file catdir.ref found
     42  if (not($found))
     43    mkref
     44  end
     45
     46  $KAPA = kapa -noX
     47  if (not($?RefConfig)) init
     48
     49  foreach sky 20.0
     50    foreach fwhm 1.0
     51      foreach psf_in PS1_V1
     52
     53        sprintf name "test/test.%02d.%02d.%s" $sky {10*$fwhm} $psf_in
     54        mkexp $name $sky $fwhm $psf_in
     55        data $name.dat
     56        read fluxIn 3 MagPredIn 5
     57        set MagRealIn = -2.5*log(fluxIn)
     58        set dMagIn = MagRealIn - MagPredIn
     59       
     60        foreach psf_out GAUSS PS1_V1
     61          foreach mode CONSTANT IMAGE_VAR MODEL_VAR
     62            runphot $name $name.$psf_out.$mode "-D LINEAR_FIT_VARIANCE_MODE $mode -D PSF_MODEL PS_MODEL_$psf_out"
     63            ckchip.mags $name.in.cmf $name.$psf_out.$mode.cmf $name.$psf_out.$mode 0.0
     64          end
     65        end
     66      end
     67    end
     68  end
    6969end
    7070
    7171macro go
    72 
    7372  mkdir test
    7473
    75   $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
    76   mkexp test/image.00 1.0
    77   $ExtraOptions = -D PSF.MODEL PS_MODEL_PGAUSS
    78   mkexp test/image.01 1.0
    79   $ExtraOptions = -D PSF.MODEL PS_MODEL_PS1_V1
    80   mkexp test/image.02 1.0
     74  local sky fwhm psf_in psf_out
     75
     76  # mkref creates refimage.* and catdir.ref
     77  file catdir.ref found
     78  if (not($found))
     79    mkref
     80  end
     81
     82  $KAPA = kapa -noX
     83  if (not($?RefConfig)) init
     84
     85  foreach sky 19.0 20.0 21.0
     86    foreach fwhm 1.0 1.3 1.6
     87      foreach psf_in GAUSS PS1_V1
     88
     89        sprintf name "test/test.%02d.%02d.%s" $sky {10*$fwhm} $psf_in
     90        # mkexp $name $sky $fwhm $psf_in
     91        data $name.dat
     92        read fluxIn 3 MagPredIn 5
     93        set MagRealIn = -2.5*log(fluxIn)
     94        set dMagIn = MagRealIn - MagPredIn
     95       
     96        foreach psf_out GAUSS PS1_V1
     97          foreach mode CONSTANT IMAGE_VAR MODEL_VAR
     98            # runphot $name $name.$psf_out.$mode "-D LINEAR_FIT_VARIANCE_MODE $mode -D PSF_MODEL PS_MODEL_$psf_out"
     99            ckchip.mags $name.in.cmf $name.$psf_out.$mode.cmf $name.$psf_out.$mode 0.0
     100          end
     101        end
     102      end
     103    end
     104  end
    81105end
    82106
    83107# create a reference database of fake stars to be used by ppSim below
    84108macro mkref
     109  if (not($?RefConfig)) init
     110
    85111  exec rm -rf catdir.ref
    86112  exec rm -f refimage.fits
    … …  
    100126# create a realistic distribution of fake stars, GAUSS PSF
    101127macro mkexp
     128  if ($0 != 5)
     129    echo "USAGE: mkexp basename sky fwhm psf_model"
     130    break
     131  end
     132
     133  local fwhm basename psf_model
     134  $basename = $1
     135  $sky = $2
     136  $fwhm = $3
     137  $psf_model = $4
     138
     139  $ExtraOptions = -D PSF.MODEL PS_MODEL_$psf_model
     140
     141  # create the raw image
     142  echo ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
     143  exec ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
     144  exec /bin/mv -f $basename.cmf $basename.in.cmf
     145
     146  # create the chip output
     147  echo ppImage $ppImageConfig -file $basename.fits $basename
     148  exec ppImage $ppImageConfig -file $basename.fits $basename
     149end
     150
     151# create a realistic distribution of fake stars, GAUSS PSF
     152macro mkexp.deep
    102153  if ($0 != 3)
    103154    echo "USAGE: mkexp basename fwhm"
    … …  
    109160  $fwhm = $2
    110161
     162  $RealOptionsDeep = $BaseOptions -exptime 100
     163
    111164  # create the raw image
    112   echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
    113   exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
     165  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
     166  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
     167  exec /bin/mv -f $basename.cmf $basename.in.cmf
     168
     169  # create the chip output
     170  echo ppImage $ppImageConfig -file $basename.fits $basename
     171  exec ppImage $ppImageConfig -file $basename.fits $basename
     172end
     173
     174# create a realistic distribution of fake stars, GAUSS PSF
     175macro mkexp.bright
     176  if ($0 != 3)
     177    echo "USAGE: mkexp basename fwhm"
     178    break
     179  end
     180
     181  local fwhm basename
     182  $basename = $1
     183  $fwhm = $2
     184
     185  # basic options for the these images (filter, location, obstype)
     186  $BaseOptions = -type OBJECT -filter r -skymags 20.0 -ra 270.70 -dec -23.70 -pa 0.0
     187  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
     188  $RealOptionsDeep = $BaseOptions -exptime 100
     189
     190  # create the raw image
     191  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
     192  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
    114193  exec /bin/mv -f $basename.cmf $basename.in.cmf
    115194
    … …  
    133212  echo psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
    134213  exec psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
     214end
     215
     216# compare two cmf files with extname Chip.psf
     217# things to compare:
     218# * completeness (which sources in (1) are not detected in (2)
     219# * positions (X_PSF, Y_PSF)
     220# * instrumental psf mags
     221# * position errors (no input errors; use a model?)
     222# * measured FWHM?
     223# * kron mags (fluxes)
     224# * etc, etc
     225macro ckchip.mags
     226  if ($0 != 5)
     227    echo "USAGE: ckchip.mags (raw) (out) (output) (zpt_off)"
     228    break
     229  end
     230
     231list pairs
     232  PSF_INST_MAG_out      M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
     233  AP_MAG_out            M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
     234end
     235
     236  load.cmf $1 Chip.psf raw
     237  load.cmf $2 Chip.psf out
     238
     239  # images generated with convolution will not have the right output positions
     240  set X_raw = int(X_PSF_raw) + 0.5
     241  set Y_raw = int(Y_PSF_raw) + 0.5
     242  set M_raw = PSF_INST_MAG_raw + $4
     243  set K_out = -2.5*log(KRON_FLUX_out)
     244  match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2
     245
     246  local i NX NY nx ny N
     247
     248  resize 1000 1000
     249
     250  clear
     251  section a0 0.0 0.0 1.0 0.5
     252  label -fn courier 14;
     253  style -pt 0 -sz 0.4
     254  show.pair 0
     255  reindex ap = AP_MAG_out using index1
     256  set dap = ap - v2
     257  plot -c red -pt 2 -sz 0.5 rv dap
     258
     259  delete -q imag_V dmag_V smag_V Dmag_V dmag_Vraw
     260  for imag -11 -6 0.2
     261    subset dmsub = delta if (rv > $imag) && (rv < $imag + 0.2)
     262    vstat -q dmsub
     263    concat $imag imag_V
     264    concat $MEAN dmag_V
     265    concat $MEDIAN Dmag_V
     266    concat $SIGMA smag_V
     267  end
     268
     269  vstat -q dmag_V
     270  $offset = $MEDIAN
     271  set dmag_V = dmag_V - $offset
     272  set Dmag_V = Dmag_V - $offset
     273
     274  section a1 0.0 0.50 1.0 0.25
     275  lim rv -0.025 0.075; label -fn courier 14; box;
     276  # plot -c black -pt 0 -sz 0.3 MagPredIn dMagIn
     277  plot -c red   -pt 7 -sz 2.0 imag_V dmag_V
     278  plot -c darkgreen -pt 3 -sz 2.0 imag_V Dmag_V
     279  # plot -c gold -pt 4 -sz 2.0 imag_V dmag_Vraw
     280  plot -c blue  -pt 2 -sz 2.0 imag_V smag_V
     281  # label -y "mean (red), median (green), unclipped mean (gold), sigma (blue) of delta"
     282  label -y "mean (red), median (green), sigma (blue) of delta"
     283  sprintf line "OFFSET: %6.3f" $offset
     284  textline -frac 0.1 0.8 -fn courier 24 "$line"
     285
     286  subset dmsub = delta if (rv < -13) && (rv > -16)
     287  vstat -q dmsub
     288  sprintf line "STDEV: %6.3f" $SIGMA
     289  textline -frac 0.1 0.7 -fn courier 24 "$line"
     290
     291
     292  set lChiNorm = log(PSF_CHISQ_out / PSF_NDOF_out)
     293  reindex chi = lChiNorm using index1
     294  section a2 0.0 0.75 1.0 0.25
     295  label -fn courier 14;
     296  lim rv chi; box;
     297  plot -c red -pt 0 -sz 0.5 rv chi
     298
     299  label -y "chisq"
     300  png -name $3.png
     301
     302  # section a1 0.0 0.5 1.0 0.5
     303  # style -pt 0 -sz 0.4
     304  # show.pair 1
     305end
     306
     307macro go.phot
     308 runphot test.00 test.00.varmode.C "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE CONSTANT"
     309 runphot test.00 test.00.varmode.I "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE IMAGE_VAR"
     310 runphot test.00 test.00.varmode.M "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_VAR"
     311 runphot test.00 test.00.varmode.S "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_SKY"
     312end
     313
     314macro go.vars
     315 dev -n varC; ckchip.mags test.00.in.cmf test.00.varmode.C.cmf test.00.varmode.C 0.0
     316 dev -n varI; ckchip.mags test.00.in.cmf test.00.varmode.I.cmf test.00.varmode.I 0.0
     317 dev -n varM; ckchip.mags test.00.in.cmf test.00.varmode.M.cmf test.00.varmode.M 0.0
     318 dev -n varS; ckchip.mags test.00.in.cmf test.00.varmode.S.cmf test.00.varmode.S 0.0
    135319end
    136320
    … …  
    347531    lim rv $word:4 $word:5; box; plot rv delta
    348532  end
    349   label -y '$word:0' -x '$word:2'
     533  $line = '$word:0' - '$word:1'
     534  label -y "$line" -x '$word:2'
    350535end
    351536
    … …  
    637822end
    638823
     824macro init
     825  # config for ppImage to generate chip, mask, weight
     826  $ppImageConfig = -recipe PPIMAGE PPIMAGE_N
     827  $ppImageConfig = $ppImageConfig -Db BACKGROUND T
     828  $ppImageConfig = $ppImageConfig -Db CHIP.FITS T
     829  $ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
     830  $ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
     831  $ppImageConfig = $ppImageConfig -Db BASE.FITS F
     832  $ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
     833  $ppImageConfig = $ppImageConfig -Db PHOTOM F
     834 
     835  # basic options for the these images (filter, location, obstype)
     836  $BaseOptions = -type OBJECT -filter r -ra 270.70 -dec -23.70 -pa 0.0
     837  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
     838 
     839  # PSF.CONVOLVE : if true, we insert delta functions (and optionally
     840  #                galaxies) and smooth the image with the psf model
     841  #                (uses a GAUSS regardless of the model). Note that
     842  #                PSF.CONVOLVE = T is faster than F, but (a) only
     843  #                allows Gauss models and (b) only yields quantized
     844  #                locations
     845
     846  # create an image with fake sources and insert the resulting cmf file into a dvodb
     847  $RefConfig = -camera SIMTEST
     848  $RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE
     849  $RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     850  $RefConfig = $RefConfig -Db PSF.CONVOLVE F
     851 
     852  # options for the reference image
     853  $RefOptions = $BaseOptions
     854  $RefOptions = $RefOptions -exptime 100.0
     855  $RefOptions = $RefOptions -seeing 1.0
     856  $RefOptions = $RefOptions -skymags 21.0 
     857  $RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS
     858  $RefOptions = $RefOptions -Df STARS.DENSITY 10.0
     859  $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
     860
     861  # basic config for ppSim with randomly distributed stars and NO galaxies
     862  $RealConfig = -camera SIMTEST
     863  $RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN
     864  $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     865  $RealConfig = $RealConfig -Db STARS.FAKE F
     866  $RealConfig = $RealConfig -Db STARS.REAL T
     867  $RealConfig = $RealConfig -Db MATCH.DENSITY F
     868  $RealConfig = $RealConfig -Db PSF.CONVOLVE F
     869  $RealConfig = $RealConfig -Df STARS.DENSITY 10.0
     870  $RealConfig = $RealConfig -Df STARS.SIGMA.LIM 1.0
     871  $RealConfig = $RealConfig -Db GALAXY.FAKE F
     872  $RealConfig = $RealConfig -Db GALAXY.GRID F
     873 
     874  # options for the repeated images
     875  $RealOptions = $BaseOptions -exptime 30.0
     876   
     877  # sample alternate options:
     878  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
     879  # $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
     880  # $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
     881  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
     882 
     883  list fwhm
     884   1.0
     885   1.1
     886   1.2
     887   1.5
     888  end
     889end
     890
    639891if ($SCRIPT)
    640892  fulltest 4
Note: See TracChangeset for help on using the changeset viewer.