From b8d576ac83f411b82181da91b9978992f0890351 Mon Sep 17 00:00:00 2001 From: Juan Carlos Diaz-Velez Date: Wed, 3 Apr 2024 10:14:40 -0500 Subject: [PATCH] Added rotation test --- src/CMakeLists.txt | 20 +++++ src/include/iter-lhreco-proj/illh-utils.h | 2 +- src/private/iter-lhreco/illh-utils.cc | 2 +- src/private/iter-lhreco/multi-llh.cc | 6 +- src/private/test/rotation.cc | 96 +++++++++++++++++++++++ 5 files changed, 121 insertions(+), 5 deletions(-) create mode 100644 src/private/test/rotation.cc diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bd09051..e7ce2e6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -33,6 +33,12 @@ 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 @@ -40,6 +46,20 @@ target_link_libraries(multi-llh ${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 $ + ) + if(IS_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/scripts) install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/scripts diff --git a/src/include/iter-lhreco-proj/illh-utils.h b/src/include/iter-lhreco-proj/illh-utils.h index f82bb5c..902317c 100755 --- a/src/include/iter-lhreco-proj/illh-utils.h +++ b/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); } diff --git a/src/private/iter-lhreco/illh-utils.cc b/src/private/iter-lhreco/illh-utils.cc index 9bbef59..c3d8ece 100755 --- a/src/private/iter-lhreco/illh-utils.cc +++ b/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); diff --git a/src/private/iter-lhreco/multi-llh.cc b/src/private/iter-lhreco/multi-llh.cc index 0fbfc47..8d39e14 100755 --- a/src/private/iter-lhreco/multi-llh.cc +++ b/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]; diff --git a/src/private/test/rotation.cc b/src/private/test/rotation.cc new file mode 100644 index 0000000..3cc2739 --- /dev/null +++ b/src/private/test/rotation.cc @@ -0,0 +1,96 @@ +/** + * cra-llh-ahlers + * + * @version $Id: $ + * + * @date: $Date: $ + * + * @author Juan Carlos Diaz-Velez + * +*/ + + +#include +#include + +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 < 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; +} +