Changeset 38153 for trunk/Ohana/src/dvolens

Timestamp:

Apr 18, 2015, 6:33:51 AM (11 years ago)

Author:

eugene

Message:

merge changes from eam branch ipp-20150405: add R7 to lensing, fix R5,R6,R7 definitions and mean photometry; calculate E1,E2 in dvolens; fix precision for dvopsps FW objects; fix histograms at axis; new colormap options; require primary skycell in diff mean photometry

Location:

trunk

Files:

• . (modified) (1 prop)
• Ohana/src/dvolens/Makefile (modified) (2 diffs)
• Ohana/src/dvolens/include/dvolens.h (modified) (2 diffs)
• Ohana/src/dvolens/src/myIndex.c (copied) (copied from branches/eam_branches/ipp-20150405/Ohana/src/dvolens/src/myIndex.c )
• Ohana/src/dvolens/src/update_objects_catalog.c (modified) (11 diffs)

trunk/Ohana/src/dvolens/Makefile

@@ -29 +29 @@
 $(SRC)/help.$(ARCH).o            \
 $(SRC)/initialize.$(ARCH).o      \
+$(SRC)/myIndex.$(ARCH).o         \
 $(SRC)/update_objects.$(ARCH).o  \
 $(SRC)/update_objects_catalog.$(ARCH).o  \
@@ -44 +45 @@
 $(SRC)/help.$(ARCH).o            \
 $(SRC)/initialize.$(ARCH).o      \
+$(SRC)/myIndex.$(ARCH).o         \
 $(SRC)/update_objects.$(ARCH).o  \
 $(SRC)/update_objects_catalog.$(ARCH).o  \

trunk/Ohana/src/dvolens/include/dvolens.h

@@ -4 +4 @@
 # include <signal.h>
 # include <pthread.h>
