Changeset 31272 for branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/StarOps.c

Timestamp:

Apr 12, 2011, 6:06:39 AM (15 years ago)

Author:

eugene

Message:

fix the multifilter analysis

File:

• branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/StarOps.c (modified) (16 diffs)

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/StarOps.c

@@ -29 +29 @@
 float getMrel (Catalog *catalog, off_t meas, int cat) {
 
-  int Nsec, Nsecfilt, code;
+  int Nsec, Nsecfilt, photcode;
   int ave;
   float value;
 
   ave = catalog[cat].measure[meas].averef;
-  code = catalog[cat].measure[meas].photcode;
-
-  // XXX DEP : this should use a flag associated with secfilt
-  if (catalog[cat].average[ave].flags & STAR_BAD) return (NAN);  
+  photcode = catalog[cat].measure[meas].photcode;
+
+  int ecode = GetPhotcodeEquivCodebyCode (photcode);
+  Nsec = GetPhotcodeNsec(ecode);
+  Nsecfilt = GetPhotcodeNsecfilt ();
+
+  // is this star OK?
+  if (catalog[cat].secfilt[Nsecfilt*ave+Nsec].flags & STAR_BAD) return (NAN);  
  
-  Nsec = GetPhotcodeNsec(code);
-  Nsecfilt = GetPhotcodeNsecfilt ();
-
   value = catalog[cat].secfilt[Nsecfilt*ave+Nsec].M;
   return (value);
@@ -52 +53 @@
   float Msys, Mcal, Mmos, Mgrid;
   StatType stats;
-  int Nsec, Nsecfilt, ecode;
-
-  Nsecfilt = GetPhotcodeNsecfilt ();
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
   Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0;
 
@@ -64 +64 @@
 
         int thisCode = photcodes[Ns][0].code;
-        Nsec = GetPhotcodeNsec(thisCode);
+        int Nsec = GetPhotcodeNsec(thisCode);
 
         /* calculate the average value for a single star */
-        // XXX this flag should be set by secfilt, not average
-        if (catalog[i].average[j].flags & STAR_BAD) continue;  
+
+        // skip bad stars
+        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
         m = catalog[i].average[j].measureOffset;
 
@@ -74 +75 @@
         for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
           // skip measurements that do not match the current photcode
-          ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
+          int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
           if (ecode != thisCode) { continue; }
 
@@ -123 +124 @@
 
         // when performing the grid analysis, STAR_TOOFEW will be set to 1;
-        // XXX DEP this makes no sense: need a separate flag for each secfilt
         if (N <= STAR_TOOFEW) { /* too few measurements */
-          catalog[i].average[j].flags |= ID_STAR_FEW;
+          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
           Nfew ++;
         } else {
-          catalog[i].average[j].flags &= ~ID_STAR_FEW;
+          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
         }       
 
@@ -408 +408 @@
     for (i = Ntot = 0; i < Ncatalog; i++) {
       for (j = 0; j < catalog[i].Naverage; j++) {
-        if (catalog[i].average[j].flags & STAR_BAD) continue;
+        if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue;
         Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
         if (Xm == -1) continue;
@@ -434 +434 @@
         mark = (dM > MaxScatter) || (Xm == NAN_S_SHORT) || (Chisq > MaxChisq);
         if (mark) {
-          catalog[i].average[j].flags |= ID_STAR_POOR;
+          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR;
           Ndel ++;
         } else {
-          catalog[i].average[j].flags &= ~ID_STAR_POOR;
+          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR;
         }
         Nave ++;
@@ -492 +492 @@
 
         int thisCode = photcodes[Ns][0].code;
-
+        int Nsec = GetPhotcodeNsec(thisCode);
+        
         /* skip bad stars to prevent them from becoming good (on inner sample) */
-        if (catalog[i].average[j].flags & STAR_BAD) continue;  
+        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
 
         /* accumulate list of valid measurements */
@@ -582 +583 @@
 }
 
-StatType statsStarN (Catalog *catalog, int Ncatalog, int seccode) {
+StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode) {
 
   off_t j, k, m, Ntot;
@@ -589 +590 @@
   float Mcal, Mmos, Mgrid;
   StatType stats;
+  int N1, N2, N3, N4, N0;
+  N1 = N2 = N3 = N4 = N0 = 0;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
 
   Ntot = 0;
@@ -603 +608 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) { N1++; continue;  }
       m = catalog[i].average[j].measureOffset;
-
-      int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
-      if (ecode == seccode) continue;
 
       N = 0;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+        int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
+        if (ecode != seccode) { N0++; continue;}
         Mcal = getMcal  (m, i);
-        if (isnan(Mcal)) continue;
+        if (isnan(Mcal)) { N2++; continue;}
         Mmos = getMmos  (m, i);
-        if (isnan(Mmos)) continue;
+        if (isnan(Mmos)) { N3++; continue; }
         Mgrid = getMgrid (m, i);
-        if (isnan(Mgrid)) continue;
+        if (isnan(Mgrid)) { N4++; continue;}
         N++;
       }
@@ -626 +630 @@
   }
 
+  fprintf (stderr, "N1: %d, N2: %d, N3: %d, N4: %d, N0: %d\n", N1, N2, N3, N4, N0);
   liststats (list, dlist, n, &stats);
   free (list);
@@ -655 +660 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
         
       Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
@@ -694 +699 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
 
       dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
@@ -733 +738 @@
     for (i = 0; i < Ncatalog; i++) {
       for (j = 0; j < catalog[i].Naverage; j++) {
-        if (catalog[i].average[j].flags & STAR_BAD) continue;  
+        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
         dMrel = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
         bin = dMrel / 0.00025;
@@ -773 +778 @@
     for (i = 0; i < Ncatalog; i++) {
       for (j = 0; j < catalog[i].Naverage; j++) {
-        if (catalog[i].average[j].flags & STAR_BAD) continue;
+        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
         xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].M;
         value    = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;