diff --git a/cmake/MacroRootDict.cmake b/cmake/MacroRootDict.cmake index 2422372d4..ba51925dd 100644 --- a/cmake/MacroRootDict.cmake +++ b/cmake/MacroRootDict.cmake @@ -519,4 +519,10 @@ MACRO(COMPILELIB dependency) set(dirs_included ${dirs_to_include} PARENT_SCOPE) set(library_added ${libname}) set(library_added ${library_added} PARENT_SCOPE) -ENDMACRO() \ No newline at end of file + + # define REST_*Lib e.g. REST_DetectorLib using library name: RestDetector -> REST_DetectorLib + string(REGEX REPLACE "^Rest(.+)$" "REST_\\1Lib" DEFINE_VARIABLE_NAME ${libname}) + message(STATUS "Adding compile definition: ${DEFINE_VARIABLE_NAME}") + add_compile_definitions(${DEFINE_VARIABLE_NAME}) + +ENDMACRO() diff --git a/cmake/Testing.cmake b/cmake/Testing.cmake index 41fe65f74..12a40b7fa 100644 --- a/cmake/Testing.cmake +++ b/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}") diff --git a/source/framework/tools/inc/TRestPhysics.h b/source/framework/tools/inc/TRestPhysics.h index 4f4938066..c0f68c177 100644 --- a/source/framework/tools/inc/TRestPhysics.h +++ b/source/framework/tools/inc/TRestPhysics.h @@ -25,8 +25,10 @@ #include +#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); diff --git a/source/framework/tools/inc/TRestStringHelper.h b/source/framework/tools/inc/TRestStringHelper.h index 2dfa9c8af..c143eb229 100644 --- a/source/framework/tools/inc/TRestStringHelper.h +++ b/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 Split(std::string in, std::string separator, bool allowBlankString = false, bool removeWhiteSpaces = false, int startPos = -1); -std::vector StringToElements(std::string in, std::string separator, bool allowBlankString = false, - bool removeWhiteSpaces = false, int starPos = -1); +std::vector StringToElements(std::string in, std::string separator); +std::vector 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); diff --git a/source/framework/tools/inc/TRestTools.h b/source/framework/tools/inc/TRestTools.h index 460e01586..7c80306af 100644 --- a/source/framework/tools/inc/TRestTools.h +++ b/source/framework/tools/inc/TRestTools.h @@ -54,13 +54,19 @@ class TRestTools { Int_t skipLines = 0); template - static int ReadBinaryTable(std::string fName, std::vector>& data, Int_t columns); + static int ReadBinaryTable(std::string fName, std::vector>& 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 - static T GetMaxValueFromTable(std::vector> data, Int_t column); + static T GetMaxValueFromTable(const std::vector>& data, Int_t column = -1); template - static T GetMinValueFromTable(std::vector> data, Int_t column); + static T GetMinValueFromTable(const std::vector>& data, Int_t column = -1); template static T GetLowestIncreaseFromTable(std::vector> data, Int_t column); @@ -71,6 +77,9 @@ class TRestTools { template static int ExportASCIITable(std::string fname, std::vector>& data); + template + static int ExportBinaryTable(std::string fname, std::vector>& data); + static Int_t isValidFile(const std::string& path); static bool fileExists(const std::string& filename); static bool isRootFile(const std::string& filename); diff --git a/source/framework/tools/src/TRestPhysics.cxx b/source/framework/tools/src/TRestPhysics.cxx index 7c3c87341..f9d06640a 100644 --- a/source/framework/tools/src/TRestPhysics.cxx +++ b/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`. /// diff --git a/source/framework/tools/src/TRestStringHelper.cxx b/source/framework/tools/src/TRestStringHelper.cxx index 73dbbb0d7..b8ea5b19a 100644 --- a/source/framework/tools/src/TRestStringHelper.cxx +++ b/source/framework/tools/src/TRestStringHelper.cxx @@ -187,13 +187,10 @@ std::vector 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 REST_StringHelper::StringToElements(std::string in, string separator, - bool allowBlankString, bool removeWhiteSpaces, - int startPos) { -std: +/// +std::vector REST_StringHelper::StringToElements(std::string in, string separator) { vector result; - vector vec_str = - REST_StringHelper::Split(in, separator, allowBlankString, removeWhiteSpaces, startPos); + vector 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 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 REST_StringHelper::StringToElements(std::string in, string headChar, string separator, + string tailChar) { + std::vector result; + size_t startPos = in.find(headChar); + size_t endPos = in.find(tailChar); + if (startPos == string::npos || endPos == string::npos) { + return result; + } + std::vector 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"); } diff --git a/source/framework/tools/src/TRestTools.cxx b/source/framework/tools/src/TRestTools.cxx index 574c8aeef..c82fe1bd1 100644 --- a/source/framework/tools/src/TRestTools.cxx +++ b/source/framework/tools/src/TRestTools.cxx @@ -160,7 +160,34 @@ template int TRestTools::ExportASCIITable(std::string fname, std::vector>& data); /////////////////////////////////////////////// -/// \brief Reads a binary file containing a fixed-columns table with values +/// \brief Writes the contents of the vector table given as argument to `fname` as a binary file. +/// Allowed types are Int_t, Float_t and Double_t. +/// +template +int TRestTools::ExportBinaryTable(std::string fname, std::vector>& data) { + ofstream file(fname, ios::out | ios::binary); + if (!file.is_open()) { + ferr << "Unable to open file for writting : " << fname << endl; + return 1; + } + + for (int n = 0; n < data.size(); n++) + for (int m = 0; m < data[n].size(); m++) { + file.write((char*)&data[n][m], sizeof(T)); + } + file.close(); + + return 0; +} + +template int TRestTools::ExportBinaryTable(std::string fname, std::vector>& data); +template int TRestTools::ExportBinaryTable(std::string fname, + std::vector>& data); +template int TRestTools::ExportBinaryTable(std::string fname, + std::vector>& data); + +/////////////////////////////////////////////// +/// \brief Reads a binary file containning a fixed-columns table with values /// /// This method will open the file fName. This file should contain a /// table with numeric values of the type specified inside the syntax < >. @@ -169,7 +196,7 @@ template int TRestTools::ExportASCIITable(std::string fname, /// /// \code /// std::vector > fvec; -/// ReadBinaryFile( "myfile.bin", fvec, 6); +/// ReadBinaryTable( "myfile.bin", fvec, 6); /// \endcode /// /// The values on the table will be loaded in the matrix provided through the @@ -178,11 +205,20 @@ template int TRestTools::ExportASCIITable(std::string fname, template int TRestTools::ReadBinaryTable(string fName, std::vector>& data, Int_t columns) { if (!TRestTools::isValidFile((string)fName)) { - cout << "TRestTools::ReadBinaryTable. Error." << endl; - cout << "Cannot open file : " << fName << endl; + ferr << "TRestTools::ReadBinaryTable. Error." << endl; + ferr << "Cannot open file : " << fName << endl; return 0; } + if (columns == -1) { + columns = GetBinaryFileColumns(fName); + if (columns == -1) { + ferr << "TRestTools::ReadBinaryTable. Format extension error." << endl; + ferr << "Please, specify the number of columns at the method 3rd argument" << endl; + return 0; + } + } + std::ifstream fin(fName, std::ios::binary); fin.seekg(0, std::ios::end); const size_t num_elements = fin.tellg() / sizeof(T); @@ -197,9 +233,10 @@ int TRestTools::ReadBinaryTable(string fName, std::vector>& data, } std::vector dataArray(columns); + fin.read(reinterpret_cast(&dataArray[0]), columns * sizeof(T)); while (fin.good()) { - fin.read(reinterpret_cast(&dataArray[0]), columns * sizeof(T)); data.push_back(dataArray); + fin.read(reinterpret_cast(&dataArray[0]), columns * sizeof(T)); } return 1; } @@ -211,50 +248,118 @@ template int TRestTools::ReadBinaryTable(string fName, std::vector(string fName, std::vector>& data, Int_t columns); +/////////////////////////////////////////////// +/// \brief It identifies if the filename extension follows the formatting ".Nxyzf", where the +/// number of columns is `xyz`, and the last character is the type of data (f/d/i), float, +/// double and integer respectively. +/// +Bool_t TRestTools::IsBinaryFile(string fname) { + if (GetBinaryFileColumns(fname) > 0) return true; + return false; +} + +/////////////////////////////////////////////// +/// \brief It extracts the number of columns from the filename extension given by argument. +/// The file will containing a binary formatted table with a fixed number of rows and columns. +/// +/// The filename extension will be : ".Nxyzf", where the number of columns is `xyz`, and the +/// last character is the type of data (f/d/i), float, double and integer respectively. +/// +int TRestTools::GetBinaryFileColumns(string fname) { + string extension = GetFileNameExtension(fname); + if (extension.find("N") != 0) { + ferr << "Wrong filename extension." << endl; + ferr << "Cannot guess the number of columns" << endl; + return -1; + } + + size_t pos = extension.find("i"); + if (pos != string::npos) { + return StringToInteger(extension.substr(1, pos - 1)); + } + + pos = extension.find("f"); + if (pos != string::npos) { + return StringToInteger(extension.substr(1, pos - 1)); + } + + pos = extension.find("d"); + if (pos != string::npos) { + return StringToInteger(extension.substr(1, pos - 1)); + } + + ferr << "Format " << ToUpper(extension) << " not recognized" << endl; + return -1; +} + /////////////////////////////////////////////// /// \brief It returns the maximum value for a particular `column` from the table given by -/// argument. +/// argument. If no column is specified in the arguments, then it gets the maximum from the +/// full table. /// /// This method is available for tables of type Float_t, Double_t and Int_t. /// template -T TRestTools::GetMaxValueFromTable(std::vector> data, Int_t column) { - if (data.size() == 0 || data[0].size() <= column) return 0; - T maxValue = data[0][column]; - for (int n = 0; n < data.size(); n++) - if (maxValue < data[n][column]) maxValue = data[n][column]; +T TRestTools::GetMaxValueFromTable(const std::vector>& data, Int_t column) { + if (data.size() == 0) return 0; + if (column != -1 && data[0].size() <= column) return 0; + + T maxValue = data[0][0]; + if (column == -1) { + for (int n = 0; n < data.size(); n++) + for (int c = 0; c < data[n].size(); c++) + if (maxValue < data[n][c]) maxValue = data[n][c]; + } else { + maxValue = data[0][column]; + for (int n = 0; n < data.size(); n++) + if (maxValue < data[n][column]) maxValue = data[n][column]; + } + return maxValue; } -template Int_t TRestTools::GetMaxValueFromTable(std::vector> data, Int_t column); +template Int_t TRestTools::GetMaxValueFromTable(const std::vector>& data, + Int_t column); -template Float_t TRestTools::GetMaxValueFromTable(std::vector> data, +template Float_t TRestTools::GetMaxValueFromTable(const std::vector>& data, Int_t column); -template Double_t TRestTools::GetMaxValueFromTable(std::vector> data, +template Double_t TRestTools::GetMaxValueFromTable(const std::vector>& data, Int_t column); /////////////////////////////////////////////// /// \brief It returns the minimum value for a particular `column` from the table given by -/// argument. +/// argument. If no column is specified in the arguments, then it gets the minimum from the +/// full table. /// /// This method is available for tables of type Float_t, Double_t and Int_t. /// template -T TRestTools::GetMinValueFromTable(std::vector> data, Int_t column) { - if (data.size() == 0 || data[0].size() <= column) return 0; - T minValue = data[0][column]; - for (int n = 0; n < data.size(); n++) - if (minValue > data[n][column]) minValue = data[n][column]; +T TRestTools::GetMinValueFromTable(const std::vector>& data, Int_t column) { + if (data.empty()) return 0; + if (column != -1 && data[0].size() <= column) return 0; + + T minValue = data[0][0]; + if (column == -1) { + for (int n = 0; n < data.size(); n++) + for (int c = 0; c < data[n].size(); c++) + if (minValue > data[n][c]) minValue = data[n][c]; + } else { + minValue = data[0][column]; + for (int n = 0; n < data.size(); n++) + if (minValue > data[n][column]) minValue = data[n][column]; + } + return minValue; } -template Int_t TRestTools::GetMinValueFromTable(std::vector> data, Int_t column); +template Int_t TRestTools::GetMinValueFromTable(const std::vector>& data, + Int_t column); -template Float_t TRestTools::GetMinValueFromTable(std::vector> data, +template Float_t TRestTools::GetMinValueFromTable(const std::vector>& data, Int_t column); -template Double_t TRestTools::GetMinValueFromTable(std::vector> data, +template Double_t TRestTools::GetMinValueFromTable(const std::vector>& data, Int_t column); /////////////////////////////////////////////// @@ -487,6 +592,20 @@ std::pair TRestTools::SeparatePathAndName(string fullname) { return result; } +/////////////////////////////////////////////// +/// \brief Gets the file extension as the substring found after the lastest "." +/// +/// Input: "/home/jgalan/abc.txt" Output: "txt" +/// +string TRestTools::GetFileNameExtension(string fullname) { + int pos = fullname.find_last_of('.', -1); + + if (pos != -1) { + return fullname.substr(pos + 1, fullname.size() - pos - 1); + } + return fullname; +} + /////////////////////////////////////////////// /// \brief Returns the input string but without multiple slashes ("/") /// @@ -766,6 +885,8 @@ std::string TRestTools::DownloadRemoteFile(string url) { if (out == 1024) { warning << "Retrying download in 5 seconds" << endl; std::this_thread::sleep_for(std::chrono::seconds(5)); + } else if (attempts < 10) { + success << "Download suceeded after " << 10 - attempts << " attempts" << endl; } attempts--; } while (out == 1024 && attempts > 0);