Added rotation test

src/CMakeLists.txt

@@ -33,13 +33,33 @@ add_executable(multi-llh
             private/iter-lhreco/illh-utils.cc
             )
 
+add_executable(rotation-test
+            private/test/rotation.cc
+            private/iter-lhreco/illh-utils.cc
+            )
+
+
 
 
 target_link_libraries(multi-llh 
                 ${CFITSIO_LIBRARIES} 
                 ${HEALPIX-CXX_LIBRARIES} 
                 ${Boost_LIBRARIES})
 
+target_link_libraries(rotation-test PRIVATE
+                ${CFITSIO_LIBRARIES} 
+                ${HEALPIX-CXX_LIBRARIES} 
+                ${Boost_LIBRARIES})
+
+
+enable_testing()
+
+# define tests
+add_test(
+  NAME rotation_test
+  COMMAND $<TARGET_FILE:rotation-test>
+  )
+
 
 if(IS_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/scripts)
        install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/scripts

src/include/iter-lhreco-proj/illh-utils.h

@@ -154,7 +154,7 @@ namespace illh
      * @return rotated pixel id
      */
     unsigned 
-    eq2log_idx(unsigned i, unsigned timeidx, double lat, double lon, unsigned nTimesteps, SkyMap& CRmap);
+    eq2loc_idx(unsigned i, unsigned timeidx, double lat, double lon, unsigned nTimesteps, SkyMap& CRmap);
 
 
 }

src/private/iter-lhreco/illh-utils.cc

