Changeset 35282

Timestamp:

Mar 9, 2013, 6:33:50 AM (13 years ago)

Author:

eugene

Message:

working on high-speed pm

File:

• branches/eam_branches/ipp-20130306/Ohana/src/relastro/src/hpm_objects.c (modified) (9 diffs)

branches/eam_branches/ipp-20130306/Ohana/src/relastro/src/hpm_objects.c

@@ -1 +1 @@
 # include "relastro.h"
 
-# define NEXT_I { i++; continue; }
+# define NEXT_I { if (Ngroup < 2) slowMoving[ni] = TRUE; newI = TRUE; i++; continue; }
 # define NEXT_J { j++; continue; }
 
-int high_speed_objects (SkyRegion *region, Catalog *catalog) {
+int hpm_objects (SkyRegion *region, Catalog *catalog) {
 
   off_t i, j, m, J, ni, nj, *N1;
@@ -15 +15 @@
   Coords tcoords;
   Catalog catalogOut;
+  Catalog testcat;
 
   int Nsecfilt;
   char filename[1024];
 
+  // XXX are we saving these in an hpm dvodb?
   snprintf (filename, 1024, "%s/%s.cpt", HIGH_SPEED_DIR, region[0].name);
   fprintf (stderr, "%s\n",filename);
@@ -49 +51 @@
   REALLOCATE (catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt);
 
-  // high-speed between different surveys (easier case):
-  // we have two sets of photcodes (A) and (B) and are looking for objects
-  // with detections in only (A) and separately only (B).
+  // testcat is used to determine the fit for a single object group
+  // objects which do not have a high-quality testcat fit are not kept
+  dvo_catalog_init (&testcat, TRUE);
+  testcat.Naverage = 1; // this is fixed -- only one obj in testcat
+  NMEASURE_TEST = 1000;
+  REALLOCATE (catalogOut.average, Average, 1);
+  REALLOCATE (catalogOut.measure, Measure, NMEASURE_TEST);
+  REALLOCATE (catalogOut.secfilt, SecFilt, Nsecfilt);
 
   // we need at least 2 objects if we are going to match anything...
@@ -59 +66 @@
   ALLOCATE (slowMoving, int, catalog[0].Naverage);
   memset (slowMoving, 0, catalog[0].Naverage*sizeof(int));
-
-  // record to which photcode group the object belongs:
-  NgroupA = 0;
-  ALLOCATE (groupA, int, catalog[0].Naverage);
-  memset (groupA, 0, catalog[0].Naverage*sizeof(int));
-
-  NgroupB = 0;
-  ALLOCATE (groupB, int, catalog[0].Naverage);
-  memset (groupB, 0, catalog[0].Naverage*sizeof(int));
 
   if (VERBOSE) fprintf (stderr, "checking "OFF_T_FMT" objects\n",  catalog[0].Naverage);
@@ -84 +82 @@
     }
 
-    // do any of the measures for this object match group A?
-    m = catalog[0].average[i].measureOffset;
-    foundA = FALSE;
-    for (j = 0; !foundA && (j < catalog[0].average[i].Nmeasure); j++, m++) {
-
-      if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupA, NphotcodesGroupA)) {
-        foundA = TRUE;
-      }
-    }
-
-    // do any of the measures for this object match group B?
-    m = catalog[0].average[i].measureOffset;
-    foundB = FALSE;
-    for (j = 0; !foundB && (j < catalog[0].average[i].Nmeasure); j++, m++) {
-          
-      if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupB, NphotcodesGroupB)) {
-        foundB = TRUE;
-      }
-    }
-
-    // object found in both - mark as slow moving
-    if (foundA && foundB) {
+    // selection criteria:
+    // (Nps1 > 2) && (Trange < 180)
+    if (catalog[0].average[i].Trange > MAX_TRANGE) {
       slowMoving[i] = TRUE;
       Nslow ++;
       continue;
     }
-
-    // Apply additional constraints:
-    if (foundA && !foundB) {
-        if (applyConstraintsA(catalog, i)) {
-            groupA[i] = TRUE;
-            NgroupA++;
-            continue;
-        } else {
-            NgroupAbad ++;
-            // fprintf (stderr, "skip "OFF_T_FMT" (group A)\n",  i);
-        }
-        continue;
-    } 
-    if (foundB && !foundA) {
-        if (applyConstraintsB(catalog, i)) {
-            groupB[i] = TRUE;
-            NgroupB++;
-            continue;
-        } else {
-            NgroupBbad++;
-        }
-        continue;
-    }
-    // this object does not have a detection matching the contraints from either gorupA or groupB -- skip as if slow
-    slowMoving[i] = TRUE;
-    Ninvalid ++;
-  }
-
-  fprintf (stderr, OFF_T_FMT" slow, "OFF_T_FMT" invalid, ("OFF_T_FMT" group A, "OFF_T_FMT" group B), "OFF_T_FMT" total objects; "OFF_T_FMT" in group A, "OFF_T_FMT" in group B\n",  Nslow,  Ninvalid, NgroupAbad, NgroupBbad, catalog[0].Naverage,  NgroupA,  NgroupB);
+    // count the PS1 detections via explicit photcode ranges?
+    // XXX this is a total hard-wired hack...
+    int Nps1 = 0;
+    for (j = 0; i < 5; j++) {
+      Nps1 += catalog[0].secfilt[Nsecfilt*i+j].Ncode;
+    }
+    if (Nps1 < MIN_PS1_DET) {
+      slowMoving[i] = TRUE;
+      Nslow ++;
+      continue;
+    }      
+  }
+
+  fprintf (stderr, OFF_T_FMT" slow, "OFF_T_FMT" total objects; "OFF_T_FMT" possible fast\n",  Nslow,  catalog[0].Naverage,  catalog[0].Naverage - Nslow);
     
   // double loop over unmarked objects (sorted in RA / X)
