IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 13, 2017, 10:53:48 AM (9 years ago)
Author:
eugene
Message:

merge EAM development branch changes for DR2 into trunk (add PS1_V6 dvo format; change Mcal to McalPSF, McalAPER; change opihi int vectors to 64bit)

Location:
trunk/Ohana
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

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

    r39926 r40291  
    22
    33# if (0)
    4 # define TEST_OBJ_ID 0x0000000e
    5 # define TEST_CAT_ID 0x000076ee
     4# define TEST_OBJ_ID 0x000098f9
     5# define TEST_CAT_ID 0x00005ad3
    66# else
    77# define TEST_OBJ_ID 0
    … …  
    102102// 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
    103103// 4) some reference photcode of some kind can be specified as fixed and have a high weight
     104
     105// Although I calculate McalAPER for exposures, I am only using McalPSF for chips.  Note
     106// in StarOps.c:setMcalOutput I am setting measure->McalAPER to image->McalPSF for chips
     107// and warps, but not stacks
    104108
    105109// set mean of chip measurements (selected by photcode range for now):
    … …  
    224228      // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
    225229      Mmos = Mgrid = 0;
    226       Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value
     230      Mcal = measureT[k].McalPSF; // check that this is zero for loaded REF value
    227231    } else {
    228       // getMcal_alt returns image[].Mcal modified by flatcorr(image.photom_map_id,x,y)
    229       // NOTE: getMcal_alt does not include measure.Mflat
    230       Mcal  = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd);
     232      // getMcal returns image[].Mcal; note the flat-field correction is stored in measure.Mflat
     233      Mcal  = getMcal (meas, cat, MAG_CLASS_PSF);
    231234      if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
    232235      Mmos  = getMmos  (meas, cat);
    … …  
    268271    }
    269272
    270     float Mkron, dMkron;
     273    // NOTE: we need to calculate the averge chip Kron on each pass to be able to calibrate the stacks
     274    float Mkron = PhotCatTiny (&measureT[k], MAG_CLASS_KRON);
     275    float dMkron;
    271276    if (isSetMrelFinal) {
    272       Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
    273277      dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR);
    274       if (CHECK_VALID_MAG(Mkron, dMkron)) {
    275         int Nkron = results->kronData[Nsec].Nlist;
    276         results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
    277         results->kronData[Nsec].errlist[Nkron] = dMkron;
    278         results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    279         results->kronData[Nsec].ranking[Nkron] = measureRank[k];
    280         results->kronData[Nsec].measSeq[Nkron] = k;
    281         results->kronData[Nsec].msklist[Nkron] = 0;
    282         results->kronData[Nsec].Nlist ++;
    283       }
     278    } else {
     279      dMkron = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
     280    }
     281    if (CHECK_VALID_MAG(Mkron, dMkron)) {
     282      int Nkron = results->kronData[Nsec].Nlist;
     283      results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
     284      results->kronData[Nsec].errlist[Nkron] = dMkron;
     285      results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
     286      results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
     287      results->kronData[Nsec].measSeq[Nkron] = k;
     288      results->kronData[Nsec].msklist[Nkron] = 0;
     289      results->kronData[Nsec].Nlist ++;
    284290    }
    285291
    … …  
    377383        }
    378384      }
    379       secfilt[Nsec].M      = Mpsf + ZP;
    380       secfilt[Nsec].dM    = 0.6;
    381       secfilt[Nsec].Mchisq = 0.0;
    382       secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
     385      secfilt[Nsec].MpsfChp  = Mpsf + ZP;
     386      secfilt[Nsec].dMpsfChp = 0.6;
     387      secfilt[Nsec].Mchisq   = 0.0;
     388      secfilt[Nsec].flags   |= ID_SECF_USE_SYNTH;
    383389      continue;
    384390    }
    … …  
    386392    // if too few valid measurements meet the minimum criteria, go to the next entry
    387393    StatType *psfstats = &results->psfstats;
    388     int Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
    389     if (Nranking < Nminmeas) {
     394    int NrankingPSF = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
     395    if (NrankingPSF < Nminmeas) {
    390396      secfilt[Nsec].flags |= ID_OBJ_FEW;
    391397    } else {
    392       secfilt[Nsec].M    = psfstats->mean;
    393       secfilt[Nsec].dM   = psfstats->error;
    394       secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
    395     }
    396     int minRankPSF = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
     398      secfilt[Nsec].MpsfChp  = psfstats->mean;
     399      secfilt[Nsec].dMpsfChp = psfstats->error;
     400      secfilt[Nsec].Mchisq   = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
     401    }
    397402
    398403    // when running -averages, we have no information about the images, so we cannot set this
    … …  
    401406    }
    402407
     408    StatType *kronstats = &results->kronstats;
     409    int NrankingKRON = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
     410    if (NrankingKRON) {
     411      secfilt[Nsec].MkronChp  = kronstats->mean;
     412      secfilt[Nsec].dMkronChp = kronstats->error;
     413      secfilt[Nsec].sMkronChp = kronstats->sigma;
     414      secfilt[Nsec].NusedKron = NrankingKRON;
     415    }
     416
    403417    if (isSetMrelFinal) {
    404418      if ((average[0].objID == TEST_OBJ_ID) && (average[0].catID == TEST_CAT_ID)) {
    … …  
    406420      }
    407421
    408       // mark the measurements matching this ranking
    409       markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
    410 
    411       if (Nranking) {
    412         secfilt[Nsec].Mstdev = psfstats->sigma;
    413         secfilt[Nsec].Nused  = psfstats->Nmeas;
    414         secfilt[Nsec].Mmax   = psfstats->max;
    415         secfilt[Nsec].Mmin   = psfstats->min;
     422      if (NrankingPSF) {
     423        secfilt[Nsec].sMpsfChp = psfstats->sigma;
     424        secfilt[Nsec].Nused    = psfstats->Nmeas;
     425        secfilt[Nsec].Mmax     = psfstats->max;
     426        secfilt[Nsec].Mmin     = psfstats->min;
    416427      }
    417428      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
    … …  
    421432      MAX_NOT_NAN (psfQfPerfMax, secfilt[Nsec].psfQfPerfMax);
    422433
     434      // mark the measurements matching this ranking
     435      int minRankPSF = (NrankingPSF > 0) ? results->psfData[Nsec].ranking[0] : 10;
     436      markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
     437
     438      int minRankKron = (NrankingKRON > 0) ? results->kronData[Nsec].ranking[0] : 10;
     439      markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
     440
    423441      StatType *apstats = &results->apstats;
    424       Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
    425       if (Nranking) {
    426         secfilt[Nsec].Map     = apstats->mean;
    427         secfilt[Nsec].dMap    = apstats->error;
    428         secfilt[Nsec].sMap    = apstats->sigma;
    429         secfilt[Nsec].NusedAp = Nranking;
    430       }
    431       int minRankAper = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10;
     442      int NrankingAPER = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
     443      if (NrankingAPER) {
     444        secfilt[Nsec].MapChp  = apstats->mean;
     445        secfilt[Nsec].dMapChp = apstats->error;
     446        secfilt[Nsec].sMapChp = apstats->sigma;
     447        secfilt[Nsec].NusedAp = NrankingAPER;
     448      }
     449      int minRankAper = (NrankingAPER > 0) ? results->aperData[Nsec].ranking[0] : 10;
    432450      markMeasureByRanking (&results->aperData[Nsec], measure, minRankAper, ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER);
    433 
    434       StatType *kronstats = &results->kronstats;
    435       Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
    436       if (Nranking) {
    437         secfilt[Nsec].Mkron     = kronstats->mean;
    438         secfilt[Nsec].dMkron    = kronstats->error;
    439         secfilt[Nsec].sMkron    = kronstats->sigma;
    440         secfilt[Nsec].NusedKron = Nranking;
    441       }
    442       int minRankKron = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10;
    443       markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
    444451
    445452      // does this object appear extended in > 50% of measurements?
    … …  
    518525// 2) select the BEST detections per filter (regardless of PRIMARY)
    519526// 3) apply the zero point and AB->Jy transformations
     527
     528// I calculate McalAPER and McalPSF independently for stacks.  I use McalAPER for Mkron
     529// and Map, and McalPSF for Mpsf.  Note in StarOps.c:setMcalOutput I am setting
     530// measure->McalAPER to image->McalPSF for chips and warps, but not stacks
     531
    520532int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr) {
    521533
    … …  
    528540  SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt];
    529541
    530   off_t k;
    531 
    532   float Mcal = 0, Mmos = 0, Mgrid = 0, Finst = 0;
     542  float McalPSF = 0, McalAPER = 0, Mmos = 0, Mgrid = 0;
    533543
    534544  // set the primary projection cell and skycell for this coordinate
    … …  
    546556  int haveStackObject = FALSE;
    547557
    548   int Ns;
    549   for (Ns = 0; Ns < Nphotcodes; Ns++) {
     558  for (int Ns = 0; Ns < Nphotcodes; Ns++) {
    550559
    551560    int thisCode = photcodes[Ns][0].code;
    … …  
    556565    int haveStack = FALSE;
    557566
    558     float psfQFbest = 0.0;
    559 
    560     off_t stackBestMeasure = -1;
    561     off_t stackPrimaryMeasure = -1;
    562 
    563567    int isBad = FALSE;
    564568    int isSuspect = FALSE;
    565569
    566     int Nstack = 0; // number for this photcode
    567     int NstackDet = 0; // number for this photcode
    568 
     570    int NstackMeas = 0; // number of stack measurements for this photcode
     571    int NstackDet  = 0; // number of stack detections for this photcode (not forced)
     572
     573    int   primaryLevelMax = 0;
     574    off_t primaryEntryMax = -1;
     575    float primaryValueMax = 0.0;
     576
     577    int   bestLevelMax = 0;
     578    off_t bestEntryMax = -1;
     579    float bestValueMax = 0.0;
     580
     581    int Nprimary = 0;
     582
     583    // reset all stack-related values for this secfilt:
    569584    secfilt[Nsec].stackBestOff = -1;
    570585    secfilt[Nsec].stackPrmryOff = -1;
    571 
    572     off_t meas = measureOffset;
    573     for (k = 0; k < Nmeasure; k++, meas++) {
     586    secfilt[Nsec].FpsfStk   = NAN;
     587    secfilt[Nsec].dFpsfStk  = NAN;
     588    secfilt[Nsec].FkronStk  = NAN;
     589    secfilt[Nsec].dFkronStk = NAN;
     590    secfilt[Nsec].FapStk    = NAN;
     591    secfilt[Nsec].dFapStk   = NAN;
     592    secfilt[Nsec].MpsfStk   = NAN;
     593    secfilt[Nsec].MkronStk  = NAN;
     594    secfilt[Nsec].MapStk    = NAN;
     595    secfilt[Nsec].Nstack    = 0;
     596    secfilt[Nsec].NstackDet = 0;
     597    secfilt[Nsec].flags    &= ~ID_SECF_STACK_FLAGS; // reset the stack flags
     598
     599    for (off_t k = 0; k < Nmeasure; k++) {
    574600
    575601      // only examine gpc1 stack data
    … …  
    581607      if (code->equiv != thisCode) { continue; }
    582608
    583       Nstack ++;
    584       if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++;
     609      NstackMeas ++;
     610      if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++; // 0x04 = PM_SOURCE_MODE2_MATCHED
    585611
    586612      // clear this bit for all measurements
    … …  
    590616      haveStack = TRUE;
    591617      haveStackObject = TRUE;
     618
     619      // if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad);
     620
     621      int isPrimary = FALSE;
    592622     
    593       // ** find the PRIMARY measurement
    594 
    595       // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement is from the
    596       // primary skycell for this position
    597       if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) {
    598         stackPrimaryMeasure = k;
     623      // ** is this a PRIMARY measurement? (there may be more than one)
     624
     625      // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement
     626      // is from the primary skycell for this position. (note that MatchImageSkycellID
     627      // requires the entry in the full Measure table, not just this object)
     628      if (MatchImageSkycellID (measureOffset + k, cat, tessID, projectID, skycellID)) {
    599629        measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
    600         secfilt[Nsec].stackPrmryOff = meas;
    601         myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range");
    602       }
    603 
    604       // ** now choose the BEST measurements (may also be PRIMARY)
    605 
    606       // ensure that we at least have a single best measure
    607       if (stackBestMeasure == -1) stackBestMeasure = k;
    608 
    609       // choose the best psfQFperf value for the BEST measurement
    610       if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > psfQFbest)) {
    611         psfQFbest = measure[k].psfQFperf;
    612         stackBestMeasure = k;
    613       }
    614       // ... UNLESS psfQFperf > 0.98 for the primary, in which case just use the primary.
    615       if ((measure[k].dbFlags & ID_MEAS_STACK_PRIMARY) && isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.98)) {
    616         psfQFbest = 1000; // force this to be the best entry
    617         stackBestMeasure = k;
    618       }
    619 
    620       if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad);
    621 
     630        isPrimary = TRUE;
     631        Nprimary ++;
     632      }
     633
     634      // soften the error floor (can dM be 0.0?)
     635      int SNvalue = isfinite(measure[k].dM) ? 1.0 / hypot (measure[k].dM, MIN_ERROR) : NAN;
     636      int psfQFperfAboveLimit = isfinite(SNvalue) && isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.95);
     637
     638      // ** determine the BEST level
     639      int bestLevel = 0;
     640      if (!isPrimary && !psfQFperfAboveLimit) bestLevel = 1;
     641      if ( isPrimary && !psfQFperfAboveLimit) bestLevel = 2;
     642      if (!isPrimary &&  psfQFperfAboveLimit) bestLevel = 3;
     643      if ( isPrimary &&  psfQFperfAboveLimit) bestLevel = 4;
     644
     645      // here is the rule for choosing the best value:
     646      float bestValue = (bestLevel > 2) ? SNvalue : measure[k].psfQFperf;
     647
     648      // if we have not reached this level before, set the new level
     649      if (bestLevelMax < bestLevel) {
     650        bestValueMax = bestValue;
     651        bestLevelMax = bestLevel;
     652        bestEntryMax = k;
     653      }
     654      // if we have reached this level before, set the new value if we beat the old one
     655      if ((bestLevelMax == bestLevel) && (bestValueMax < bestValue)) {
     656        bestValueMax = bestValue;
     657        bestEntryMax = k;
     658      }
     659      myAssert (bestEntryMax > -1, "this should not happen");
     660
     661      // ** determine the PRIMARY level
     662      int primaryLevel = 0;
     663      if ( isPrimary && !psfQFperfAboveLimit) primaryLevel = 1;
     664      if ( isPrimary &&  psfQFperfAboveLimit) primaryLevel = 2;
     665
     666      // here is the rule for choosing the PRIMARY value:
     667      float primaryValue = 0.0;
     668      if (primaryLevel == 1) { primaryValue = measure[k].psfQFperf; }
     669      if (primaryLevel == 2) { primaryValue = SNvalue; }
     670
     671      // if we have not reached this level before, set the new level
     672      if (isPrimary && (primaryLevelMax < primaryLevel)) {
     673        primaryValueMax = primaryValue;
     674        primaryLevelMax = primaryLevel;
     675        primaryEntryMax = k;
     676      }
     677      // if we have reached this level before, set the new value if we beat the old one
     678      if (isPrimary && (primaryLevelMax == primaryLevel) && (primaryValueMax < primaryValue)) {
     679        primaryValueMax = primaryValue;
     680        primaryEntryMax = k;
     681      }
     682    }
     683
     684    // if we do not have a stack measurement for this photcode, skip everything below
     685    if (!haveStack) continue;
     686 
     687    // now we have a BEST and a PRIMARY entry (may be the same entry)
     688
     689    off_t meas = bestEntryMax;
     690    off_t measSeq = meas + measureOffset;
     691
     692    myAssert (meas >= 0, "this should not happen");
     693
     694    // measurements without an image are either external reference photometry or
     695    // data for which the associated image has not been loaded (probably because of
     696    // overlaps).  we only want measurements associated with stack images in this loop
     697   
     698    // match measurement to its image (this is just a check, right?)
     699    if (getImageEntry (measSeq, cat) < 0) {
    622700      // measurements without an image are either external reference photometry or
    623701      // data for which the associated image has not been loaded (probably because of
    624       // overlaps).  we only want measurements associated with stack images in this loop
    625 
    626       // match measurement to its image (this is just a check, right?)
    627       if (getImageEntry (meas, cat) < 0) {
    628         // measurements without an image are either external reference photometry or
    629         // data for which the associated image has not been loaded (probably because of
    630         // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
    631         Mmos = Mgrid = 0;
    632         Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
    633       } else {
    634         Mcal  = getMcal_alt  (meas, cat, flatcorr, measure[k].Xccd, measure[k].Yccd);
    635         if (isnan(Mcal))  SKIP_THIS_MEAS_STACK(Ncal);
    636         Mmos  = getMmos  (meas, cat);
    637         if (isnan(Mmos))  SKIP_THIS_MEAS_STACK(Nmos);
    638         Mgrid = getMgrid (meas, cat);
    639         if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid);
    640       }
    641 
    642       // NOTE: negative and insignificant flux values are allowed, but not NAN flux values
    643       Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF);
    644       if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst);
    645 
    646       // data quality assessment
    647       isBad |= (measure[k].photFlags & code->photomBadMask);
    648       isBad |= (measure[k].psfQF < 0.85);
    649       isBad |= isnan(measure[k].psfQF);
    650       isBad |= measure[k].dM > 0.2; // S/N < 5.0
    651 
    652       isSuspect |= (measure[k].photFlags & code->photomPoorMask);
    653       isSuspect |= (measure[k].psfQFperf < 0.85);
    654     }
    655 
    656     if (!haveStack) continue;
    657 
    658     if (!isSuspect && !isBad) {
    659       NstackGood ++;
    660     }
    661     if (isSuspect && !isBad) {
    662       NstackSuspect ++;
    663     }
    664 
    665     // measurements which are bad will not have a valid stack entry and are skipped
    666     k = (stackBestMeasure >= 0) ? stackBestMeasure : stackPrimaryMeasure;
    667     if (k < 0) continue;
    668 
    669     // we are now populating stackDetectID not stack Image ID in secfilt
    670     // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
    671     // ID = measure[k].extID; // for the stack, this is the stackDetectID
    672 
    673     // get the zero point for the selected image
    674     float zp = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (Mcal + Mmos + Mgrid);
     702      // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
     703      McalPSF  = measure[meas].McalPSF; // check that this is zero for loaded REF value
     704      McalAPER = McalPSF; // check that this is zero for loaded REF value
     705      Mmos     = 0.0;
     706      Mgrid    = 0.0;
     707    } else {
     708      McalPSF   = getMcal  (measSeq, cat, MAG_CLASS_PSF);
     709      McalAPER  = USE_MCAL_PSF_FOR_STACK_APER ? getMcal  (measSeq, cat, MAG_CLASS_PSF) : getMcal  (measSeq, cat, MAG_CLASS_KRON);
     710      Mmos      = getMmos  (measSeq, cat);
     711      Mgrid     = getMgrid (measSeq, cat);
     712    }
     713   
     714    // ** Here is the math to relate mag,zp to flux
    675715
    676716    // flux_cgs : erg sec^1 cm^-2 Hz^-1
    … …  
    703743    // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56)
    704744
     745    // get the zero point for the selected image
     746    // Use a different zero point for the PSF-like and APERTURE-like magnitudes
     747    float zpPSF  = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalPSF  + Mmos + Mgrid);
     748    float zpAPER = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalAPER + Mmos + Mgrid);
     749
    705750    // zpFactor to go from instrumental flux to Janskies
    706     float zpFactor = pow(10.0, -0.4*zp + 3.56);
     751    float zpFactorPSF  = pow(10.0, -0.4*zpPSF  + 3.56);
     752    float zpFactorAPER = pow(10.0, -0.4*zpAPER + 3.56);
    707753
    708754    // need to put in AB mag factor to get to Janskies (or uJy?)
    709     secfilt[Nsec].FpsfStk   = zpFactor * measure[k].FluxPSF; 
    710     secfilt[Nsec].dFpsfStk  = zpFactor * measure[k].dFluxPSF;
    711     secfilt[Nsec].FkronStk  = zpFactor * measure[k].FluxKron;
    712     secfilt[Nsec].dFkronStk = zpFactor * measure[k].dFluxKron;
    713     secfilt[Nsec].FapStk    = zpFactor * measure[k].FluxAp;
     755    secfilt[Nsec].FpsfStk   = zpFactorPSF  * measure[meas].FluxPSF; 
     756    secfilt[Nsec].dFpsfStk  = zpFactorPSF  * measure[meas].dFluxPSF;
     757    secfilt[Nsec].FkronStk  = zpFactorAPER * measure[meas].FluxKron;
     758    secfilt[Nsec].dFkronStk = zpFactorAPER * measure[meas].dFluxKron;
     759    secfilt[Nsec].FapStk    = zpFactorAPER * measure[meas].FluxAp;
    714760
    715761    // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245).  we are going to
    716762    // use PSF flux error instead here:
    717     // secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxAp;
    718     secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxPSF;
     763    // secfilt[Nsec].dFapStk   = zpFactorAPER * measure[meas].dFluxAp;
     764    secfilt[Nsec].dFapStk   = zpFactorAPER * measure[meas].dFluxPSF;
    719765
    720766    // Jy to AB mags
    721     secfilt[Nsec].MpsfStk   = (measure[k].FluxPSF  > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN;
    722     secfilt[Nsec].MkronStk  = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN;
    723     secfilt[Nsec].MapStk    = (measure[k].FluxAp   > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN;
    724 
    725     secfilt[Nsec].stackBestOff = k + measureOffset;
     767    secfilt[Nsec].MpsfStk   = (measure[meas].FluxPSF  > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk)  : NAN;
     768    secfilt[Nsec].MkronStk  = (measure[meas].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN;
     769    secfilt[Nsec].MapStk    = (measure[meas].FluxAp   > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk)   : NAN;
     770
     771    // record the measurement which gave the best value
     772    secfilt[Nsec].stackBestOff = (bestEntryMax == -1) ? -1 : bestEntryMax + measureOffset;
    726773    myAssert (secfilt[Nsec].stackBestOff <= catalog[0].Nmeasure, "stackBestOff out of range");
    727774
    728     secfilt[Nsec].Nstack    = Nstack;
     775    // record the selected primary measurement
     776    secfilt[Nsec].stackPrmryOff = (primaryEntryMax == -1) ? -1 : primaryEntryMax + measureOffset;
     777    myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range");
     778
     779    // this is the selected measurement used by secfilt[]
     780    measure[meas].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
     781    if (Nprimary) {
     782      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY;
     783      if (Nprimary > 1) {
     784        secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY_MULTIPLE;
     785      }
     786    }
     787
     788    // stack measurement used for 'best' was a detection (not forced from the other bands)
     789    if ((measure[bestEntryMax].photFlags2 & 0x00000004) == 0) {
     790      secfilt[Nsec].flags |= ID_SECF_STACK_BESTDET;
     791    }
     792
     793    // stack measurement used for 'primary' was a detection (not forced from the other bands)
     794    if ((primaryEntryMax >= 0) && ((measure[primaryEntryMax].photFlags2 & 0x00000004) == 0)) {
     795      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMDET;
     796    }
     797
     798    secfilt[Nsec].flags |= ID_SECF_HAS_PS1_STACK;
     799
     800    // NOTE: negative and insignificant flux values are allowed, but not NAN flux values
     801    // float Finst = PhotFluxInst (&measure[meas], MAG_CLASS_PSF);
     802   
     803    // data quality assessment
     804    isBad |= (measure[meas].photFlags & photcodes[Ns][0].photomBadMask);
     805    isBad |= (measure[meas].psfQF < 0.85);
     806    isBad |= isnan(measure[meas].psfQF);
     807    isBad |= measure[meas].dM > 0.2; // S/N < 5.0
     808
     809    isSuspect |= (measure[meas].photFlags & photcodes[Ns][0].photomPoorMask);
     810    isSuspect |= (measure[meas].psfQFperf < 0.85);
     811
     812    if (!isSuspect && !isBad) {
     813      NstackGood ++;
     814    }
     815    if (isSuspect && !isBad) {
     816      NstackSuspect ++;
     817    }
     818
     819    secfilt[Nsec].Nstack    = NstackMeas;
    729820    secfilt[Nsec].NstackDet = NstackDet;
    730821
    731     // this is the measurement used by secfilt[]
    732     measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
    733     if (k == stackPrimaryMeasure) {
    734       secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY;
    735     }
    736 
    737     // stack measurement used for 'best' was a detection (not forced from the other bands)
    738     if ((measure[k].photFlags2 & 0x00000004) == 0) {
    739       secfilt[Nsec].flags |= ID_SECF_STACK_BESTDET;
    740     }
    741     // stack measurement used for 'primary' was a detection (not forced from the other bands)
    742     if ((stackPrimaryMeasure >= 0) && ((measure[stackPrimaryMeasure].photFlags2 & 0x00000004) == 0)) {
    743       secfilt[Nsec].flags |= ID_SECF_STACK_PRIMDET;
    744     }
    745 
    746     secfilt[Nsec].flags |= ID_SECF_HAS_PS1_STACK;
    747822  } // Nsecfilt loop
    748823
    … …  
    755830
    756831  int PrimaryIsBest = TRUE;
    757   for (Ns = 0; Ns < Nphotcodes; Ns++) {
     832  for (int Ns = 0; Ns < Nphotcodes; Ns++) {
    758833    if (!(secfilt[Ns].flags & ID_SECF_HAS_PS1_STACK)) continue; // no stack detection in PS1, nothing is best
    759834    if (secfilt[Ns].flags & ID_SECF_STACK_PRIMARY) continue;    // primary stack detection is best
    … …  
    780855// * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr
    781856// analysis is done on flux, not mags (as the faintest objects will be nearly insignificant)
     857
     858// Although I calculate McalAPER for exposures, I am only using McalPSF for warps..
     859// Note in StarOps.c:setMcalOutput I am setting measure->McalAPER to image->McalPSF for
     860// chips and warps, but not stacks
     861
    782862int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
    783863  OHANA_UNUSED_PARAM(flatcorr);
    … …  
    857937      // data for which the associated image has not been loaded (probably because of
    858938      // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
    859       Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
     939      Mcal = measure[k].McalPSF; // check that this is zero for loaded REF value
    860940    } else {
    861       // use getMcal not getMcal_alt?
    862       Mcal  = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd);
    863       // Mcal  = getMcal (meas, cat);
     941      Mcal  = getMcal (meas, cat, MAG_CLASS_PSF);
    864942      if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
    865943    }
Note: See TracChangeset for help on using the changeset viewer.