IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 30, 2013, 5:42:28 AM (13 years ago)
Author:
eugene
Message:

adding CFF input for psphotFullForce; adding XGAL output element for psphotFullForce

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CMF.c.in

    r36158 r36162  
    13441344}
    13451345
    1346 // XXX this layout is still the same as PS1_DEV_1
    1347 bool pmSourcesWrite_CMF_@CMFMODE@_GAL (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
     1346// XXX where should I record the number of columns??
     1347bool pmSourcesWrite_CMF_@CMFMODE@_XGAL (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
    13481348{
    1349 
    1350     psArray *table;
    1351     psMetadata *row;
    1352     psF32 *PAR, *dPAR;
    1353     psF32 xPos, yPos;
    1354     psF32 xErr, yErr;
    1355     char name[64];
     1349    bool status = false;
     1350
     1351    // perform full non-linear fits / extended source analysis?
     1352    if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) {
     1353        psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n");
     1354        return true;
     1355    }
    13561356
    13571357    // create a header to hold the output data
    … …  
    13641364    sources = psArraySort (sources, pmSourceSortByFlux);
    13651365
    1366     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
    1367     int nParamMax = 0;
     1366    psArray *table = psArrayAllocEmpty (sources->n);
     1367
    13681368    for (int i = 0; i < sources->n; i++) {
    1369         // this is the source associated with this image
     1369
    13701370        pmSource *thisSource = sources->data[i];
    13711371
    … …  
    13731373        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    13741374
     1375        // if we did not fit the galaxy model, modelFits will be NULL
    13751376        if (source->modelFits == NULL) continue;
    1376         for (int j = 0; j < source->modelFits->n; j++) {
    1377             pmModel *model = source->modelFits->data[j];
    1378             assert (model);
    1379             nParamMax = PS_MAX (nParamMax, model->params->n);
    1380         }
    1381     }
    1382 
    1383     @>PS1_DV2@ pmChip *chip = readout->parent->parent;
    1384 
    1385     table = psArrayAllocEmpty (sources->n);
    1386 
    1387     // we write out all sources, regardless of quality.  the source flags tell us the state
    1388     for (int i = 0; i < sources->n; i++) {
    1389 
    1390         pmSource *thisSource = sources->data[i];
    1391 
    1392         // this is the "real" version of this source
    1393         pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    1394 
    1395         // XXX if no model fits are saved, write out modelEXT?
    1396         if (source->modelFits == NULL) continue;
    1397 
    1398         // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
    1399         for (int j = 0; j < source->modelFits->n; j++) {
    1400 
    1401             // choose the convolved EXT model, if available, otherwise the simple one
    1402             pmModel *model = source->modelFits->data[j];
    1403             assert (model);
    1404 
    1405             // pmSourceExtFitPars *extPars = source->extFitPars->data[j];
    1406             // assert (extPars);
    1407 
    1408             // skip models which were not actually fitted
    1409             if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
    1410 
    1411             PAR = model->params->data.F32;
    1412             dPAR = model->dparams->data.F32;
    1413             xPos = PAR[PM_PAR_XPOS];
    1414             yPos = PAR[PM_PAR_YPOS];
    1415             xErr = dPAR[PM_PAR_XPOS];
    1416             yErr = dPAR[PM_PAR_YPOS];
    1417 
    1418             @>PS1_DV2@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    1419             @>PS1_DV2@ float posAngle = 0.0;
    1420             @>PS1_DV2@ float pltScale = 0.0;
    1421             @>PS1_DV2@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
    1422             @>PS1_DV2@ double raPos = ptSky.r*PS_DEG_RAD;
    1423             @>PS1_DV2@ double decPos = ptSky.d*PS_DEG_RAD;
    1424             @>PS1_DV2@ posAngle *= PS_DEG_RAD;
    1425             @>PS1_DV2@ pltScale *= PS_DEG_RAD*3600.0;
    1426 
    1427             float kronFlux = source->moments ? source->moments->KronFlux : NAN;
    1428             float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN;
    1429 
    1430             row = psMetadataAlloc ();
    1431 
    1432             // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
    1433             // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
    1434             // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment"
    1435             psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
    1436             psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
    1437             psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
    1438             psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
    1439             psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
    1440             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT",           0, "EXT model ra coordinate",                    raPos);
    1441             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT",          0, "EXT model dec coordinate",                   decPos);
    1442             @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN;
    1443             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX",    0, "EXT fit instrumental counts",                instFlux);
    1444 
    1445             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
    1446             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
    1447 
    1448             @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN;
    1449             @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration",   calMag);
    1450             @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Magnitude using supplied calibration",   model->chisq);
    1451             @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Magnitude using supplied calibration",   model->nDOF);
    1452 
    1453             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE",   0, "position angle at source (degrees)",         posAngle);
    1454             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE",   0, "plate scale at source (arcsec/pixel)",       pltScale);
    1455 
    1456             // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XX",       0, "second moment in x",                      extPars->Mxx);
    1457             // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XY",       0, "second moment in x,y",                    extPars->Mxy);
    1458             // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_YY",       0, "second moment in y",                      extPars->Myy);
    1459             // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_R1",       0, "first radial moment",                     extPars->Mrf);
    1460             // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_RH",       0, "half radial moment",                      extPars->Mrh);
    1461 
    1462             psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_INST_MAG",     0, "PSF fit instrumental magnitude",             source->psfMag);
    1463             psMetadataAddF32 (row, PS_LIST_TAIL, "AP_MAG",           0, "PSF-sized aperture magnitude",               source->apMag);
    1464             psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_MAG",         0, "Kron Mag",                                   kronMag);
    1465 
    1466             @PS1_DV0,PS1_DV1@ psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
    1467             @ALL,!PS1_DV0,!PS1_DV1@ psMetadataAddS32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
    1468             psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
    1469 
    1470             // XXX these should be major and minor, not 'x' and 'y'
    1471             if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) {
    1472                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]);
    1473                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  PAR[PM_PAR_SIGMA]); // this is not fitted
    1474                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    PAR[PM_PAR_THETA]);
    1475                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_LENGTH]);
    1476                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            NAN); // this is not fitted, so error is NAN
    1477                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_THETA]);
    1478             } else {
    1479                 if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY])) {
    1480                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",     0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]);
    1481                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",     0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]);
    1482                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",         0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]);
    1483                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]);
    1484                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]);
    1485                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",     0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]);
    1486                 } else {
    1487                     psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
    1488                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
    1489                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);
    1490                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    axes.theta);
    1491                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_SXX]);
    1492                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            dPAR[PM_PAR_SYY]);
    1493                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_SXY]);
    1494                 }
    1495             }
    1496 
    1497             // write out the other generic parameters
    1498             for (int k = 0; k < nParamMax; k++) {
    1499                 if (k == PM_PAR_I0) continue;
    1500                 if (k == PM_PAR_SKY) continue;
    1501                 if (k == PM_PAR_XPOS) continue;
    1502                 if (k == PM_PAR_YPOS) continue;
    1503                 if (k == PM_PAR_SXX) continue;
    1504                 if (k == PM_PAR_SXY) continue;
    1505                 if (k == PM_PAR_SYY) continue;
    1506 
    1507                 snprintf (name, 64, "EXT_PAR_%02d", k);
    1508 
    1509                 if (k < model->params->n) {
    1510                     psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
    1511                 } else {
    1512                     psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN);
    1513                 }
    1514             }
    1515 
    1516             // optionally, write out the covariance matrix values
    1517             // XXX do I need to pad this to match the biggest covar matrix?
    1518             if (model->covar) {
    1519                 for (int iy = 0; iy < model->covar->numCols; iy++) {
    1520                     for (int ix = iy; ix < model->covar->numCols; ix++) {
    1521                         snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
    1522                         psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);
    1523 
    1524                     }
    1525                 }                   
    1526             }
     1377
     1378        // if we did not fit the galaxy model, galaxyFits will also be NULL
     1379        if (source->galaxyFits == NULL) continue;
     1380
     1381        pmModel *model = source->modelFits->data[0];
     1382        if (!model) return false;
     1383
     1384        // X,Y coordinates are stored with the model parameters
     1385        psF32 *PAR = model->params->data.F32;
     1386
     1387        psMetadata *row = psMetadataAlloc ();
     1388
     1389        // we write out the x,y positions so people can link to the psf either way (position or ID)
     1390        psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index", source->seq);
     1391        psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT",            0, "model x coordinate",             PAR[PM_PAR_XPOS]);
     1392        psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT",            0, "model y coordinate",             PAR[PM_PAR_YPOS]);
     1393        psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX",             0, "number of pixels for fits",      source->galaxyFits->nPix);
     1394
     1395        psVector *Io = source->galaxyFits->Io;
     1396        psVector *dIo = source->galaxyFits->dIo;
     1397        psVector *chisq = source->galaxyFits->chisq;
     1398
     1399        // we write a series of columns for each of the elements (maybe I should save them as vectors instead?)
     1400        for (int j = 0; j < Io->n; j++) {
     1401
     1402            char field[16];
     1403
     1404            snprintf (field, 16, "CHI_%02d", j);
     1405            psMetadataAddF32 (row, PS_LIST_TAIL, field, 0, "chisq of model fit", chisq->data.F32[j]);
     1406
     1407            snprintf (field, 16, "I0_%02d", j);
     1408            psMetadataAddF32 (row, PS_LIST_TAIL, field, 0, "chisq of model fit", Io->data.F32[j]);
     1409
     1410            snprintf (field, 16, "D_I0_%02d", j);
     1411            psMetadataAddF32 (row, PS_LIST_TAIL, field, 0, "chisq of model fit", dIo->data.F32[j]);
     1412
    15271413            psArrayAdd (table, 100, row);
    15281414            psFree (row);
    … …  
    15421428    }
    15431429
    1544     psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     1430    psTrace ("pmFPAfile", 5, "writing galaxy data %s\n", extname);
    15451431    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    1546         psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
     1432        psError(psErrorCodeLast(), false, "writing galaxy data %s\n", extname);
    15471433        psFree (outhead);
    15481434        psFree(table);
Note: See TracChangeset for help on using the changeset viewer.