IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39464


Ignore:
Timestamp:
Mar 14, 2016, 1:40:34 PM (10 years ago)
Author:
eugene
Message:

remove long-deprecated function setMrelCatalog (someday rename setMrelCatalog_alt back); add photFlagsUpper,Lower calculations; save NwarpOK, save median stargal separator

Location:
trunk/Ohana/src/relphot
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/include/relphot.h

    r39345 r39464  
    132132  double *psfqfperf_list;     // psfqfperf for all filters
    133133  double *stargal_list;       // stargal for all filters
     134  uint32_t *photflag_list;      // photflags for all filters
     135
     136  int Nstargal;
     137  int Nphotflags;
    134138
    135139  int   *havePS1;          // this secfilt has synthetic mags

  • trunk/Ohana/src/relphot/src/StarOps.c

    r39457 r39464  
    117117    results->haveUbercal[i] = 0;
    118118
     119    results->Nstargal   = 0;
     120    results->Nphotflags = 0;
     121
    119122    results->tessID[i]    = 0;
    120123    results->projID[i]    = 0;
    … …  
    143146    ALLOCATE (results->psfqfperf_list, double, Nmax);
    144147    ALLOCATE (results->stargal_list,   double, Nmax);
     148    ALLOCATE (results->photflag_list,  uint32_t, Nmax);
    145149
    146150    ALLOCATE (results->havePS1,   int, results->Nsecfilt);
    … …  
    183187    results->psfqfperf_list = NULL;
    184188    results->stargal_list   = NULL;
     189    results->photflag_list  = NULL;
    185190
    186191    results->havePS1 = NULL;
    … …  
    214219  FREE (results->psfqfperf_list);
    215220  FREE (results->stargal_list);
     221  FREE (results->photflag_list);
    216222
    217223  FREE (results->havePS1);
    … …  
    423429    continue; }
    424430
    425 
    426 // setMrel_catalog is used in 3 different contexts:
    427 
    428 // * during the main relphot iterations:
    429 // ** operations only apply to measureT / averageT
    430 // ** we are applying the analysis to the bright subset catalog
    431 // ** we skip any stars found to be bad (STAR_BAD)
    432 
    433 // * during the final pass
    434 // ** operations are applied to all objects
    435 // ** we have special tests for PS1 extended data, 2MASS extended & good data, & synthetic photometry
    436 
    437 // * during the stand-alone '-averages' mode (relphot_objects)
    438 // ** no image data is loaded
    439 // ** we only use the provided calibration (measure.Mcal)
    440 // ** we skip outlier rejection of measurements
    441 
    442 // set the Mrel values for the specified catalog
    443 // NOTE: here 'catalog' is a pointer to a specific catalog, not the root of the array
    444 int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
    445   OHANA_UNUSED_PARAM(catalog);
    446   OHANA_UNUSED_PARAM(Nc);
    447   OHANA_UNUSED_PARAM(pass);
    448   OHANA_UNUSED_PARAM(flatcorr);
    449   OHANA_UNUSED_PARAM(results);
    450   OHANA_UNUSED_PARAM(Nsecfilt);
    451 
    452   fprintf (stderr, "??? should you still be using the old version of setMrel_catalog??\n");
    453   exit (3);
    454 
    455 # if (0)
    456 
    457   off_t j, k, m, ID;
    458   int N;
    459   float Msys, Mcal= 0, Mmos = 0, Mgrid = 0;
    460 
    461   StatType stats, apstats, kronstats;
    462   liststats_setmode (&psfstats, STATMODE);
    463   liststats_setmode (&apstats, STATMODE);
    464   liststats_setmode (&kronstats, STATMODE);
    465 
    466   double *list      = results->list;
    467   double *dlist     = results->dlist;
    468   double *wlist     = results->wlist;
    469   double *aplist    = results->aplist;
    470   double *kronlist  = results->kronlist;
    471   double *dkronlist = results->dkronlist;
    472 
    473   double *psfqf_list     = results->psfqf_list;
    474   double *psfqfperf_list = results->psfqfperf_list;
    475   double *stargal_list   = results->stargal_list;
    476 
    477   SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
    478 
    479   int isSetMrelFinal = (pass >= 0);
    480 
    481   char *primaryCell = NULL;
    482   if (isSetMrelFinal) {
    483     ALLOCATE (primaryCell, char, DVO_MAX_PATH);
    484   }
    485 
    486   for (j = 0; j < catalog[Nc].Naverage; j++) {
    487     // option for a test print
    488     if (FALSE && (catalog[Nc].average[j].objID == 0x7146) && (catalog[Nc].average[j].catID == 0x49d8)) {
    489       fprintf (stderr, "test obj\n");
    490       print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure);
    491     }
    492 
    493     if (isSetMrelFinal) {
    494       // set the name of the primary skycell (this is used in a strcmp to match the skycells in stack detections)
    495       // XXX : this whole function is deprecated
    496       // BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
    497     }
    498 
    499     int GoodPS1 = FALSE;
    500     int Good2MASS = FALSE;
    501     int Galaxy2MASS = FALSE;
    502 
    503     int NextPS1 = 0;
    504     int NpsfPS1 = 0;
    505 
    506     int Nmeas = 0; // number of meas used for psfqf, psfqfperf, and stargal
    507 
    508     int NstackGood = 0;
    509     int NstackSuspect = 0;
    510 
    511     int haveStackObject = FALSE;
    512 
    513     int Ns;
    514     for (Ns = 0; Ns < Nphotcodes; Ns++) {
    515 
    516       int thisCode = photcodes[Ns][0].code;
    517       int Nsec = GetPhotcodeNsec(thisCode);
    518 
    519       /* calculate the average mag in this SEC photcode for a single star */
    520 
    521       /* star/photcodes already calibrated */
    522       if ( isSetMrelFinal && catalog[Nc].found_t[Nsecfilt*j+Nsec]) continue; 
    523      
    524       // skip bad stars
    525       if (!isSetMrelFinal && (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD)) continue;
    526 
    527       int Ncode = 0;
    528       int Next = 0;
    529       int haveSynth = FALSE;
    530       int haveStack = FALSE;
    531       int havePS1   = FALSE;
    532 
    533       // need to find the measurement closest to the center of its skycell, as well as the
    534       // closest for the subset of primary projection cells
    535 
    536       float stackCenterOffsetMin = 1e9;
    537       off_t stackCenterMeasureMin = -1;
    538 
    539       float stackPrimaryOffsetMin = 1e9;
    540       off_t stackPrimaryMeasureMin = -1;
    541 
    542       int forceSynth = FALSE;
    543       int forceSynthEntry = -1;
    544 
    545       int haveUbercal = FALSE;
    546 
    547       int minUbercalDist = 1000;
    548    
    549       int isBad = FALSE;
    550       int isSuspect = FALSE;
    551 
    552       N = 0;
    553       m = catalog[Nc].averageT[j].measureOffset;
    554       for (k = 0; k < catalog[Nc].averageT[j].Nmeasure; k++, m++) {
    555 
    556         // skip measurements that do not match the current photcode
    557         PhotCode *code = GetPhotcodebyCode (catalog[Nc].measureT[m].photcode);
    558         if (!code) continue;
    559         if (code->equiv != thisCode) { continue; }
    560         Ncode ++;
    561 
    562         catalog[Nc].measure[m].dbFlags &= ~ID_MEAS_STACK_PRIMARY;
    563 
    564         if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad);
    565 
    566         if (getImageEntry (m, Nc) < 0) {
    567           // measurements without an image are either external reference photometry or
    568           // data for which the associated image has not been loaded (probably because of
    569           // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
    570           Mmos = Mgrid = 0;
    571           Mcal = catalog[Nc].measureT[m].Mcal; // check that this is zero for loaded REF value
    572         } else {
    573           Mcal  = getMcal  (m, Nc, flatcorr, catalog);
    574           if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
    575           Mmos  = getMmos  (m, Nc);
    576           if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
    577           Mgrid = getMgrid (m, Nc);
    578           if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
    579         }
    580 
    581         // skip some absurd values NAN, < 0.0, > 30.0
    582         Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
    583         if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
    584         if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
    585         if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
    586 
    587         int myUbercalDist = getUbercalDist(m,Nc);
    588         minUbercalDist = MIN(minUbercalDist, myUbercalDist);
    589 
    590         if (isSetMrelFinal) {
    591           // special options for PS1 data
    592           if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) {
    593             // count the extended detections
    594             if (!isnan(catalog[Nc].measure[m].Map)) {
    595               float dMagAp = catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map;
    596               float SigmaAp = hypot(0.1, 2.5*catalog[Nc].measure[m].dM);
    597               // XXX this is still quite ad hoc, but at least it:
    598               // (a) converges to 0.1 mag offset at the bright end
    599               // (b) converges to 0.5 mag offset at the faint end (dM = 0.2)
    600               if (dMagAp > SigmaAp) {
    601                 Next ++;
    602                 NextPS1 ++;
    603               } else {
    604                 NpsfPS1 ++;
    605               }
    606             }
    607             havePS1 = TRUE;
    608           }
    609 
    610           // gpc1 stack data
    611           if ((catalog[Nc].measure[m].photcode >= 11000) && (catalog[Nc].measure[m].photcode <= 11400)) {
    612             haveStack = TRUE;
    613             haveStackObject = TRUE;
    614 
    615             unsigned int stackImageID;
    616 
    617             // which stack image should we use for the mean value?
    618             // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell
    619             if (MatchImageName (m, Nc, primaryCell)) {
    620               float stackPrimaryOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID);
    621               if (stackPrimaryOffset < stackPrimaryOffsetMin) {
    622                 stackPrimaryOffsetMin = stackPrimaryOffset;
    623                 // stackPrimaryIDmin = stackImageID;
    624                 stackPrimaryMeasureMin = m;
    625               }
    626             }
    627 
    628             // get the center distance for the generic case:
    629             float stackCenterOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID);
    630             if (stackCenterOffset < stackCenterOffsetMin) {
    631               stackCenterOffsetMin = stackCenterOffset;
    632               // stackCenterIDmin = stackImageID;
    633               stackCenterMeasureMin = m;
    634             }
    635 
    636             // data quality assessment
    637             isBad |= (catalog[Nc].measure[m].photFlags & code->photomBadMask);
    638             isBad |= (catalog[Nc].measure[m].psfQF < 0.85);
    639             isBad |= isnan(catalog[Nc].measure[m].psfQF);
    640             isBad |= catalog[Nc].measure[m].dM > 0.2; // S/N < 5.0
    641 
    642             isSuspect |= (catalog[Nc].measure[m].photFlags & code->photomPoorMask);
    643             isSuspect |= (catalog[Nc].measure[m].psfQFperf < 0.85);
    644 
    645             // mean psfQF, psfQFperf, and stargal values are calculated per object (all filters, stacks and exposures)
    646             psfqf_list[Nmeas]     = catalog[Nc].measure[m].psfQF;
    647             psfqfperf_list[Nmeas] = catalog[Nc].measure[m].psfQFperf;
    648             stargal_list[Nmeas]   = catalog[Nc].measure[m].extNsigma;
    649             Nmeas++;
    650 
    651             continue; // do not include stack measurements in the per-exposure statistics
    652           }
    653 
    654           // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
    655           if ((catalog[Nc].measure[m].photcode >= 2011) && (catalog[Nc].measure[m].photcode <= 2013)) {
    656             if (catalog[Nc].measure[m].photFlags & 0x00c00000) {
    657               Next ++;
    658               Galaxy2MASS = TRUE;
    659             }
    660             if (pass == 0) {
    661               if (catalog[Nc].measure[m].photFlags & 0x00000007) {
    662                 Good2MASS = TRUE;
    663               } else {
    664                 // detections without one of these bits should only be used in PASS_1
    665                 SKIP_THIS_MEAS(Nbad);
    666               }
    667             }
    668           }
    669 
    670           // Blindly accepth the SYNTH mags if we are above saturation, otherwise,
    671           // ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars)
    672           if ((catalog[Nc].measure[m].photcode >= 3001) && (catalog[Nc].measure[m].photcode <= 3005)) {
    673             // something of a hack: force object to use synth values if synth mags >>
    674             // saturation (3pi instrumental mags < -15)
    675             float MaxMagForceSynth = NAN;
    676             switch (catalog[Nc].measure[m].photcode) {
    677               case 3001:
    678                 MaxMagForceSynth = 13.64;
    679                 break;
    680               case 3002:
    681                 MaxMagForceSynth = 13.76;
    682                 break;
    683               case 3003:
    684                 MaxMagForceSynth = 13.74;
    685                 break;
    686               case 3004:
    687                 MaxMagForceSynth = 12.94;
    688                 break;
    689               case 3005:
    690                 MaxMagForceSynth = 12.01;
    691                 break;
    692             }
    693             if (catalog[Nc].measureT[m].M < MaxMagForceSynth) {
    694               forceSynth = TRUE;
    695               forceSynthEntry = N;
    696             } else {
    697               if (pass < 4) {
    698                 SKIP_THIS_MEAS(Nbad);
    699               }
    700               haveSynth = TRUE;
    701             }
    702           }
    703 
    704           // Map (aplist) and Mkron (kronlist,dkronlist) are used to calculate mean mags per filter
    705           float Map = PhotAper (&catalog[Nc].measure[m]);
    706           aplist[N] = Map - Mcal - Mmos - Mgrid;
    707 
    708           float Mkron = PhotKron (&catalog[Nc].measure[m]);
    709           kronlist[N] = Mkron - Mcal - Mmos - Mgrid;
    710           dkronlist[N] = catalog[Nc].measure[m].dMkron;
    711 
    712           // mean psfQF, psfQFperf, and stargal values are calculated per object (all filters, stacks and exposures)
    713           psfqf_list[Nmeas]     = catalog[Nc].measure[m].psfQF;
    714           psfqfperf_list[Nmeas] = catalog[Nc].measure[m].psfQFperf;
    715           stargal_list[Nmeas]   = catalog[Nc].measure[m].extNsigma;
    716           Nmeas++;
    717         }
    718 
    719         // dlist gives the error per measurement, wlist gives the weight
    720         // we can modify the error and weight in a few ways:
    721         // 1) MIN_ERROR guarantees a floor
    722         // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
    723         // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
    724         // 4) some reference photcode of some kind can be specified as fixed and have a high weight
    725         dlist[N] = MAX (hypot(catalog[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR);
    726         wlist[N] = 1.0;
    727         list[N]  = Msys - Mcal - Mmos - Mgrid;
    728 
    729         // up-weight the ubercal values (or convergence can take a long time...)
    730         if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    731           haveUbercal = TRUE;
    732           wlist[N] = 10.0;
    733         }
    734 
    735         // tie down reference photometry if the -refcode (code) option is selected
    736         // eg, -refcode g_SDSS
    737         // this probably makes no sense in the context of multifilter analysis
    738         if (refPhotcode) {
    739           if (code->code == refPhotcode->code) {
    740             wlist[N] = 100.0;
    741           }
    742         }
    743         N++;
    744       }
    745 
    746       int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
    747 
    748       if (haveStack && (N > 1)) {
    749         // fprintf (stderr, "multiple stack values for %10.6f %10.6f\n", catalog[Nc].averageT[j].R, catalog[Nc].averageT[j].D);
    750       }
    751 
    752       if (haveStack) {
    753         if (!isSuspect && !isBad) {
    754           NstackGood ++;
    755         }
    756         if (isSuspect && !isBad) {
    757           NstackSuspect ++;
    758         }
    759       }
    760 
    761       // when performing the grid analysis, STAR_TOOFEW should be set to 1;
    762       if (N < Nminmeas) { /* too few measurements */
    763         // fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[Nc].averageT[j].R, catalog[Nc].averageT[j].D, catalog[Nc].measureT[catalog[Nc].averageT[j].measureOffset].imageID, N, STAR_TOOFEW);
    764         catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_OBJ_FEW;
    765         if (Ncode == 0) {
    766           results->Ncode ++;
    767         } else {
    768           results->Nfew ++;
    769         }
    770         continue;
    771       } else {
    772         catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_OBJ_FEW;
    773       }
    774 
    775       if (forceSynth) {
    776         // use the single SYNTH value instead of the other mags here
    777         myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error");
    778         list[0]  = list[forceSynthEntry];
    779         dlist[0] = dlist[forceSynthEntry];
    780         wlist[0] = wlist[forceSynthEntry];
    781         N = 1;
    782       }
    783       liststats (list, dlist, wlist, N, &stats);
    784 
    785       catalog[Nc].secfilt[Nsecfilt*j+Nsec].M      = stats.mean;
    786       catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM     = stats.error;
    787       catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mchisq = (stats.Nmeas > 1) ? stats.chisq : NAN;
    788 
    789       // when running -averages, we have no information about the images, so we cannot set this
    790       if (minUbercalDist > -1) {
    791         catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
    792       }
    793 
    794       if (isSetMrelFinal) {
    795         catalog[Nc].found_t[Nsecfilt*j+Nsec] = TRUE;
    796 
    797         catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision)
    798         catalog[Nc].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
    799         catalog[Nc].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;
    800 
    801         catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
    802         catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;
    803 
    804         // NOTE : use the modified weight for apmags as well as psf mags
    805         liststats (aplist, dlist, wlist, N, &apstats);
    806         catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map  = apstats.mean;
    807 
    808         liststats (kronlist, dkronlist, wlist, N, &kronstats);
    809         catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mkron  = kronstats.mean;
    810         catalog[Nc].secfilt[Nsecfilt*j+Nsec].dMkron = kronstats.error;
    811 
    812         if (haveStack) {
    813           m  = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;
    814 
    815           // we are now populating stackDetectID not stack Image ID in secfilt
    816           // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
    817           ID = catalog[Nc].measure[m].extID; // for the stack, this is the stackDetectID
    818 
    819           // get the zero point for the selected image
    820           float zp = PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]) - (Mcal + Mmos + Mgrid);
    821 
    822           // flux_cgs : erg sec^1 cm^-2 Hz^-1
    823           // mag_inst : -2.5 log (cts/sec)
    824           // mag_inst : -2.5 log (flux_inst)
    825           // flux_inst = ten(-0.4*mag_inst)
    826 
    827           // mag_AB = -2.5 log (flux_cgs) - 48.6 (~by definition) [~Vega flux in V-band]
    828           // flux_cgs = ten(-0.4*(mag_AB + 48.6))
    829 
    830           // flux_AB : ten(-0.4*mag_AB)
    831 
    832           // flux_cgs = ten(-0.4*48.6) * flux_AB
    833           // flux_AB  = ten(+0.4*48.6) * flux_cgs
    834            
    835           // flux_Jy : flux_cgs * 10^23
    836 
    837           // flux_AB = ten(+0.4*48.6) * ten(-23) * flux_Jy
    838 
    839           // mag_AB = mag_inst + ZP
    840 
    841           // flux_inst = ten(-0.4*(mag_AB - ZP)) = ten(0.4*ZP) * flux_AB
    842 
    843           // flux_AB = flux_inst * ten(-0.4*ZP)
    844 
    845           // flux_inst * ten(-0.4*ZP) = ten(+0.4*48.6 - 23) * flux_Jy
    846 
    847           // flux_inst = flux_Jy * ten(0.4*ZP + 0.4*48.6 - 23)
    848           // flux_inst = flux_Jy * ten(0.4*ZP - 3.56)
    849           // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56)
    850 
    851           // zpFactor to go from instrumental flux to Janskies
    852           float zpFactor = pow(10.0, -0.4*zp + 3.56);
    853 
    854           // need to put in AB mag factor to get to Janskies (or uJy?)
    855           catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxPSF   = zpFactor * catalog[Nc].measure[m].FluxPSF; 
    856           catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxPSF  = zpFactor * catalog[Nc].measure[m].dFluxPSF;
    857           catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxKron  = zpFactor * catalog[Nc].measure[m].FluxKron;
    858           catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron;
    859 
    860           catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackDetectID = ID;
    861 
    862           catalog[Nc].measure[m].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
    863           if (stackPrimaryMeasureMin >= 0) {
    864             catalog[Nc].measure[m].dbFlags |= ID_MEAS_STACK_PRIMARY;
    865           }
    866         }
    867 
    868         // NOTE: for 2MASS measurements, Next should be 1, as should N
    869         if ((Next > 0) && (Next > 0.5*N)) {
    870           catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT;
    871         }
    872 
    873         switch (pass) {
    874           case 0:
    875             catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
    876             if (havePS1) GoodPS1 = TRUE;
    877             break;
    878           case 1:
    879             catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
    880             if (havePS1) GoodPS1 = TRUE;
    881             break;
    882           case 2:
    883             catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
    884             if (havePS1) GoodPS1 = TRUE;
    885             break;
    886           case 3:
    887             catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3;
    888             break;
    889           case 4:
    890             catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4;
    891             break;
    892         }
    893         if (haveSynth) {
    894           catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH;
    895         }       
    896         if (haveStack) {
    897           catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_HAS_STACK;
    898         }       
    899         if (havePS1) {
    900           catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_HAS_PS1;
    901         }       
    902         if (haveUbercal) {
    903           catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL;
    904         }       
    905       }
    906     }
    907 
    908     if (isSetMrelFinal) {
    909       if (pass == 0) {
    910         DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
    911         // we attempt to set a few flags here; reset those bits before trying:
    912         catalog[Nc].average[j].flags &= ~flagBits;
    913       }
    914 
    915       if (NextPS1 && (NextPS1 > NpsfPS1)) {
    916         catalog[Nc].average[j].flags |= ID_OBJ_EXT;
    917       }
    918       if (GoodPS1) {
    919         catalog[Nc].average[j].flags |= ID_OBJ_GOOD;
    920       }
    921       if (Galaxy2MASS) {
    922         catalog[Nc].average[j].flags |= ID_OBJ_EXT_ALT;
    923       }
    924       if (Good2MASS) {
    925         catalog[Nc].average[j].flags |= ID_OBJ_GOOD_ALT;
    926       }
    927 
    928       if (NstackGood >= 2) {
    929         catalog[Nc].average[j].flags |= ID_OBJ_GOOD_STACK;
    930       } else if (NstackGood + NstackSuspect >= 2) {
    931         catalog[Nc].average[j].flags |= ID_OBJ_SUSPECT_STACK;
    932       } else if (haveStackObject) {
    933         catalog[Nc].average[j].flags |= ID_OBJ_BAD_STACK;
    934       }
    935 
    936       // use NAN values for catalog[Nc].average[j].psfQF as signal
    937       if (!isfinite(catalog[Nc].average[j].psfQF)) {
    938         // do something here
    939         float psfQFmax     = 0.0;
    940         float psfQFperfmax = 0.0;
    941         float stargalmax   = 0.0;
    942         int seq;
    943         for (seq = 0; seq < Nmeas; seq++) {
    944           if (isfinite(psfqf_list[seq])) {
    945             psfQFmax     = MAX (psfqf_list[seq], psfQFmax);
    946           }
    947           if (isfinite(psfqfperf_list[seq])) {
    948             psfQFperfmax = MAX (psfqfperf_list[seq], psfQFperfmax);
    949           }
    950           if (isfinite(stargal_list[seq])) {
    951             stargalmax   = MAX (stargal_list[seq], stargalmax);
    952           }
    953         }
    954 
    955         catalog[Nc].average[j].psfQF     = psfQFmax;
    956         catalog[Nc].average[j].psfQFperf = psfQFperfmax;
    957         catalog[Nc].average[j].stargal   = stargalmax;
    958       }
    959     }
    960   }
    961   if (primaryCell) free (primaryCell);
    962   return (TRUE);
    963 # endif
    964 }
    965431
    966432int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {

  • trunk/Ohana/src/relphot/src/relphot_objects.c

    r39377 r39464  
    7474     
    7575      for (j = 0; j < catalog.Naverage; j++) {
    76         catalog.average[j].flags &= ~photomBits; // reset all except astrometry bits
     76        catalog.average[j].flags    &= ~photomBits; // reset all except astrometry bits
    7777        catalog.average[j].psfQF     = NAN;     // reset (will be re-calculated here)
    7878        catalog.average[j].psfQFperf = NAN;     // reset (will be re-calculated here)
    7979        catalog.average[j].stargal   = NAN;     // reset (will be re-calculated here)
     80        catalog.average[j].photFlagsUpper = 0;  // reset (will be re-calculated here)
     81        catalog.average[j].photFlagsLower = 0;  // reset (will be re-calculated here)
     82        catalog.average[j].NwarpOK        = 0;  // reset (will be re-calculated here)
    8083        for (k = 0; k < Nsecfilt; k++) {
    8184          dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k], SECFILT_RESET_ALL);

  • trunk/Ohana/src/relphot/src/setMrelCatalog.c

    r39457 r39464  
    99
    1010int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags);
     11void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag);
    1112
    1213# define UBERCAL_WEIGHT 100.0
    … …  
    183184        if (isfinite(measure[k].extNsigma) && isfinite(measure[k].psfQF) && (measure[k].psfQF > 0.85)) {
    184185          stargalmax = MAX (stargalmax, measure[k].extNsigma);
     186          results->stargal_list[results->Nstargal] = measure[k].extNsigma;
     187          results->Nstargal ++;
    185188        }
     189
     190        // for all GPC1 measurements, track the photflags
     191        results->photflag_list[results->Nphotflags] = measure[k].photFlags;
     192        results->Nphotflags ++;
    186193      }
    187194
    … …  
    454461      average[0].flags |= ID_OBJ_GOOD_ALT;
    455462    }
     463
     464    GetPhotFlagStats (&average[0].photFlagsUpper, &average[0].photFlagsLower, results->photflag_list, results->Nphotflags);
     465
     466    // calculate the median of the stargal (ExtNsigma) values
     467    dsort(results->stargal_list, results->Nstargal);
     468    double stargal_median;
     469    if (results->Nstargal % 2) {
     470      int Ncenter = results->Nstargal / 2;
     471      stargal_median = results->stargal_list[Ncenter];
     472    } else {
     473      int Ncenter = results->Nstargal / 2 - 1;
     474      stargal_median = 0.5*(results->stargal_list[Ncenter] + results->stargal_list[Ncenter + 1]);
     475    }
     476    average[0].stargal   = stargal_median;
     477
    456478    average[0].psfQF     = psfQfMax;
    457479    average[0].psfQFperf = psfQfPerfMax;
    … …  
    756778  }
    757779
     780  average->NwarpOK = 0; // count the number of warps with psfQF > 0.0
     781
    758782  // assign measurements to the photcode lists
    759783  for (k = 0; k < average->Nmeasure; k++) {
    … …  
    774798
    775799    results->Nmeas[Nsec] ++;
    776     if (measure[k].psfQFperf > 0.85) results->NmeasGood[Nsec] ++;
     800    if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.85)) results->NmeasGood[Nsec] ++;
     801    if (isfinite(measure[k].psfQF)     && (measure[k].psfQF     > 0.00)) average->NwarpOK ++;
    777802
    778803    // use primary skycell for DIFF, stack skycell for WARP
    … …  
    960985}
    961986
     987void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag) {
     988
     989  int i, bit;
     990
     991  uint32_t outputUpper = 0;
     992  uint32_t outputLower = 0;
     993
     994  uint32_t flag = 0x1; // start with the flag on the 0 bit
     995  for (bit = 0; bit < 32; bit++) {
     996
     997    int Nraised = 0;
     998    for (i = 0; i < Nphotflag; i++) {
     999      if (photflag_list[i] & flag) Nraised ++;
     1000    }
     1001
     1002    float fraction = Nraised / (float) Nphotflag;
     1003
     1004    if (fraction < 0.25) {
     1005      // do nothing (no bits raised)
     1006    }
     1007    if ((fraction >= 0.25) && (fraction < 0.50)) {
     1008      outputLower |= flag;
     1009    }
     1010    if ((fraction >= 0.50) && (fraction < 0.75)) {
     1011      outputUpper |= flag;
     1012    }
     1013    if (fraction >= 0.75) {
     1014      outputLower |= flag;
     1015      outputUpper |= flag;
     1016    }
     1017    flag = (flag << 1);
     1018  }
     1019  *photFlagUpper = outputUpper;
     1020  *photFlagLower = outputLower;
     1021}

  • trunk/Ohana/src/relphot/src/setMrelFinal.c

    r39377 r39464  
    9494    catalog[0].average[i].psfQFperf = NAN;      // force recalculation below
    9595    catalog[0].average[i].stargal   = NAN;      // force recalculation below
     96    catalog[0].average[i].photFlagsUpper = 0;   // reset (will be re-calculated here)
     97    catalog[0].average[i].photFlagsLower = 0;   // reset (will be re-calculated here)
     98    catalog[0].average[i].NwarpOK        = 0;   // reset (will be re-calculated here)
    9699  }
    97100
Note: See TracChangeset for help on using the changeset viewer.