+
+# ifndef MAX_INT
+# define MAX_INT 2147483647
+# endif
+
+typedef struct {
+  int minID;
+  int maxID;
+  int *index;
+  int Nindex;
+  int NINDEX;
+} myIndexType;
 
 typedef enum {
@@ -63 +75 @@
 int           client_logger_message   PROTO((char *format,...));
 
+myIndexType *myIndexAlloc ();
+int myIndexFree (myIndexType *myIndex);
+void myIndexInit (myIndexType *myIndex);
+int myIndexUpdateLimits (myIndexType *myIndex, int value);
+int myIndexSetRange (myIndexType *myIndex);
+int myIndexSetEntry (myIndexType *myIndex, int value, int entry);
+int myIndexGetEntry (myIndexType *myIndex, int value);

trunk/Ohana/src/dvolens/src/update_objects_catalog.c

@@ -1 +1 @@
 # include "dvolens.h"
+# define SCALE 0.001
+float MagToFlux (float Mag); // in libdvo, but not exposed?
 
 int update_objects_catalog (Catalog *catalog) {
@@ -11 +13 @@
   int Nsecfilt = GetPhotcodeNsecfilt ();
   
+  float *Mxx_obj = NULL;
+  float *Mxy_obj = NULL;
+  float *Myy_obj = NULL;
+
+  ALLOCATE (Mxx_obj, float, Nsecfilt);
+  ALLOCATE (Mxy_obj, float, Nsecfilt);
+  ALLOCATE (Myy_obj, float, Nsecfilt);
+
   // I think we have already allocated lensobj
   REALLOCATE (catalog->lensobj, Lensobj, catalog->Naverage * Nsecfilt);
@@ -19 +29 @@
   int Nlensobj = 0;
 
+  myIndexType *measureIndex = myIndexAlloc();
+
   for (i = 0; i < catalog->Naverage; i++) {
     Average *average = &catalog->average[i];
@@ -27 +39 @@
     for (j = 0; j < Nsecfilt; j++) {
       dvo_lensobj_init (&lensobj[j], TRUE); // init accumulated values to 0
+      Mxx_obj[j] = 0.0;
+      Mxy_obj[j] = 0.0;
+      Myy_obj[j] = 0.0;
     }
 
     // reset the average values
-
     off_t Loff = average->lensingOffset;
     off_t Moff = average->measureOffset;
+
+    // I have Nmeasure entries for this object.  I need to match measure to lensing (by imageID)
+    // I could generated a sorted list (imageID, measureSeq) and use bisection to find the desired imageID
+    // I could make an index measureSeq[imageID-imageIDmin]
+    
+    myIndexInit (measureIndex);
+
+    // generate an index for these measure entries (based on Mj and imageID)
+    for (Mj = 0; Mj < average->Nmeasure; Mj++) {
+      myIndexUpdateLimits (measureIndex, catalog->measure[Moff + Mj].imageID);
+    }
+    myIndexSetRange (measureIndex);
+    for (Mj = 0; Mj < average->Nmeasure; Mj++) {
+      // if there are duplicates (multiple measures from the same image), 
+      // myIndexSetEntry will only select one
+      myIndexSetEntry (measureIndex, catalog->measure[Moff + Mj].imageID, Mj);
+    }
 
     // loop over the lensing measurements.  for each one, I need to find the corresponding measurement (make an index in lensing?)
@@ -41 +72 @@
       if (!isfinite(lensing->X11_sm_obj)) continue;
 
+      // pointer from lensing entry to corresponding measure entry (keep updated?)
+      // XX ALT int Mseq = lensing->measureSeq;
+      // XX ALT myAssert (Mseq < average->Nmeasure, "oops");
+      // XX ALT Measure *measure = &catalog->measure[Moff + Mseq];
+      // XX ALT myAssert (measure->imageID == lensing->imageID, "oops, deux");
+
+      Mj = myIndexGetEntry (measureIndex, lensing->imageID);
+      myAssert (Mj > -1, "missing index");
+
+      Measure *measure = &catalog->measure[Moff + Mj];
+      myAssert (measure->imageID == lensing->imageID, "oops, deux");
+
+# if (0)      
       Measure *measure = NULL;
       int found = FALSE;
@@ -54 +98 @@
         abort();
       }
+# endif
       
       // skip measurements that do not match the current photcode
@@ -67 +112 @@
       if (measure->psfQFperf < 0.85) continue; 
       
-
+      // I should probably only use lensing entries which correspond to warps 
+      // included in the mean warp flux (setMrelCatalog.c)
+      if ((measure->dbFlags & ID_MEAS_WARP_USED) == 0) continue;
+
+      Mxx_obj[Nsec] += measure->Mxx;
+      Mxy_obj[Nsec] += measure->Mxy;
+      Myy_obj[Nsec] += measure->Myy;
+      
       lensobj[Nsec].X11_sm_obj += lensing->X11_sm_obj;
       lensobj[Nsec].X12_sm_obj += lensing->X12_sm_obj;
@@ -92 +144 @@
       lensobj[Nsec]. E2_sh_psf += lensing-> E2_sh_psf;
 
-      lensobj[Nsec]. F_ApR5 += lensing-> F_ApR5;
-      lensobj[Nsec].dF_ApR5 += SQ(lensing->dF_ApR5);
-      lensobj[Nsec].sF_ApR5 += SQ(lensing->sF_ApR5);
-      lensobj[Nsec].fF_ApR5 += lensing->fF_ApR5;
-
-      lensobj[Nsec]. F_ApR6 += lensing-> F_ApR6;
-      lensobj[Nsec].dF_ApR6 += SQ(lensing->dF_ApR6);
-      lensobj[Nsec].sF_ApR6 += SQ(lensing->sF_ApR6);
-      lensobj[Nsec].fF_ApR6 += lensing->fF_ApR6;
+      // relphot sets measure->Mcal (setMcalOutput.c, called by setMrelFinal.c)
+      float Mcal = code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure->Mcal;
+      float Fcal = 3630.8 * MagToFlux(Mcal);
+
+      // lensing->F_ApR5, etc need to be in units of DN/sec
+      // Fcal * lensing->F_ApR5 is in Jy
+
+      lensobj[Nsec]. F_ApR5 +=    Fcal * lensing-> F_ApR5;
+      lensobj[Nsec].dF_ApR5 += SQ(Fcal * lensing->dF_ApR5);
+      lensobj[Nsec].sF_ApR5 += SQ(Fcal * lensing->sF_ApR5);
+      lensobj[Nsec].fF_ApR5 +=           lensing->fF_ApR5;
+
+      lensobj[Nsec]. F_ApR6 +=    Fcal * lensing-> F_ApR6;
+      lensobj[Nsec].dF_ApR6 += SQ(Fcal * lensing->dF_ApR6);
+      lensobj[Nsec].sF_ApR6 += SQ(Fcal * lensing->sF_ApR6);
+      lensobj[Nsec].fF_ApR6 +=           lensing->fF_ApR6;
+
+      lensobj[Nsec]. F_ApR7 +=    Fcal * lensing-> F_ApR7;
+      lensobj[Nsec].dF_ApR7 += SQ(Fcal * lensing->dF_ApR7);
+      lensobj[Nsec].sF_ApR7 += SQ(Fcal * lensing->sF_ApR7);
+      lensobj[Nsec].fF_ApR7 +=           lensing->fF_ApR7;
 
       lensobj[Nsec].Nmeas   ++;
@@ -139 +203 @@
       lensobj[j].sF_ApR5 /= Nmeas;
       lensobj[j].fF_ApR5 /= Nmeas;
+
       lensobj[j]. F_ApR6 /= Nmeas;
       lensobj[j].dF_ApR6 /= Nmeas;
@@ -144 +209 @@
       lensobj[j].fF_ApR6 /= Nmeas;
 
+      lensobj[j]. F_ApR7 /= Nmeas;
+      lensobj[j].dF_ApR7 /= Nmeas;
+      lensobj[j].sF_ApR7 /= Nmeas;
+      lensobj[j].fF_ApR7 /= Nmeas;
+
       lensobj[j].dF_ApR5 = sqrt(lensobj[j].dF_ApR5);
       lensobj[j].sF_ApR5 = sqrt(lensobj[j].sF_ApR5);
       lensobj[j].dF_ApR6 = sqrt(lensobj[j].dF_ApR6);
       lensobj[j].sF_ApR6 = sqrt(lensobj[j].sF_ApR6);
+      lensobj[j].dF_ApR7 = sqrt(lensobj[j].dF_ApR7);
+      lensobj[j].sF_ApR7 = sqrt(lensobj[j].sF_ApR7);
 
       // somewhere in here I should calculate the values of E1, E2, and gamma
+
+      float e0 = Mxx_obj[j] + Myy_obj[j];
+      lensobj[j].E1 = (Mxx_obj[j] - Myy_obj[j]) / e0;
+      lensobj[j].E2 = 2.0 * Mxy_obj[j] / e0;
     }
 
@@ -159 +235 @@
   }
 
+  free (Mxx_obj);
+  free (Mxy_obj);
+  free (Myy_obj);
+
+  myIndexFree (measureIndex);
+
   catalog->Nlensobj = Nlensobj;
   catalog->Nlensobj_disk = 0;