[genvectorx] Add Invariant Masses test

Author: mdessole

math/experimental/genvectorx/test/CMakeLists.txt

@@ -8,14 +8,22 @@
 
 
     ROOT_EXECUTABLE(testGenvectorSYCL testGenVectorSYCL.cxx) # LIBRARIES GenVectorSYCL)
-    if (oneapi)
-        target_link_options(testGenvectorSYCL PUBLIC -E)
-    endif()
     target_link_libraries(testGenvectorSYCL PUBLIC GenVectorSYCL )
     add_sycl_to_root_target(TARGET testGenvectorSYCL 
     SOURCES testGenVectorSYCL.cxx 
     COMPILE_DEFINITIONS ROOT_MATH_SYCL
     DEPENDENCIES GenVectorSYCL
     )
+
+
+    ROOT_EXECUTABLE(testGenvectorSYCLInvMasses testGenVectorSYCLInvMasses.cxx) # LIBRARIES GenVectorSYCL)
+    target_link_libraries(testGenvectorSYCLInvMasses PUBLIC GenVectorSYCL )
+    add_sycl_to_root_target(TARGET testGenvectorSYCLInvMasses 
+    SOURCES testGenvectorSYCLInvMasses.cxx 
+    COMPILE_DEFINITIONS ROOT_MATH_SYCL
+    DEPENDENCIES GenVectorSYCL
+    )
+
     ROOT_ADD_TEST(test-genvector-genvectorsycl COMMAND testGenvectorSYCL)
+    ROOT_ADD_TEST(test-genvector-genvectorsyclinvmasses COMMAND testGenvectorSYCLInvMasses)
 

math/experimental/genvectorx/test/testGenVectorSYCLInvMasses.cxx

@@ -0,0 +1,173 @@
+#include "MathX/Vector3D.h"
+#include "MathX/Point3D.h"
+
+#include "MathX/Vector2D.h"
+#include "MathX/Point2D.h"
+
+#include "MathX/EulerAngles.h"
+
+#include "MathX/Transform3D.h"
+#include "MathX/Translation3D.h"
+
+#include "MathX/Rotation3D.h"
+#include "MathX/RotationX.h"
+#include "MathX/RotationY.h"
+#include "MathX/RotationZ.h"
+#include "MathX/Quaternion.h"
+#include "MathX/AxisAngle.h"
+#include "MathX/RotationZYX.h"
+
+#include "MathX/LorentzRotation.h"
+#include "MathX/PtEtaPhiM4D.h"
+#include "MathX/LorentzVector.h"
+
+#include "MathX/VectorUtil.h"
+
+#include "MathX/GenVectorX/AccHeaders.h"
+
+#include <sycl/sycl.hpp>
+
+#include <vector>
+
+using namespace ROOT::ROOT_MATH_ARCH;
+using namespace ROOT::ROOT_MATH_ARCH::VectorUtil;
+
+typedef LorentzVector<PtEtaPhiM4D<double>> vec4d;
+
+int compare(double v1, double v2, double scale = 1.0)
+{
+   //  ntest = ntest + 1;
+
+   // numerical double limit for epsilon
+   double eps = scale * std::numeric_limits<double>::epsilon();
+   int iret = 0;
+   double delta = v2 - v1;
+   double d = 0;
+   if (delta < 0)
+      delta = -delta;
+   if (v1 == 0 || v2 == 0) {
+      if (delta > eps) {
+         iret = 1;
+      }
+   }
+   // skip case v1 or v2 is infinity
+   else {
+      d = v1;
+
+      if (v1 < 0)
+         d = -d;
+      // add also case when delta is small by default
+      if (delta / d > eps && delta > eps)
+         iret = 1;
+   }
+
+   //    if (iret == 0)
+   //       std::cout << ".";
+   //    else {
+   //       int pr = std::cout.precision(18);
+   //       std::cout << "\nDiscrepancy in " << name << "() : " << v1 << " != " << v2 << " discr = " << int(delta / d /
+   //       eps)
+   //                 << "   (Allowed discrepancy is " << eps << ")\n";
+   //       std::cout.precision(pr);
+   //       // nfail = nfail + 1;
+   //    }
+   return iret;
+}
+
+template <class Transform>
+bool IsEqual(const Transform &t1, const Transform &t2, unsigned int size)
+{
+   // size should be an enum of the Transform class
+   std::vector<double> x1(size);
+   std::vector<double> x2(size);
+   t1.GetComponents(x1.begin(), x1.end());
+   t2.GetComponents(x2.begin(), x2.end());
+   bool ret = true;
+   unsigned int i = 0;
+   while (ret && i < size) {
+      // from TMath::AreEqualRel(x1,x2,2*eps)
+      bool areEqual =
+         std::abs(x1[i] - x2[i]) < std::numeric_limits<double>::epsilon() * (std::abs(x1[i]) + std::abs(x2[i]));
+      ret &= areEqual;
+      i++;
+   }
+   return ret;
+}
+
+vec4d *GenVectors(int n)
+{
+  vec4d *vectors = new vec4d[n];
+
+  // generate n -4 momentum quantities
+  for (int i = 0; i < n; ++i)
+  {
+    // fill vectors
+    vectors[i] = {1., 1., 1., 1.};
+  }
+
+  return vectors;
+}
+
+int testInvariantMasses(int N)
+{
+   int iret_host = 0;
+   std::cout << "testing Invariant Masses Computation \t:\n";
+
+   sycl::default_selector device_selector;
+   sycl::queue queue(device_selector);
+
+   auto v1 = GenVectors(N);
+   auto v2 = GenVectors(N);
+   double* invMasses = new double[N];
+
+   std::cout << "sycl::queue check - selected device:\n"
+                << queue.get_device().get_info<sycl::info::device::name>()
+                << std::endl;
+
+   {
+
+      vec4d* d_v1 = sycl::malloc_device<vec4d>(N, queue);
+      vec4d* d_v2 = sycl::malloc_device<vec4d>(N, queue);
+      double* d_invMasses = sycl::malloc_device<double>(N, queue);
+
+      queue.memcpy(d_v1, v1, N * sizeof(vec4d));
+      queue.memcpy(d_v2, v2, N * sizeof(vec4d));
+      queue.wait();
+
+      queue.submit([&](sycl::handler &cgh) {
+         cgh.parallel_for(sycl::range<1>(N), [=](sycl::id<1> indx) {
+            vec4d v = d_v2[indx]+d_v2[indx];
+            d_invMasses[indx] = v.M();
+
+         });
+      });
+
+      queue.wait();
+      queue.memcpy(invMasses, d_invMasses, N * sizeof(double));
+      queue.wait();
+   }
+
+   for (int i = 0; i < N; i++){
+      iret_host += (std::abs(invMasses[i] - 2.) > 1e-5);
+   }
+      
+   if (iret_host == 0)
+      std::cout << "\tOK\n";
+   else
+      std::cout << "\t FAILED\n";
+
+   return iret_host;
+
+}
+
+
+int main()
+{
+   int n = 128;
+   int ret = testInvariantMasses(n);
+   if (ret)
+      std::cerr << "test FAILED !!! " << std::endl;
+   else
+      std::cout << "test OK " << std::endl;
+   return ret;
+}