Few TRestTools/TRestPhysics/TRestStringHelper upgrades motivated by TRestAxionOpticsMirror

cmake/Testing.cmake

@@ -19,10 +19,8 @@ endmacro()
 macro(ADD_LIBRARY_TEST)
     if (TEST)
         message(STATUS "Adding tests at ${CMAKE_CURRENT_SOURCE_DIR}")
-        set(TESTING_EXECUTABLE testRestGeant4)
 
         get_filename_component(DIR_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
-
         string(SUBSTRING ${DIR_NAME} 0 1 FIRST_LETTER)
         string(TOUPPER ${FIRST_LETTER} FIRST_LETTER)
         string(REGEX REPLACE "^.(.*)" "${FIRST_LETTER}\\1" DIR_NAME_CAPITALIZED "${DIR_NAME}")

source/framework/tools/inc/TRestPhysics.h

@@ -25,8 +25,10 @@
 
 #include <iostream>
 
+#include "TMatrixD.h"
 #include "TString.h"
 #include "TVector3.h"
+#include "TVectorD.h"
 
 /// This namespace serves to define physics constants and other basic physical operations
 namespace REST_Physics {
@@ -68,7 +70,15 @@ TVector3 MoveToPlane(TVector3 pos, TVector3 dir, TVector3 n, TVector3 a);
 TVector3 MoveByDistance(TVector3 pos, TVector3 dir, Double_t d);
 TVector3 MoveByDistanceFast(TVector3 pos, TVector3 dir, Double_t d);
 
-TVector3 GetPlaneVectorIntersection(TVector3 pos, TVector3 dir, TVector3 n, TVector3 a);
+TVector3 GetPlaneVectorIntersection(const TVector3& pos, const TVector3& dir, TVector3 const& n,
+                                    TVector3 const& a);
+
+TMatrixD GetConeMatrix(const TVector3& d, const Double_t& cosTheta);
+
+Double_t GetConeVectorIntersection(const TVector3& pos, const TVector3& dir, const TVector3& d,
+                                   const TVector3& v, const Double_t& cosTheta);
+Double_t GetConeVectorIntersection(const TVector3& pos, const TVector3& dir, const TMatrixD& M,
+                                   const TVector3& axis, const TVector3& v);
 
 Double_t DistanceToAxis(const TVector3& axisPoint, const TVector3& axisVector, const TVector3& point);
 

source/framework/tools/inc/TRestStringHelper.h

@@ -34,14 +34,16 @@ Float_t StringToFloat(std::string in);
 Double_t StringToDouble(std::string in);
 Int_t StringToInteger(std::string in);
 std::string IntegerToString(Int_t n);
+std::string DoubleToString(Double_t d);
 Bool_t StringToBool(std::string in);
 Long64_t StringToLong(std::string in);
 TVector3 StringTo3DVector(std::string in);
 TVector2 StringTo2DVector(std::string in);
 std::vector<std::string> Split(std::string in, std::string separator, bool allowBlankString = false,
                                bool removeWhiteSpaces = false, int startPos = -1);
-std::vector<double> StringToElements(std::string in, std::string separator, bool allowBlankString = false,
-                                     bool removeWhiteSpaces = false, int starPos = -1);
+std::vector<double> StringToElements(std::string in, std::string separator);
+std::vector<double> StringToElements(std::string in, std::string headChar, std::string separator,
+                                     std::string tailChar);
 std::string RemoveWhiteSpaces(std::string in);
 std::string Replace(std::string in, std::string thisString, std::string byThisString, size_t fromPosition = 0,
                     Int_t N = 0);

source/framework/tools/inc/TRestTools.h

@@ -54,13 +54,19 @@ class TRestTools {
                               Int_t skipLines = 0);
 
     template <typename T>
-    static int ReadBinaryTable(std::string fName, std::vector<std::vector<T>>& data, Int_t columns);
+    static int ReadBinaryTable(std::string fName, std::vector<std::vector<T>>& data, Int_t columns = -1);
+
+    static Bool_t IsBinaryFile(std::string fname);
+
+    static std::string GetFileNameExtension(std::string fullname);
+
+    static int GetBinaryFileColumns(std::string fname);
 
     template <typename T>
-    static T GetMaxValueFromTable(std::vector<std::vector<T>> data, Int_t column);
+    static T GetMaxValueFromTable(const std::vector<std::vector<T>>& data, Int_t column = -1);
 
     template <typename T>
-    static T GetMinValueFromTable(std::vector<std::vector<T>> data, Int_t column);
+    static T GetMinValueFromTable(const std::vector<std::vector<T>>& data, Int_t column = -1);
 
     template <typename T>
     static T GetLowestIncreaseFromTable(std::vector<std::vector<T>> data, Int_t column);
@@ -71,6 +77,9 @@ class TRestTools {
     template <typename T>
     static int ExportASCIITable(std::string fname, std::vector<std::vector<T>>& data);
 
+    template <typename T>
+    static int ExportBinaryTable(std::string fname, std::vector<std::vector<T>>& data);
+
     static Int_t isValidFile(const std::string& path);
     static bool fileExists(const std::string& filename);
     static bool isRootFile(const std::string& filename);

source/framework/tools/src/TRestPhysics.cxx

@@ -78,10 +78,104 @@ Double_t DistanceToAxis(const TVector3& axisPoint, const TVector3& axisVector, c
 /// and  moving in direction `dir` and the plane defined by its normal vector `n` and the point `a`. This is
 /// equivalent to move/translate the position `pos` to the plane.
 ///
-TVector3 GetPlaneVectorIntersection(TVector3 pos, TVector3 dir, TVector3 n, TVector3 a) {
+TVector3 GetPlaneVectorIntersection(const TVector3& pos, const TVector3& dir, const TVector3& n,
+                                    const TVector3& a) {
     return MoveToPlane(pos, dir, n, a);
 }
 
+///////////////////////////////////////////////
+/// \brief It returns the cone matrix M = d^T x d - cosTheta^2 x I, extracted from the document
+/// by "David Eberly, Geometric Tools, Redmond WA 98052, Intersection of a Line and a Cone".
+///
+TMatrixD GetConeMatrix(const TVector3& d, const Double_t& cosTheta) {
+    double cAxis[3];
+    d.GetXYZ(cAxis);
+
+    TVectorD coneAxis(3, cAxis);
+
+    TMatrixD M(3, 3);
+    M.Rank1Update(coneAxis, coneAxis);
+
+    double cT2 = cosTheta * cosTheta;
+    TMatrixD gamma(3, 3);
+    gamma.UnitMatrix();
+    gamma *= cT2;
+
+    M -= gamma;
+    return M;
+}
+
+///////////////////////////////////////////////
+/// \brief This method will find the intersection of the trajectory defined by the vector starting at
+/// `pos` and  moving in direction `dir` and the cone defined by its axis vector `d` and the vertex`v`.
+/// The cosine of the angle defining the cone should be also given inside the `cosTheta` argument.
+///
+/// This method will return `t`, which is the value the particle position, `pos`, needs to be displaced
+/// by the vector, `dir`, to get the particle at the surface of the cone. If the particle does not
+/// cross the cone, then the value returned will be zero (no particle displacement).
+//
+/// This method is based on the document by "David Eberly, Geometric Tools, Redmond WA 98052,
+/// Intersection of a Line and a Cone".
+///
+Double_t GetConeVectorIntersection(const TVector3& pos, const TVector3& dir, const TVector3& d,
+                                   const TVector3& v, const Double_t& cosTheta) {
+    TMatrixD M = GetConeMatrix(d, cosTheta);
+    return GetConeVectorIntersection(pos, dir, M, d, v);
+}
+
+///////////////////////////////////////////////
+/// \brief This method will find the intersection of the trajectory defined by the vector starting at `pos`
+/// and moving in direction `dir` and the cone defined by its characteristic matrix `M`, which is built
+/// using the cone axis vector `d` as `d^T x d`,  and the vertex`v`. The resulting TVector3 will be the
+/// position of the particle placed at the cone surface.
+///
+/// This method will return `t`, which is the value the particle position, `pos`, needs to be displaced
+/// by the vector, `dir`, to get the particle at the surface of the cone. If the particle does not
+/// cross the cone, then the value returned will be zero (no particle displacement).
+///
+/// This method is based on the document by "David Eberly, Geometric Tools, Redmond WA 98052,
+/// Intersection of a Line and a Cone".
+///
+Double_t GetConeVectorIntersection(const TVector3& pos, const TVector3& dir, const TMatrixD& M,
+                                   const TVector3& axis, const TVector3& v) {
+    double u[3];
+    dir.GetXYZ(u);
+    TMatrixD U(3, 1, u);
+    TMatrixD Ut(1, 3, u);
+
+    double delta[3];
+    TVector3 deltaV = pos - v;
+    deltaV.GetXYZ(delta);
+    TMatrixD D(3, 1, delta);
+    TMatrixD Dt(1, 3, delta);
+
+    TMatrixD C2 = Ut * M * U;
+    Double_t c2 = C2[0][0];
+
+    TMatrixD C1 = Ut * M * D;
+    Double_t c1 = C1[0][0];
+
+    TMatrixD C0 = Dt * M * D;
+    Double_t c0 = C0[0][0];
+
+    Double_t root = c1 * c1 - c0 * c2;
+    if (root < 0) return 0;
+
+    Double_t t1 = (-c1 + TMath::Sqrt(root)) / c2;
+    Double_t t2 = (-c1 - TMath::Sqrt(root)) / c2;
+
+    // The projections along the cone axis. If positive then the solution
+    // gives the cone intersection with the side defined by `axis`
+    Double_t h1 = t1 * dir.Dot(axis) + axis.Dot(deltaV);
+    Double_t h2 = t2 * dir.Dot(axis) + axis.Dot(deltaV);
+
+    // We use it to select the root we are interested in
+    if (h2 > 0)
+        return t2;
+    else
+        return t1;
+}
+
 ///////////////////////////////////////////////
 /// \brief This method transports a position `pos` by a distance `d` in the direction defined by `dir`.
 ///

source/framework/tools/src/TRestStringHelper.cxx

@@ -187,13 +187,10 @@ std::vector<string> REST_StringHelper::Split(std::string in, string separator, b
 /// \brief Convert the input string into a  vector of double elements
 ///
 /// e.g. Input: "1,2,3,4", Output: {1.,2.,3.,4.}
-std::vector<double> REST_StringHelper::StringToElements(std::string in, string separator,
-                                                        bool allowBlankString, bool removeWhiteSpaces,
-                                                        int startPos) {
-std:
+///
+std::vector<double> REST_StringHelper::StringToElements(std::string in, string separator) {
     vector<double> result;
-    vector<string> vec_str =
-        REST_StringHelper::Split(in, separator, allowBlankString, removeWhiteSpaces, startPos);
+    vector<string> vec_str = REST_StringHelper::Split(in, separator);
     for (unsigned int i = 0; i < vec_str.size(); i++) {
         double temp = REST_StringHelper::StringToDouble(vec_str[i]);
         result.push_back(temp);
@@ -202,6 +199,31 @@ std::vector<double> REST_StringHelper::StringToElements(std::string in, string s
     return result;
 }
 
+///////////////////////////////////////////////
+/// \brief Convert the input string `in` into a vector of double elements
+///
+/// Called as `StringToElements( in, "[", ",", "]" );` will get the
+/// elements from a string with the following format "[a,b,c]" where a,b,c
+/// are double numbers.
+///
+std::vector<double> REST_StringHelper::StringToElements(std::string in, string headChar, string separator,
+                                                        string tailChar) {
+    std::vector<double> result;
+    size_t startPos = in.find(headChar);
+    size_t endPos = in.find(tailChar);
+    if (startPos == string::npos || endPos == string::npos) {
+        return result;
+    }
+    std::vector<string> values = Split(in.substr(startPos + 1, endPos - startPos - 1), ",");
+
+    for (unsigned int i = 0; i < values.size(); i++) {
+        double temp = REST_StringHelper::StringToDouble(values[i]);
+        result.push_back(temp);
+    }
+
+    return result;
+}
+
 ///////////////////////////////////////////////
 /// \brief Returns the input string removing all white spaces.
 ///
@@ -475,6 +497,11 @@ Int_t REST_StringHelper::StringToInteger(string in) {
 ///
 string REST_StringHelper::IntegerToString(Int_t n) { return Form("%d", n); }
 
+///////////////////////////////////////////////
+/// \brief Gets a string from a double
+///
+string REST_StringHelper::DoubleToString(Double_t d) { return Form("%4.2lf", d); }
+
 Bool_t REST_StringHelper::StringToBool(std::string in) {
     return (ToUpper(in) == "TRUE" || ToUpper(in) == "ON");
 }