diff --git a/cmake/Modules/LauraExternalDependencies.cmake b/cmake/Modules/LauraExternalDependencies.cmake index ac62c81..b973cee 100644 --- a/cmake/Modules/LauraExternalDependencies.cmake +++ b/cmake/Modules/LauraExternalDependencies.cmake @@ -1,32 +1,34 @@ # Laura++ needs ROOT, we just require that it is a version that exports CMake targets properly # We may depend on other things in the future, e.g. JSON or YAML for easier configuration of particle properties etc. if(DEFINED ENV{ROOTSYS}) list(APPEND CMAKE_PREFIX_PATH $ENV{ROOTSYS}) endif() if(LAURA_BUILD_ROOFIT_TASK) find_package(ROOT 6.14 REQUIRED COMPONENTS EG RooFitCore RooFit) else() find_package(ROOT 6.14 REQUIRED COMPONENTS EG RooFitCore) endif() +find_package(nlohmann_json REQUIRED) + #message(STATUS "ROOT include directories: ${ROOT_INCLUDE_DIRS}") #message(STATUS "ROOT libraries: ${ROOT_LIBRARIES}") #message(STATUS "ROOT definitions: ${ROOT_DEFINITIONS}") #message(STATUS "ROOT CXX flags: ${ROOT_CXX_FLAGS}") #message(STATUS "ROOT CC flags: ${ROOT_CC_FLAGS}") #message(STATUS "ROOT use file: ${ROOT_USE_FILE}") # Don't want to do this because it uses old-style CMake #include(${ROOT_USE_FILE}) if(EXISTS "${ROOT_DIR}/RootMacros.cmake") message(STATUS "Laura++: Including ROOT macros module: ${ROOT_DIR}/RootMacros.cmake") include(${ROOT_DIR}/RootMacros.cmake) elseif(EXISTS "${ROOT_DIR}/modules/RootNewMacros.cmake") message(STATUS "Laura++: Including ROOT macros module: ${ROOT_DIR}/modules/RootNewMacros.cmake") include(${ROOT_DIR}/modules/RootNewMacros.cmake) else() message(WARNING "Laura++: Cannot locate ROOT macros module in ${ROOT_DIR}") endif() diff --git a/inc/Lau1DCubicSpline.hh b/inc/Lau1DCubicSpline.hh index adb0544..dbfff33 100644 --- a/inc/Lau1DCubicSpline.hh +++ b/inc/Lau1DCubicSpline.hh @@ -1,277 +1,298 @@ /* Copyright 2015 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file Lau1DCubicSpline.hh \brief File containing declaration of Lau1DCubicSpline class. */ /*! \class Lau1DCubicSpline \brief Class for defining a 1D cubic spline based on a set of knots. Class for defining a 1D cubic spline based on a set of knots. Interpolation between the knots is performed either by one of two types of cubic spline (standard or Akima) or by linear interpolation. The splines are defined by a piecewise cubic function which, between knots i and i+1, has the form f_i(x) = (1 - t)*y_i + t*y_i+1 + t*(1 - t)(a*(1 - t) + b*t) where t = (x - x_i)/(x_i+1 - x_i), a = k_i *(x_i+1 - x_i) - (y_i+1 - y_i), b = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i) and k_i is (by construction) the first derivative at knot i. f(x) and f'(x) are continuous at the internal knots by construction. For the standard splines, f''(x) is required to be continuous at all internal knots placing n-2 constraints on the n parameters, k_i. The final two constraints are set by the boundary conditions. At each boundary, the function may be: (i) Clamped : f'(x) = C at the last knot (ii) Natural : f''(x) = 0 at the last knot (iii) Not a knot : f'''(x) continuous at the second last knot The algorithms used in these splines can be found on: http://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm For the Akima splines, the values of k_i are determined from the slopes of the four nearest segments (a_i-1, a_i, a_i+1 and a_i+2) as k_i = ( | a_i+2 - a_i+1 | a_i + | a_i - a_i-1 | a_i+1 ) / ( | a_i+2 - a_i+1 | + | a_i - a_i-1 | ) and as k_i = ( a_i + a_i+1 ) / 2 in the special case a_i-1 == a_i != a_i+1 == a_i+2. Boundary conditions are specified by the relations a_2 - a_1 = a_1 - a_0 = a_0 - a_-1 and a_n-1 - a_n-2 = a_n - a_n-1 = a_n+1 - a_n The algorithms used in these splines can be found in: J.ACM vol. 17 no. 4 pp 589-602 */ #ifndef LAU_1DCUBICSPLINE #define LAU_1DCUBICSPLINE #include #include +#include + #include "Rtypes.h" #include "TF1.h" #include "TFitResultPtr.h" #include "LauPrint.hh" class TH1; class LauAbsRValue; class LauParameter; class Lau1DCubicSpline final { public: //! Define the allowed interpolation types enum class SplineType { StandardSpline, /*!< standard cubic splines with f''(x) continuous at all internal knots */ AkimaSpline, /*!< Akima cubic splines with f'(x) at each knot defined locally by the positions of only five knots */ LinearInterpolation /*! Linear interpolation between each pair of knots */ }; //! Define the allowed boundary condition types /*! These are only supported by standard splines */ enum class BoundaryType { Clamped, /*!< clamped boundary - f'(x) = C */ Natural, /*!< natural boundary - f''(x) = 0 */ NotAKnot /*!< 'not a knot' boundary - f'''(x) continuous at second last knot */ }; //! Constructor /*! \param [in] xs the x-values of the knots \param [in] ys the y-values of the knots \param [in] type the type of spline (Standard, Akima, Linear) \param [in] leftBound the left-hand boundary condition \param [in] rightBound the right-hand boundary condition \param [in] dydx0 the gradient at the left-hand boundary of a clamped spline \param [in] dydxn the gradient at the right-hand boundary of a clamped spline */ Lau1DCubicSpline(const std::vector& xs, const std::vector& ys, const SplineType type = SplineType::StandardSpline, const BoundaryType leftBound = BoundaryType::NotAKnot, const BoundaryType rightBound = BoundaryType::NotAKnot, const Double_t dydx0 = 0.0, const Double_t dydxn = 0.0); //! Evaluate the function at given point /*! \param [in] x the x-coordinate \return the value of the spline at x */ Double_t evaluate(const Double_t x) const; //! Update the y-values of the knots /*! \param [in] ys the y-values of the knots */ void updateYValues(const std::vector& ys); //! Update the y-values of the knots /*! \param [in] ys the y-values of the knots */ void updateYValues(const std::vector& ys); //! Update the y-values of the knots /*! \param [in] ys the y-values of the knots */ void updateYValues(const std::vector& ys); //! Update the type of interpolation to perform /*! \param [in] type the type of interpolation */ void updateType(const SplineType type); //! Update the boundary conditions for the spline /*! \param [in] leftBound the left-hand boundary condition \param [in] rightBound the right-hand boundary condition \param [in] dydx0 the gradient at the left-hand boundary of a clamped spline \param [in] dydxn the gradient at the right-hand boundary of a clamped spline */ void updateBoundaryConditions(const BoundaryType leftBound, const BoundaryType rightBound, const Double_t dydx0 = 0.0, const Double_t dydxn = 0.0); //! Get the number of knots std::size_t getnKnots() const {return nKnots_;} //! Get y values const std::vector& getYValues() const {return y_;} //! Get x values const std::vector& getXValues() const {return x_;} //! Get the coefficients of a given spline segment in the form a + bx + cx^2 + dx^3 /*! \param [in] segIndex refers to the index of the knot at the left end of the segment (segIndex = 0 gets the coefficients of the the segment between x_0 and x_1) \param [in] normalise if true, the coefficients returned will be normalised by the integral of the complete spline (defaults to false) \return coefficients {a, b, c, d} */ std::array getCoefficients(const std::size_t segIndex, const bool normalise = false) const; //! Get the integral over all the spline segments Double_t integral() const; //! Make a TF1 showing the spline with its current knot values /*! \param [in] normalise whether or not you want the spline normalised \return 1D function object */ TF1* makeTF1(const bool normalise = false) const; //! Fit the the normalisation of the spline to a TH1 /*! \param [in] hist the histogram to be fitted \param [in] printLevel the level of printout desired from fit \return a TF1 fit to the histogram */ TF1* normaliseToTH1(TH1* hist, const LauOutputLevel printLevel = LauOutputLevel::Standard) const; //! Fit the spline to a TH1 /*! \param [in] hist the histogram to be fitted \param [in] printLevel the level of printout desired from fit \param [in] doWL If true do a weighted likelihood fit, else chi2 \return a shared-ownership smart pointer to the fit results */ TFitResultPtr fitToTH1(TH1* hist, const LauOutputLevel printLevel = LauOutputLevel::Standard, const bool doWL = false); //! Set spline variables to those read from a json file /* \param[in] filename to be used \param[in] name of spline, if empty no name is used \return the newly constructed spline */ //static Lau1DCubicSpline readFromJson(const TString filename, const TString splineName = ""); //! Write the spline to a json file using nlohmann json /*! \param[in] filename to be used \param[in] name of spline, if empty, no name is used \param[in] if true and filename already exists, append the spline to the file */ void writeToJson(const TString filename, const TString splineName = "", const bool append = false) const; private: //! Initialise the class void init(); //! Calculate the first derivative at each knot void calcDerivatives(); //! Calculate the first derivatives according to the standard method void calcDerivativesStandard(); //! Calculate the first derivatives according to the Akima method void calcDerivativesAkima(); //! The number of knots in the spline - const std::size_t nKnots_; + std::size_t nKnots_; //! The x-value at each knot - const std::vector x_; + std::vector x_; //! The y-value at each knot std::vector y_; //! The first derivative at each knot std::vector dydx_; //! The 'a' coefficients used to determine the derivatives std::vector a_; //! The 'b' coefficients used to determine the derivatives std::vector b_; //! The 'c' coefficients used to determine the derivatives std::vector c_; //! The 'd' coefficients used to determine the derivatives std::vector d_; //! The type of interpolation to be performed SplineType type_; //! The left-hand boundary condition on the spline BoundaryType leftBound_; //! The right-hand boundary condition on the spline BoundaryType rightBound_; //! The gradient at the left boundary for a clamped spline Double_t dydx0_; //! The gradient at the right boundary for a clamped spline Double_t dydxn_; + // JSON serialisation stuff + + // map SplineType values to JSON as strings + NLOHMANN_JSON_SERIALIZE_ENUM( SplineType, { + {SplineType::StandardSpline, "StandardSpline"}, + {SplineType::AkimaSpline, "AkimaSpline"}, + {SplineType::LinearInterpolation, "LinearInterpolation"}, + }) + + // map BoundaryType values to JSON as strings + NLOHMANN_JSON_SERIALIZE_ENUM( BoundaryType, { + {BoundaryType::Clamped, "Clamped"}, + {BoundaryType::Natural, "Natural"}, + {BoundaryType::NotAKnot, "NotAKnot"}, + }) + + // enable serialisation of this class + NLOHMANN_DEFINE_TYPE_INTRUSIVE(Lau1DCubicSpline, nKnots_, x_, y_, type_, leftBound_, rightBound_, dydx0_, dydxn_) + ClassDef(Lau1DCubicSpline, 0) // Class for defining a 1D cubic spline }; #endif diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 31fca2d..1347550 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,54 +1,56 @@ # Set the include directory (unfortunately this old-style stuff is necessary for the ROOT_GENERATE_DICTIONARY macro) include_directories(${PROJECT_SOURCE_DIR}/inc) +get_target_property(NLOHMANN_JSON_INCLUDEDIR nlohmann_json::nlohmann_json INTERFACE_INCLUDE_DIRECTORIES) +include_directories(${NLOHMANN_JSON_INCLUDEDIR}) # Use glob to find the headers and sources file(GLOB LAURA_HEADERS ${PROJECT_SOURCE_DIR}/inc/*.hh) file(GLOB LAURA_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/*.cc) if (NOT LAURA_BUILD_ROOFIT_TASK) list(REMOVE_ITEM LAURA_HEADERS ${PROJECT_SOURCE_DIR}/inc/LauRooFitTask.hh) list(REMOVE_ITEM LAURA_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/LauRooFitTask.cc) endif() # Generate the rootcint file set(LAURA_LINKDEF ${PROJECT_SOURCE_DIR}/inc/Laura++_LinkDef.h) set(LAURA_DICTIONARY_ROOT G__Laura++) set(LAURA_DICTIONARY ${LAURA_DICTIONARY_ROOT}.cxx) if (LAURA_BUILD_ROOFIT_TASK) ROOT_GENERATE_DICTIONARY( ${LAURA_DICTIONARY_ROOT} ${LAURA_HEADERS} LINKDEF ${LAURA_LINKDEF} OPTIONS -DDOLAUROOFITTASK ) else() ROOT_GENERATE_DICTIONARY( ${LAURA_DICTIONARY_ROOT} ${LAURA_HEADERS} LINKDEF ${LAURA_LINKDEF} ) endif() # Build the shared library add_library(Laura++ SHARED ${LAURA_SOURCES} ${LAURA_DICTIONARY}) set_target_properties(Laura++ PROPERTIES OUTPUT_NAME Laura++) set_target_properties(Laura++ PROPERTIES VERSION ${CMAKE_PROJECT_VERSION} SOVERSION ${CMAKE_PROJECT_VERSION_MAJOR}) set_target_properties(Laura++ PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR}) target_include_directories(Laura++ PUBLIC $ $) -target_link_libraries(Laura++ ROOT::Core ROOT::Hist ROOT::Matrix ROOT::Physics ROOT::Minuit ROOT::EG ROOT::Tree ROOT::RooFitCore) +target_link_libraries(Laura++ nlohmann_json::nlohmann_json ROOT::Core ROOT::Hist ROOT::Matrix ROOT::Physics ROOT::Minuit ROOT::EG ROOT::Tree ROOT::RooFitCore) if (LAURA_BUILD_ROOFIT_TASK) target_link_libraries(Laura++ ROOT::RooFit) endif() # Install the libraries install( TARGETS Laura++ EXPORT "LauraTargets" LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} ) # Install the pcm and rootmap files generated by ROOT_GENERATE_DICTIONARY install( FILES ${CMAKE_CURRENT_BINARY_DIR}/libLaura++.rootmap ${CMAKE_CURRENT_BINARY_DIR}/libLaura++_rdict.pcm DESTINATION ${CMAKE_INSTALL_LIBDIR} )