@@ -175 +139 @@
   RADIUS2 = SQ(RADIUS);
 
-  // mark (exclude) objects with both sets of target photcodes
+  // group is a list of objects that are within a clump.  this set will be tested
+  // via a clipped fit to the measurements after the group is identified
+  Ngroup = 0;
+  NGROUP = 100;
+  ALLOCATE (group, int, NGROUP);
+
+  // in the loop below, we need to do a bunch of things when we go to the next main object
+  newI = TRUE;
+
+  // mark (skip) objects with both sets of target photcodes
+  // the loop below is attempting to find associations of multiple objects which have 
+  // passed the cuts above.  the index i is following the primary object of interest
+  // the index j is used to explore possible near neighbors.  
+  // When we go to the next object 'i', Nmatch is reset 
   for (i = j = 0; (i < catalog[0].Naverage) && (j < catalog[0].Naverage);) {
 
     ni = N1[i];
     nj = N1[j];
+
+    if (newI) {
+      Nmatch = 0; 
+      Ngroup = 1;
+      group[0] = ni;
+      newI = FALSE;
+    }
 
     XVERB  = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
@@ -190 +174 @@
     if (slowMoving[nj]) NEXT_J;
 
-    // i => groupA, j => groupB
-    if (!groupA[ni]) NEXT_I;
-    if (!groupB[nj]) NEXT_J;
-
     if (!finite(X1[i]) || !finite(Y1[i])) NEXT_I;
     if (!finite(X1[j]) || !finite(Y1[j])) NEXT_J;
@@ -209 +189 @@
     for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) {     
       if (J == i) continue;  // avoid auto-matches
+      nj = N1[J];
 
       dX = X1[i] - X1[J];
 
-      nj = N1[J];
-      if (!groupB[nj]) continue;
+      if (slowMoving[nj]) continue;
 
       XVERB  = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
@@ -226 +206 @@
       if (dR > RADIUS2) continue;
 
-      /*** a match is found (just print it for the moment) ***/
-      /*Define a vector of matched indices and set averages in new catalogue*/
-      Nepoch=2;
-      FIT_MODE = FIT_PM_ONLY;
-      nv[0]=ni; /*THESE SHOULD BE CHANGED FOR MULTIPLE EPOCHS AS SHOULD nv*/
-      nv[1]=nj;
-
-      catalogOut.average[Nmatch]=catalog[0].average[nv[0]];
-      /*Loop over index values and set measurements in new catalogue*/
-      Nmatchmeasobj=0;
-      catalogOut.average[Nmatch].measureOffset=Nmatchmeas;
-      for(l=0;l<Nepoch;l++) {
-          m = catalog[0].average[nv[l]].measureOffset;
-          for(i1=0;i1<catalog[0].average[nv[l]].Nmeasure;i1++) {
-              catalogOut.measure[Nmatchmeas]=catalog[0].measure[m+i1];
-              /*Set offset RA and Dec wrt correct average value*/
-              catalogOut.measure[Nmatchmeas].dR=catalog[0].measure[m+i1].dR+3600.0*(catalog[0].average[nv[0]].R-catalog[0].average[nv[l]].R);
-              catalogOut.measure[Nmatchmeas].dD=catalog[0].measure[m+i1].dD+3600.0*(catalog[0].average[nv[0]].D-catalog[0].average[nv[l]].D);
-              catalogOut.measure[Nmatchmeas].averef = Nmatch;
-              Nmatchmeasobj++;
-              Nmatchmeas++;
-
-              if (Nmatchmeas == NMEASURE - 1) {
-                NMEASURE += 10000;
-                REALLOCATE (catalogOut.measure, Measure, NMEASURE);
-              }
-          }
+      /*** a match is found ***/
+      group[Ngroup] = nj;
+      Ngroup ++;
+      CHECK_REALLOCATE (group, int, NGROUP, Ngroup, 100);
+    }
+
+    if (Ngroup < 2) NEXT_I;
+
+    // we now have spatially associated group of objects.  now we need to see if the set of
+    // measurements can be fitted reasonably with a proper motion (& parallax?)
+
+    // the mean object will start with info from the primary object
+    // remember: testcat.Naverage = 1 -- does not change
+    testcat.average[0] = catalog[0].average[group[0]];
+    testcat.average[0].measureOffset = 0;
+    Nmatchmeas = 0;
+    I = group[0];
+    for (J = 0; J < Ngroup; J++) {
+      J = group[J];
+      m = catalog[0].average[Nj].measureOffset;
+      for (k = 0; k < catalog[0].average[Nj].Nmeasure; k++) {
+        testcat.measure[Nmatchmeas] = catalog[0].measure[m+k];
+        /* Set offset RA and Dec wrt correct average value*/
+        testcat.measure[Nmatchmeas].dR = catalog[0].measure[m+k].dR + 3600.0*(catalog[0].average[I].R - catalog[0].average[J].R);
+        testcat.measure[Nmatchmeas].dD = catalog[0].measure[m+k].dD + 3600.0*(catalog[0].average[I].D - catalog[0].average[J].D);
+        testcat.measure[Nmatchmeas].averef = 0;
+        Nmatchmeas++;
+        CHECK_REALLOCATE (testcat.measure, Measure, NMEASURE_TEST, Nmatchmeas, 1000);
       }
-      catalogOut.average[Nmatch].Nmeasure = Nmatchmeasobj;
-      Nmatch ++;
-
-      if (Nmatch == NAVERAGE - 1) {
-        NAVERAGE += 1000;
-        REALLOCATE (catalogOut.average, Average, NAVERAGE);
-        REALLOCATE(catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt);
+    }
+    testcat.average[0].Nmeasure = Nmatchmeas;
+
+    // we have now accumulated the measurements for this group, let's try a fit
+    // this needs to be a (fairly robust) clipped fit or we will have a hard time
+    // distinguishing a bad fit from a fit with 1 or 2 bad points
+    FIT_MODE = FIT_PM_ONLY;
+    UpdateObjects (&testcat, 1);
+
+    // logic for keeping the fit:
+    if (testcat.average[0].ChiSqPM < XXX) good = TRUE;
+    if (testcat.average[0].Npos > 0.X * testcat.average[0].Nmeasure) good = TRUE;
+    // what else?
+
+    if (good) {
+      // save the new object on catalogOut
+      {
+        catalogOut.average[Nout] = testcat.average[0];
+        catalogOut.average[Nout].measureOffset = Nmatchmeas;
+        m = testcat.average[0].measureOffset;
+        for (k = 0; k < testcat.average[0].Nmeasure; k++) {
+          catalogOut.measure[Noutmeas] = testcat.measure[m+k];
+          catalogOut.measure[Noutmeas].averef = Nout;
+        }
+        Nout ++;
+        CHECK_REALLOCATE (catalogOut.average, Average, NAVERAGE, Nout, 100);
       }
     }
-    i++;
-  }
+    NEXT_I;
+  }
+
   catalogOut.Naverage=Nmatch;
   catalogOut.Nmeasure=Nmatchmeas;