@@ -285,7 +285,7 @@ void illh::save_iter(
  * @return rotated pixel id
  */
 unsigned 
-illh::eq2log_idx(unsigned i, unsigned timeidx, double lat, double lon, unsigned nTimesteps, SkyMap& CRmap)
+illh::eq2loc_idx(unsigned i, unsigned timeidx, double lat, double lon, unsigned nTimesteps, SkyMap& CRmap)
 { 
         vec3 v =  CRmap.pix2vec(i); 
         double clat = cos(lat); 

src/private/iter-lhreco/multi-llh.cc

@@ -442,7 +442,7 @@ int main(int argc, char* argv[])
                 {
                     DetectorPtr det = *det_it;
                     int j; 
-                    j = illh::eq2log_idx(i, timeidx, det->latitude, det->longitude, nTimesteps, CRmap);
+                    j = illh::eq2loc_idx(i, timeidx, det->latitude, det->longitude, nTimesteps, CRmap);
                     if ((*det->Emap0)[j] > 0.0 ){ 
                         nEvents += (*det->Nmap[timeidx])[j];
                         nBkg    += det->norm[timeidx]*(*det->Emap0)[j]; 
@@ -499,7 +499,7 @@ int main(int argc, char* argv[])
                     for (det_it = detectors.begin(); det_it != detectors.end(); det_it++)
                     {
                             DetectorPtr det = *det_it;
-                            int j = illh::eq2log_idx(i, timeidx, det->latitude, det->longitude, nTimesteps, CRmap);
+                            int j = illh::eq2loc_idx(i, timeidx, det->latitude, det->longitude, nTimesteps, CRmap);
                             if (det->FOV[j] && ((*det->Emap)[j] > 0.0)) { 
                                 nEvents += (*det->Nmap[timeidx])[j]; 
                                 nBkg += det->norm[timeidx]*(*det->Emap)[j]; 
@@ -640,7 +640,7 @@ int main(int argc, char* argv[])
                         { 
                             DetectorPtr det = *det_it;
                             //rotation from Equatorial (ra,dec) to local (theta,phi)
-                            int j = illh::eq2log_idx(i, timeidx, det->latitude, det->longitude, nTimesteps, CRmap);
+                            int j = illh::eq2loc_idx(i, timeidx, det->latitude, det->longitude, nTimesteps, CRmap);
 
                             // global significance 
                             double ej = (*det->Emap)[j];

src/private/test/rotation.cc

@@ -0,0 +1,96 @@
+/**
+ * cra-llh-ahlers
+ *
+ * @version $Id: $
+ *
+ * @date: $Date: $
+ *
+ * @author Juan Carlos Diaz-Velez <[email protected]>
+ *
+*/
+
+
+#include <iter-lhreco-proj/illh-utils.h>
+#include <cmath>
+
+double angular_distance(pointing p0, pointing p1)
+{
+   return acos(cos(p0.theta)*cos(p1.theta) + sin(p0.theta)*sin(p1.theta)*cos(p0.phi-p1.phi));
+}
+
+int main() {
+
+
+  double latitude = 18.0*M_PI/180.;
+  double longitude = -97.0*M_PI/180.;
+
+  unsigned nside(64);
+  unsigned npix = 12*nside*nside;
+  float pixel_area = 4*M_PI/npix;
+  float pixel_size = sqrt(pixel_area);
+
+  SkyMap CRmap;
+  CRmap.SetNside(nside, RING);
+  CRmap.fill(1.0);
+
+  double crab_ra = (5.+  (34. + 32./60.)/60)/24.*2*M_PI;
+  double crab_dec = (22. +52/3600)*M_PI/180.;
+  double aries_ra = 0;
+  double aries_dec = 0;
+
+  pointing aries(M_PI - aries_ra, aries_ra);
+  pointing crab(M_PI - crab_ra, crab_ra);
+
+  int timeidx = 0;
+  int nTimesteps = 360;
+
+  //rotation from Equatorial (ra,dec) to local (theta,phi)
+  auto integers = {0, 20, 101, 235, 300, 400, 500, 1000};
+  int j,k;
+
+  for (auto i : integers) 
+  {
+    j = illh::eq2loc_idx(i, timeidx, latitude, longitude, nTimesteps, CRmap);
+    k = illh::loc2eq_idx(j, timeidx, latitude, longitude, nTimesteps, CRmap);
+
+    // get theta phi from pixels
+    pointing pt0 = CRmap.pix2ang(i); 
+    pointing pt1 = CRmap.pix2ang(k); 
+    auto ang_dist = angular_distance(pt0,pt1);
+
+    if (ang_dist > 2*pixel_size)  // tolerance is one pixel
+    {
+       std::cout << i << ", " << k << "diff = "<< ang_dist << ": "<< pixel_size <<std::endl;
+       return 1;
+    }
+  }
+
+
+  for (auto t = 0; t < nTimesteps; t++)
+  {
+    auto ak = CRmap.ang2pix(aries); 
+    auto aj = illh::eq2loc_idx(ak, t, latitude, longitude, nTimesteps, CRmap);
+    auto ai = illh::loc2eq_idx(aj, t, latitude, longitude, nTimesteps, CRmap);
+
+    auto ck = CRmap.ang2pix(crab); 
+    auto cj = illh::eq2loc_idx(ck, t, latitude, longitude, nTimesteps, CRmap);
+    auto ci = illh::loc2eq_idx(cj, t, latitude, longitude, nTimesteps, CRmap);
+
+    // get theta phi from pixels
+    pointing pt0 = CRmap.pix2ang(ai); 
+    pointing pt1 = CRmap.pix2ang(ci); 
+
+    if (angular_distance(pt0,aries) > 2*pixel_size)  // tolerance is one pixel
+    {
+       std::cout << boost::format("Aries (%f,%f) -> (%f,%f)") % aries.theta % aries.phi % pt0.theta % pt0.phi << std::endl;
+       return 1;
+    }
+    if (angular_distance(pt0,aries) > 2*pixel_size)  // tolerance is one pixel
+    {
+       std::cout << boost::format("Crab (%f,%f) -> (%f,%f)") % crab.theta % crab.phi % pt1.theta % pt1.phi << std::endl;
+       return 1;
+    }
+  }
+  return 0;
+}
+