diff --git a/inc/Lau1DCubicSpline.hh b/inc/Lau1DCubicSpline.hh index cb039fe..633715f 100644 --- a/inc/Lau1DCubicSpline.hh +++ b/inc/Lau1DCubicSpline.hh @@ -1,260 +1,251 @@ /* 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 "LauParameter.hh" #include "Rtypes.h" #include "TF1.h" -class Lau1DCubicSpline { +class Lau1DCubicSpline final { public: //! Define the allowed interpolation types enum LauSplineType { 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 LauSplineBoundaryType { 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] 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, LauSplineType type = Lau1DCubicSpline::StandardSpline, LauSplineBoundaryType leftBound = Lau1DCubicSpline::NotAKnot, LauSplineBoundaryType rightBound = Lau1DCubicSpline::NotAKnot, Double_t dydx0 = 0.0, Double_t dydxn = 0.0); - //! Destructor - virtual ~Lau1DCubicSpline(); - //! Evaluate the function at given point /*! \param [in] x the x-coordinate \return the value of the spline at x */ Double_t evaluate(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(LauSplineType 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(LauSplineBoundaryType leftBound, LauSplineBoundaryType rightBound, Double_t dydx0 = 0.0, Double_t dydxn = 0.0); //! Return the number of knots UInt_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 spline section i in the form a + bx + cx^2 + dx^3 /*! \params [in] i refers to the left-hand index of the knot (i = 0 gets the params between x_0 and x_1) \return coefficients {a, b, c, d} */ std::array getCoefficients(const UInt_t i, 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 /*! \params [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 //! So they look as good as possible when plotted on top of one another //! Useful when a sample has been generated with a hist and fitted with a spline to compare them /*! \params [in] The histogram to be fit to \return a TF1 fit to `hist` */ TF1* normaliseToTH1(TH1* hist) const; //! Fit the spline to a TH1 //! Useful as a starting-point for fitting a spline to data (maybe the hist is taken from MC) /*! \params [in] The histogram to be fit to */ void fitToTH1(TH1* hist); private: - //! Copy constructor - not implemented - Lau1DCubicSpline( const Lau1DCubicSpline& rhs ); - - //! Copy assignment operator - not implemented - Lau1DCubicSpline& operator=(const Lau1DCubicSpline& rhs); - //! 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 UInt_t nKnots_; //! The x-value at each knot 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 LauSplineType type_; //! The left-hand boundary condition on the spline LauSplineBoundaryType leftBound_; //! The right-hand boundary condition on the spline LauSplineBoundaryType 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_; ClassDef(Lau1DCubicSpline, 0); // Class for defining a 1D cubic spline }; #endif diff --git a/inc/LauAbsDecayTimeEfficiency.hh b/inc/LauAbsDecayTimeEfficiency.hh index 8be89db..1c2edb2 100644 --- a/inc/LauAbsDecayTimeEfficiency.hh +++ b/inc/LauAbsDecayTimeEfficiency.hh @@ -1,91 +1,93 @@ /* Copyright 2021 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 LauAbsDecayTimeEfficiency.hh \brief File containing declaration of LauAbsDecayTimeEfficiency class. */ /*! \class LauAbsDecayTimeEfficiency \brief Class for defining the abstract interface for modelling decay time efficiency */ #ifndef LAU_ABS_DECAYTIME_EFFICIENCY #define LAU_ABS_DECAYTIME_EFFICIENCY #include #include "Rtypes.h" class LauAbsRValue; class LauAbsDecayTimeEfficiency { public: //! Retrieve the efficiency for a given value of the decay time /*! \param abscissa the value of the decay time \return the efficiency */ virtual Double_t getEfficiency( const Double_t abscissa ) const = 0; //! Retrieve the parameters of the efficiency model so that they can be loaded into a fit /*! \return the parameters of the efficiency model */ virtual std::vector getParameters() = 0; //! Initialise the parameter cache /*! Must be called prior to starting fitting or generation */ virtual void initialise() = 0; //! Propagate any updates to parameters and recalculate information as neeeded /*! Should be called at each fit iteration */ virtual void propagateParUpdates() = 0; //! Retrieve whether any of the parameters have changed in the latest fit iteration virtual const bool& anythingChanged() const = 0; //! Constructor LauAbsDecayTimeEfficiency() = default; //! Destructor virtual ~LauAbsDecayTimeEfficiency() = default; - //! Copy constructor (deleted) - LauAbsDecayTimeEfficiency(const LauAbsDecayTimeEfficiency& other) = delete; - //! Copy assignment operator (deleted) - LauAbsDecayTimeEfficiency& operator=(const LauAbsDecayTimeEfficiency& other) = delete; //! Move constructor (default) LauAbsDecayTimeEfficiency(LauAbsDecayTimeEfficiency&& other) = default; //! Move assignment operator (default) LauAbsDecayTimeEfficiency& operator=(LauAbsDecayTimeEfficiency&& other) = default; + protected: + //! Copy constructor (default but protected) + LauAbsDecayTimeEfficiency(const LauAbsDecayTimeEfficiency& other) = default; + //! Copy assignment operator (default but protected) + LauAbsDecayTimeEfficiency& operator=(const LauAbsDecayTimeEfficiency& other) = default; + private: ClassDef(LauAbsDecayTimeEfficiency, 0) }; #endif diff --git a/inc/LauSplineDecayTimeEfficiency.hh b/inc/LauSplineDecayTimeEfficiency.hh index 4cde051..ba7cb46 100644 --- a/inc/LauSplineDecayTimeEfficiency.hh +++ b/inc/LauSplineDecayTimeEfficiency.hh @@ -1,152 +1,375 @@ /* Copyright 2021 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 LauSplineDecayTimeEfficiency.hh \brief File containing declaration of LauSplineDecayTimeEfficiency class. */ /*! \class LauSplineDecayTimeEfficiency \brief Class for defining a spline-interpolated model for decay time efficiency */ #ifndef LAU_SPLINE_DECAYTIME_EFFICIENCY #define LAU_SPLINE_DECAYTIME_EFFICIENCY #include #include #include #include "Rtypes.h" +#include "TSystem.h" #include "Lau1DCubicSpline.hh" #include "LauAbsDecayTimeEfficiency.hh" #include "LauParameter.hh" //! Defines the allowed orders of spline polynomials enum class LauSplineOrder : std::size_t { Cubic = 3, //!< 3rd order Sixth = 6 //!< 6th order }; template -class LauAbsSplineDecayTimeEfficiency : public LauAbsDecayTimeEfficiency { +class LauSplineDecayTimeEfficiency : public LauAbsDecayTimeEfficiency { public: - //! The number of coefficients of each spline segment - static constexpr std::size_t nCoeffs { static_cast>(Order) + 1 }; - - private: - ClassDefOverride(LauAbsSplineDecayTimeEfficiency, 0) -}; - -template -class LauSplineDecayTimeEfficiency : public LauAbsSplineDecayTimeEfficiency -{}; + //! Constructor for a Cubic spline + /*! + \param [in] prefix the prefix to use when forming the names of the knot parameters + \param [in] effSpline the cubic spline to be used to represent the efficiency variation + */ + template = true> + LauSplineDecayTimeEfficiency( const TString& prefix, std::unique_ptr effSpline ) : + prefix_{prefix}, + effSpline_{std::move(effSpline)} + { + this->initialSetup(); + } + + //! Constructor for a Sixth order spline + /*! + \param [in] prefix the prefix to use when forming the names of the knot parameters + \param [in] effSpline the cubic spline being used to represent the efficiency variation in another channel + \param [in] correctionSpline the cubic spline to be used to represent the correction to effSpline to obtain the efficiency variation for this channel + */ + template = true> + LauSplineDecayTimeEfficiency( const TString& prefix, const LauSplineDecayTimeEfficiency& effSpline, std::unique_ptr correctionSpline ) : + prefix_{prefix}, + effSpline_{std::move(correctionSpline)}, + otherSpline_{std::make_unique>(effSpline)} + { + this->initialSetup(); + } + + //! Destructor + ~LauSplineDecayTimeEfficiency() = default; + //! Move constructor (default) + LauSplineDecayTimeEfficiency(LauSplineDecayTimeEfficiency&& other) = default; + //! Move assignment operator (default) + LauSplineDecayTimeEfficiency& operator=(LauSplineDecayTimeEfficiency&& other) = default; + //! Copy constructor (custom) + LauSplineDecayTimeEfficiency(const LauSplineDecayTimeEfficiency& other); + //! Copy assignment operator (deleted) + LauSplineDecayTimeEfficiency& operator=(const LauSplineDecayTimeEfficiency& other) = delete; -template <> -class LauSplineDecayTimeEfficiency : public LauAbsSplineDecayTimeEfficiency { - public: - //! Constructor - LauSplineDecayTimeEfficiency( std::unique_ptr effSpline ); + //! The number of coefficients of each spline segment + static constexpr std::size_t nCoeffs { static_cast>(Order) + 1 }; //! Retrieve the efficiency for a given value of the decay time /*! \param abscissa the value of the decay time \return the efficiency */ Double_t getEfficiency( const Double_t abscissa ) const override; //! Retrieve the parameters of the efficiency model so that they can be loaded into a fit /*! \return the parameters of the efficiency model */ std::vector getParameters() override { return params_; } //! Initialise the parameter cache /*! Must be called prior to starting fitting or generation */ void initialise() override; //! Propagate any updates to parameters and recalculate information as neeeded /*! Should be called at each fit iteration */ void propagateParUpdates() override; //! Retrieve the number of segments in the spline std::size_t nSegments() const { return effSpline_->getnKnots() - 1; } //! Retrieve the positions of the spline knots const std::vector& knotPositions() const { return effSpline_->getXValues(); } //! Fix all knots void fixKnots(); //! Float all knots void floatKnots(); //! Fix or float a specific knot /*! \param knotIndex the index of the knot to fix/float \param fixed true if knot to be fixed, false if knot to be floated */ void fixKnot( const std::size_t knotIndex, const bool fixed ); //! Retrieve the polynomial coefficients for the given spline segment /*! \param segmentIndex the index of the spline segment \return the polynomial coefficients */ - std::array getCoefficients( const std::size_t segmentIndex ) const { return effSpline_->getCoefficients(segmentIndex); } + std::array getCoefficients( const std::size_t segmentIndex ) const; + + //! Retrieve whether any of the parameters are floating in the fit + bool anythingFloating() const { return anyKnotFloating_; } //! Retrieve whether any of the parameters have changed in the latest fit iteration const bool& anythingChanged() const override { return anyKnotChanged_; } private: + //! Perform the initial setup - shared by the various constructors + void initialSetup( const LauSplineDecayTimeEfficiency* other = nullptr ); + //! Update the cached parameter values void updateParameterCache(); + //! The prefix for the parameter names + const TString prefix_; + //! The spline std::unique_ptr effSpline_; + //! The other spline + std::unique_ptr> otherSpline_; + //! The spline parameters std::vector> ownedParams_; //! The spline parameters (for giving access) std::vector params_; //! The spline values std::vector values_; + //! Are any of our knot parameters floating in the fit? + bool ourKnotsFloating_{false}; + + //! Are any of the other spline's knot parameters floating in the fit? + bool otherKnotsFloating_{false}; + //! Are any of the knot parameters floating in the fit? bool anyKnotFloating_{false}; + //! Have any of our knot parameters changed in the latest fit iteration? + bool ourKnotsChanged_{false}; + + //! Have any of the other spline's knot parameters changed in the latest fit iteration? + bool otherKnotsChanged_{false}; + //! Have any of the knot parameters changed in the lastest fit iteration? bool anyKnotChanged_{false}; ClassDefOverride(LauSplineDecayTimeEfficiency, 0) }; +templateClassImp(LauSplineDecayTimeEfficiency); + +template +void LauSplineDecayTimeEfficiency::initialSetup( const LauSplineDecayTimeEfficiency* other ) +{ + if ( not effSpline_ ) { + std::cerr<<"ERROR in LauSplineDecayTimeEfficiency::initialSetup : supplied efficiency spline pointer is null"<Exit(EXIT_FAILURE); + } + + if constexpr ( Order == LauSplineOrder::Sixth ) { + std::vector ourKnots { effSpline_->getXValues() }; + std::vector otherKnots { otherSpline_->knotPositions() }; + + if ( ourKnots != otherKnots ) { + std::cerr<<"ERROR in LauSplineDecayTimeEfficiency::initialSetup : mismatch in knot positions"<Exit(EXIT_FAILURE); + } + + // TODO - any other checks we need to do? + } + + values_ = effSpline_->getYValues(); + + const std::size_t nPars { values_.size() }; + ownedParams_.reserve( nPars ); + + if ( other ) { + for ( auto& ptr : other->ownedParams_ ) { + // TODO - maybe need a way to give these a unique prefix? + ownedParams_.emplace_back( ptr->createClone() ); + } + } else { + for( std::size_t index{0}; index < nPars; ++index ) { + ownedParams_.emplace_back( std::make_unique( Form( "%s_knot_%lu", prefix_.Data(), index ), values_[index], 0.0, 2.0, kTRUE ) ); + } + } + + if constexpr ( Order == LauSplineOrder::Sixth ) { + params_ = otherSpline_->getParameters(); + params_.reserve( 2*nPars ); + } else { + params_.reserve( nPars ); + } + + for ( auto& ptr : ownedParams_ ) { + params_.push_back( ptr.get() ); + } +} + +template +LauSplineDecayTimeEfficiency::LauSplineDecayTimeEfficiency(const LauSplineDecayTimeEfficiency& other) : + prefix_{ other.prefix_ }, + effSpline_{ std::make_unique( * other.effSpline_ ) }, + otherSpline_{ other.otherSpline_ ? std::make_unique>( *other.otherSpline_ ) : nullptr } +{ + this->initialSetup( &other ); +} + +template +void LauSplineDecayTimeEfficiency::fixKnots() +{ + for ( auto& ptr : ownedParams_ ) { + ptr->fixed(true); + } +} + +template +void LauSplineDecayTimeEfficiency::floatKnots() +{ + for ( auto& ptr : ownedParams_ ) { + ptr->fixed(false); + } +} + +template +void LauSplineDecayTimeEfficiency::fixKnot( const std::size_t knotIndex, const bool fixed ) +{ + if ( knotIndex >= ownedParams_.size() ) { + std::cerr<<"ERROR in LauSplineDecayTimeEfficiency::fixKnot : supplied knot index is out of range"<fixed(fixed); +} + +template +Double_t LauSplineDecayTimeEfficiency::getEfficiency( const Double_t abscissa ) const +{ + Double_t eff { effSpline_->evaluate( abscissa ) }; + if constexpr ( Order == LauSplineOrder::Sixth ) { + eff *= otherSpline_->getEfficiency( abscissa ); + } + return eff; +} + +template +std::array::nCoeffs> LauSplineDecayTimeEfficiency::getCoefficients( const std::size_t segmentIndex ) const +{ + if constexpr ( Order == LauSplineOrder::Cubic ) { + return effSpline_->getCoefficients(segmentIndex); + } else { + std::array ourCoeffs { effSpline_->getCoefficients(segmentIndex) }; + std::array otherCoeffs { otherSpline_->getCoefficients(segmentIndex) }; + std::array coeffs; + coeffs[0] = ourCoeffs[0] * otherCoeffs[0]; + coeffs[1] = ourCoeffs[0] * otherCoeffs[1] + ourCoeffs[1] * otherCoeffs[0]; + coeffs[2] = ourCoeffs[0] * otherCoeffs[2] + ourCoeffs[1] * otherCoeffs[1] + ourCoeffs[2] * otherCoeffs[0]; + coeffs[3] = ourCoeffs[0] * otherCoeffs[3] + ourCoeffs[1] * otherCoeffs[2] + ourCoeffs[2] * otherCoeffs[1] + ourCoeffs[3] * otherCoeffs[0]; + coeffs[4] = ourCoeffs[1] * otherCoeffs[3] + ourCoeffs[2] * otherCoeffs[2] + ourCoeffs[3] * otherCoeffs[1]; + coeffs[5] = ourCoeffs[2] * otherCoeffs[3] + ourCoeffs[3] * otherCoeffs[2]; + coeffs[6] = ourCoeffs[3] * otherCoeffs[3]; + return coeffs; + } +} + +template +void LauSplineDecayTimeEfficiency::initialise() +{ + static auto isFixed = []( const std::unique_ptr& par ){ return par->fixed(); }; + + anyKnotFloating_ = ourKnotsFloating_ = not std::all_of( ownedParams_.begin(), ownedParams_.end(), isFixed ); + + if constexpr ( Order == LauSplineOrder::Sixth ) { + otherSpline_->initialise(); + otherKnotsFloating_ = otherSpline_->anythingFloating(); + anyKnotFloating_ |= otherKnotsFloating_; + } + + this->updateParameterCache(); +} + +template +void LauSplineDecayTimeEfficiency::updateParameterCache() +{ + static auto assignValue = []( const std::unique_ptr& par ){ return par->unblindValue(); }; + + // Update our local cache + std::transform( ownedParams_.begin(), ownedParams_.end(), values_.begin(), assignValue ); + + // Update the spline itself + effSpline_->updateYValues( values_ ); +} + +template +void LauSplineDecayTimeEfficiency::propagateParUpdates() +{ + // If none of the parameters are floating there's nothing to do + if ( not anyKnotFloating_ ) { + return; + } + + // Otherwise, determine which of the floating parameters have changed (if any) and act accordingly + + if ( ourKnotsFloating_ ) { + static auto checkEquality = []( const std::unique_ptr& par, const Double_t cacheVal ){ return par->unblindValue() == cacheVal; }; + + anyKnotChanged_ = ourKnotsChanged_ = not std::equal( ownedParams_.begin(), ownedParams_.end(), values_.begin(), checkEquality); + + // Update the acceptance spline if any of the knot values have changed + if ( ourKnotsChanged_ ) { + this->updateParameterCache(); + } + } + + if constexpr ( Order == LauSplineOrder::Sixth ) { + if ( otherKnotsFloating_ ) { + otherSpline_->propagateParUpdates(); + otherKnotsChanged_ = otherSpline_->anythingChanged(); + anyKnotChanged_ |= otherKnotsChanged_; + } + } +} + #endif diff --git a/inc/Laura++_LinkDef.h b/inc/Laura++_LinkDef.h index cd9715f..7b08c57 100644 --- a/inc/Laura++_LinkDef.h +++ b/inc/Laura++_LinkDef.h @@ -1,187 +1,185 @@ /* Copyright 2013 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 */ #ifdef __CINT__ #pragma link off all globals; #pragma link off all classes; #pragma link off all functions; #pragma link C++ class Lau1DCubicSpline+; #pragma link C++ class Lau1DHistPdf+; #pragma link C++ class Lau2DAbsDP+; #pragma link C++ class Lau2DAbsDPPdf+; #pragma link C++ class Lau2DAbsHistDP+; #pragma link C++ class Lau2DAbsHistDPPdf+; #pragma link C++ class Lau2DCubicSpline+; #pragma link C++ class Lau2DHistDP+; #pragma link C++ class Lau2DHistDPPdf+; #pragma link C++ class Lau2DHistPdf+; #pragma link C++ class Lau2DSplineDP+; #pragma link C++ class Lau2DSplineDPPdf+; #pragma link C++ class LauAbsBkgndDPModel+; #pragma link C++ class LauAbsCoeffSet+; #pragma link C++ class LauAbsEffModel+; #pragma link C++ class LauAbsFitter+; #pragma link C++ class LauAbsFitModel+; #pragma link C++ class LauAbsIncohRes+; #pragma link C++ class LauAbsModIndPartWave+; #pragma link C++ class LauAbsPdf+; #pragma link C++ class LauAbsResonance+; #pragma link C++ class LauAbsRValue+; #pragma link C++ class LauArgusPdf+; #pragma link C++ class LauAsymmCalc+; #pragma link C++ class LauBelleCPCoeffSet+; #pragma link C++ class LauBelleNR+; #pragma link C++ class LauBelleSymNR+; #pragma link C++ class LauBifurcatedGaussPdf+; #pragma link C++ class LauBkgndDPModel+; #pragma link C++ class LauBlattWeisskopfFactor+; #pragma link C++ class LauBlind+; #pragma link C++ class LauBreitWignerRes+; #pragma link C++ class LauBsCartesianCPCoeffSet+; #pragma link C++ class LauBsCPFitModel+; #pragma link C++ class LauCacheData+; #pragma link C++ class LauCalcChiSq+; #pragma link C++ class LauCartesianCPCoeffSet+; #pragma link C++ class LauCartesianGammaCPCoeffSet+; #pragma link C++ class LauCategoryFlavTag+; #pragma link C++ class LauChebychevPdf+; #pragma link C++ class LauCleoCPCoeffSet+; #pragma link C++ class LauComplex+; #pragma link C++ class LauCPFitModel+; #pragma link C++ class LauCruijffPdf+; #pragma link C++ class LauCrystalBallPdf+; #pragma link C++ class LauDabbaRes+; #pragma link C++ class LauDatabasePDG+; #pragma link C++ class LauDaughters+; #pragma link C++ class LauDecayTimePdf+; #pragma link C++ class LauDPDepBifurGaussPdf+; #pragma link C++ class LauDPDepCruijffPdf+; #pragma link C++ class LauDPDepGaussPdf+; #pragma link C++ class LauDPDepMapPdf+; #pragma link C++ class LauDPDepSumPdf+; #pragma link C++ class LauEffModel+; #pragma link C++ class LauEFKLLMRes+; #pragma link C++ class LauEmbeddedData+; #pragma link C++ class LauExponentialPdf+; #pragma link C++ class LauFitDataTree+; #pragma link C++ class LauFitNtuple+; #pragma link C++ class LauFitter+; #pragma link C++ class LauFitObject+; #pragma link C++ class LauFlatteRes+; #pragma link C++ class LauFlatNR+; #pragma link C++ class LauFlavTag+; #pragma link C++ class LauFormulaPar+; #pragma link C++ class LauGaussIncohRes+; #pragma link C++ class LauGaussPdf+; #pragma link C++ class LauGenNtuple+; #pragma link C++ class LauGounarisSakuraiRes+; #pragma link C++ class LauIntegrals+; #pragma link C++ class LauDPPartialIntegralInfo+; #pragma link C++ class LauIsobarDynamics+; #pragma link C++ class LauKappaRes+; #pragma link C++ class LauKinematics+; #pragma link C++ class LauKMatrixProdPole+; #pragma link C++ class LauKMatrixProdSVP+; #pragma link C++ class LauKMatrixPropagator+; #pragma link C++ class LauKMatrixPropFactory+; #pragma link C++ class LauLASSBWRes+; #pragma link C++ class LauLASSNRRes+; #pragma link C++ class LauLASSRes+; #pragma link C++ class LauLinearPdf+; #pragma link C++ class LauMagPhaseCoeffSet+; #pragma link C++ class LauMagPhaseCPCoeffSet+; #pragma link C++ class LauMergeDataFiles+; #pragma link C++ class LauMinuit+; #pragma link C++ class LauModIndPartWaveMagPhase+; #pragma link C++ class LauModIndPartWaveRealImag+; #pragma link C++ class LauNovosibirskPdf+; #pragma link C++ class LauNRAmplitude+; #pragma link C++ class LauParameter+; #pragma link C++ class LauParametricStepFuncPdf+; #pragma link C++ class LauParamFixed+; #pragma link C++ class LauParticlePDG+; #pragma link C++ class LauPolNR+; #pragma link C++ class LauPoleRes+; #pragma link C++ class LauPolarFormFactorNR+; #pragma link C++ class LauPolarFormFactorSymNR+; #pragma link C++ class LauPolarGammaCPCoeffSet+; #pragma link C++ class LauPrint+; #pragma link C++ class LauRealImagCoeffSet+; #pragma link C++ class LauRealImagCPCoeffSet+; #pragma link C++ class LauRealImagGammaCPCoeffSet+; #pragma link C++ class LauRelBreitWignerRes+; #pragma link C++ class LauResonanceInfo+; #pragma link C++ class LauRescatteringRes+; #pragma link C++ class LauRescattering2Res+; #pragma link C++ class LauResonanceMaker+; #pragma link C++ class LauResultsExtractor+; #pragma link C++ class LauRhoOmegaMix+; #ifdef DOLAUROOFITTASK #pragma link C++ class LauRooFitTask+; #endif #pragma link C++ class LauScfMap+; #pragma link C++ class LauSigmaRes+; #pragma link C++ class LauSigmoidPdf+; #pragma link C++ class LauSimpleFitModel+; #pragma link C++ class LauSimFitCoordinator+; #pragma link C++ class LauSimFitTask+; #pragma link C++ class LauSPlot+; #pragma link C++ class LauString+; #pragma link C++ class LauSumPdf+; #pragma link C++ class LauTextFileParser+; #pragma link C++ class LauTimeDepFlavModel+; #pragma link C++ class LauTimeDepFitModel+; #pragma link C++ class LauTimeDepNonFlavModel+; #pragma link C++ class LauVetoes+; #pragma link C++ class LauWeightedSumEffModel+; #pragma link C++ class LauAbsDecayTimeCalculator+; #pragma link C++ class LauAbsDecayTimeEfficiency+; #pragma link C++ class LauAbsDecayTimeIntegrator+; -#pragma link C++ class LauAbsSplineDecayTimeEfficiency+; -//#pragma link C++ class LauAbsSplineDecayTimeEfficiency+; #pragma link C++ class LauBinnedDecayTimeEfficiency+; #pragma link C++ class LauDecayTimePhysicsModel+; #pragma link C++ class LauDecayTimeResolution+; #pragma link C++ class LauNonSmearedDecayTimeCalculator+; #pragma link C++ class LauNonSmearedBinnedEfficiencyDecayTimeIntegrator+; #pragma link C++ class LauNonSmearedSplineEfficiencyDecayTimeIntegrator+; -//#pragma link C++ class LauNonSmearedSplineEfficiencyDecayTimeIntegrator+; +#pragma link C++ class LauNonSmearedSplineEfficiencyDecayTimeIntegrator+; #pragma link C++ class LauNonSmearedUniformEfficiencyDecayTimeIntegrator+; #pragma link C++ class LauSmearedDecayTimeCalculator+; #pragma link C++ class LauSmearedBinnedEfficiencyDecayTimeIntegrator+; #pragma link C++ class LauSmearedSplineEfficiencyDecayTimeIntegrator+; -//#pragma link C++ class LauSmearedSplineEfficiencyDecayTimeIntegrator+; +#pragma link C++ class LauSmearedSplineEfficiencyDecayTimeIntegrator+; #pragma link C++ class LauSmearedUniformEfficiencyDecayTimeIntegrator+; #pragma link C++ class LauSplineDecayTimeEfficiency+; -//#pragma link C++ class LauSplineDecayTimeEfficiency+; +#pragma link C++ class LauSplineDecayTimeEfficiency+; #pragma link C++ class LauUniformDecayTimeEfficiency+; #pragma link C++ namespace LauConstants+; #pragma link C++ namespace LauDecayTime+; #pragma link C++ namespace LauRandom+; #endif diff --git a/src/Lau1DCubicSpline.cc b/src/Lau1DCubicSpline.cc index beb65d6..272ef7c 100644 --- a/src/Lau1DCubicSpline.cc +++ b/src/Lau1DCubicSpline.cc @@ -1,536 +1,532 @@ /* 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.cc \brief File containing implementation of Lau1DCubicSpline class. */ #include #include #include #include #include #include "Lau1DCubicSpline.hh" ClassImp(Lau1DCubicSpline) Lau1DCubicSpline::Lau1DCubicSpline(const std::vector& xs, const std::vector& ys, LauSplineType type, LauSplineBoundaryType leftBound, LauSplineBoundaryType rightBound, Double_t dydx0, Double_t dydxn) : nKnots_(xs.size()), x_(xs), y_(ys), type_(type), leftBound_(leftBound), rightBound_(rightBound), dydx0_(dydx0), dydxn_(dydxn) { init(); } -Lau1DCubicSpline::~Lau1DCubicSpline() -{ -} - Double_t Lau1DCubicSpline::evaluate(Double_t x) const { // do not attempt to extrapolate the spline if( xx_[nKnots_-1] ) { std::cerr << "WARNING in Lau1DCubicSpline::evaluate : function is only defined between " << x_[0] << " and " << x_[nKnots_-1] << std::endl; std::cerr << " value at " << x << " returned as 0" << std::endl; return 0.; } // first determine which 'cell' of the spline x is in // cell i runs from knot i to knot i+1 Int_t cell(0); while( x > x_[cell+1] ) { ++cell; } // obtain x- and y-values of the neighbouring knots Double_t xLow = x_[cell]; Double_t xHigh = x_[cell+1]; Double_t yLow = y_[cell]; Double_t yHigh = y_[cell+1]; if(type_ == Lau1DCubicSpline::LinearInterpolation) { return yHigh*(x-xLow)/(xHigh-xLow) + yLow*(xHigh-x)/(xHigh-xLow); } // obtain t, the normalised x-coordinate within the cell, // and the coefficients a and b, which are defined in cell i as: // // a_i = k_i *(x_i+1 - x_i) - (y_i+1 - y_i), // b_i = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i) // // where k_i is (by construction) the first derivative at knot i Double_t t = (x - xLow) / (xHigh - xLow); Double_t a = dydx_[cell] * (xHigh - xLow) - (yHigh - yLow); Double_t b = -1.*dydx_[cell+1] * (xHigh - xLow) + (yHigh - yLow); Double_t retVal = (1 - t) * yLow + t * yHigh + t * (1 - t) * ( a * (1 - t) + b * t ); return retVal; } void Lau1DCubicSpline::updateYValues(const std::vector& ys) { y_ = ys; this->calcDerivatives(); } void Lau1DCubicSpline::updateYValues(const std::vector& ys) { for (UInt_t i=0; iunblindValue(); } this->calcDerivatives(); } void Lau1DCubicSpline::updateYValues(const std::vector& ys) { for (UInt_t i=0; iunblindValue(); } this->calcDerivatives(); } void Lau1DCubicSpline::updateType(LauSplineType type) { if(type_ != type) { type_ = type; this->calcDerivatives(); } } void Lau1DCubicSpline::updateBoundaryConditions(LauSplineBoundaryType leftBound, LauSplineBoundaryType rightBound, Double_t dydx0, Double_t dydxn) { Bool_t updateDerivatives(kFALSE); if(leftBound_ != leftBound || rightBound_ != rightBound) { leftBound_ = leftBound; rightBound_ = rightBound; updateDerivatives = kTRUE; } if(dydx0_ != dydx0) { dydx0_ = dydx0; if(leftBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = kTRUE; } if(dydxn_ != dydxn) { dydxn_ = dydxn; if(rightBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = kTRUE; } if(updateDerivatives) { this->calcDerivatives(); } } std::array Lau1DCubicSpline::getCoefficients(const UInt_t i, const bool normalise) const { std::array result = {0.,0.,0.,0.}; if(i >= nKnots_-1) { std::cerr << "ERROR in Lau1DCubicSpline::getCoefficients requested for too high a knot value" << std::endl; return result; } Double_t xL = x_[i] , xH = x_[i+1]; Double_t yL = y_[i] , yH = y_[i+1]; Double_t h = xH-xL; //This number comes up a lot switch(type_) { case Lau1DCubicSpline::StandardSpline: case Lau1DCubicSpline::AkimaSpline: { Double_t kL = dydx_[i], kH = dydx_[i+1]; //a and b based on definitions from https://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline Double_t a = kL*h-(yH-yL); Double_t b =-kH*h+(yH-yL); Double_t denom = -h*h*h;//The terms have a common demoninator result[0] = -b*xL*xL*xH + a*xL*xH*xH + h*h*(xL*yH - xH*yL); result[1] = -a*xH*(2*xL+xH) + b*xL*(xL+2*xH) + h*h*(yL-yH); result[2] = -b*(2*xL+xH) + a*(xL+2*xH); result[3] = -a+b; for(auto& res : result){res /= denom;} break; } /* case Lau1DCubicSpline::AkimaSpline: // Double check the Akima description of splines (in evaluate) right now they're the same except for the first derivatives { //using fomulae from https://asmquantmacro.com/2015/09/01/akima-spline-interpolation-in-excel/ std::function m = [&](Int_t j) //formula to get the straight line gradient { if(j < 0){return 2*m(j+1)-m(j+2);} if(j >= nKnots_){return 2*m(j-1)-m(j-2);} return (y_[j+1]-y_[j]) / (x_[j+1]-x_[j]); }; auto t = [&](Int_t j) { Double_t res = 0.; //originally res was called 't' but that produced a shadow warning Double_t denom = TMath::Abs(m(j+1)-m(j)) + TMath::Abs(m(j-1)-m(j-2)); if(denom == 0){res = (m(j)-m(j-1))/2;} //Special case for when denom = 0 else { res = TMath::Abs(m(j+1)-m(j))*m(j-1) + TMath::Abs(m(j-1)-m(j-2))*m(j); res /= denom; } return res; }; //These are the p's to get the spline in the form p_k(x-xL)^k Double_t pDenom = x_[i+1]/x_[i]; //a denominator used for p[2] and p[3] std::array p = {y_[i],t(i),0.,0.}; //we'll do the last 2 below p[2] = 3*m(i)-2*t(i)-t(i+1); p[2]/= pDenom; p[3] = t(i)+t(i+1)-2*m(i); p[3]/= pDenom*pDenom; //Now finally rearranging the p's into the desired results result[0] = p[0]-p[1]*xL+p[2]*xL*xL-p[3]*xL*xL*xL; result[1] = p[1]-2*p[2]*xL+3*p[3]*xL*xL; result[2] = p[2]-3*p[3]*xL; result[3] = p[3]; break; }*/ case Lau1DCubicSpline::LinearInterpolation: { result[0] = xH*yL-xL*yH; result[1] = yH-yL; for(auto& res : result){res /= h;} break; } } if(normalise) { Double_t integral = this->integral(); for(auto& res : result){res /= integral;} } return result; } Double_t Lau1DCubicSpline::integral() const { Double_t integral = 0.; for(UInt_t iKnot = 0; iKnot < nKnots_ -1; ++iKnot) { Double_t minAbs = x_[iKnot]; Double_t maxAbs = x_[iKnot+1]; std::array coeffs = this -> getCoefficients(iKnot, false); auto integralFunc = [&coeffs](Double_t x){return coeffs[0]*x + coeffs[1]*x*x/2 + coeffs[2]*x*x*x/3 + coeffs[3]*x*x*x*x/4;}; integral += integralFunc(maxAbs); integral -= integralFunc(minAbs); } return integral; } TF1* Lau1DCubicSpline::makeTF1(const bool normalise) const { TString functionString = ""; //make a long piecewise construction of all the spline pieces for(UInt_t i = 0; i < nKnots_-1; ++i) { functionString += Form("(x>%f && x<= %f)*",x_[i],x_[i+1]);//get the bounds of this piece std::array coeffs = this->getCoefficients(i,normalise); functionString += Form("(%f + %f*x + %f*x^2 + %f*x^3)",coeffs[0],coeffs[1],coeffs[2],coeffs[3]); if(i < nKnots_ -2){functionString += " + \n";}//add to all lines except the last } TF1* func = new TF1("SplineFunction", functionString, x_.front(), x_.back()); return func; } TF1* Lau1DCubicSpline::normaliseToTH1(TH1* hist) const { //first define the fit function auto fitf = [this](Double_t* x, Double_t* par) {//there is only 1 x (the abscissa) and 1 par (a scaling of the entire thing) return this->evaluate( x[0] ) * par[0]; }; //Make the function TF1* FittedFunc = new TF1("FittedFunc",fitf,x_.front(),x_.back(),1); //Set the parameter name FittedFunc -> SetParNames("Constant"); //Set the parameter limits and default value FittedFunc -> SetParLimits(0,0.,10.); FittedFunc -> SetParameter(0,1.); hist->Fit("FittedFunc","N O V"); return FittedFunc; } void Lau1DCubicSpline::fitToTH1(TH1* hist) { auto fitf = [this](Double_t* x, Double_t* par) { this -> updateYValues( std::vector(par, par + nKnots_) ); return this -> evaluate( x[0] ); }; //Make the function TF1 FittedFunc("FittedFunc",fitf,x_.front(),x_.back(),nKnots_); const Double_t knotMax = hist->GetMaximum() * 1.5; for(UInt_t knot = 0; knot <= nKnots_ ; ++knot) { FittedFunc.SetParName(knot, Form("Knot%u",knot)); FittedFunc.SetParLimits(knot, 0., knotMax); FittedFunc.SetParameter(knot, y_[knot]); } hist->Fit("FittedFunc","N O V"); return; } void Lau1DCubicSpline::init() { if( y_.size() != x_.size()) { std::cerr << "ERROR in Lau1DCubicSpline::init : The number of y-values given does not match the number of x-values" << std::endl; std::cerr << " Found " << y_.size() << ", expected " << x_.size() << std::endl; gSystem->Exit(EXIT_FAILURE); } if( nKnots_ < 3 ) { std::cerr << "ERROR in Lau1DCubicSpline::init : The number of knots is too small" << std::endl; std::cerr << " Found " << nKnots_ << ", expected at least 3 (to have at least 1 internal knot)" << std::endl; gSystem->Exit(EXIT_FAILURE); } dydx_.assign(nKnots_,0.0); a_.assign(nKnots_,0.0); b_.assign(nKnots_,0.0); c_.assign(nKnots_,0.0); d_.assign(nKnots_,0.0); this->calcDerivatives(); } void Lau1DCubicSpline::calcDerivatives() { switch ( type_ ) { case Lau1DCubicSpline::StandardSpline : this->calcDerivativesStandard(); break; case Lau1DCubicSpline::AkimaSpline : this->calcDerivativesAkima(); break; case Lau1DCubicSpline::LinearInterpolation : //derivatives not needed for linear interpolation break; } } void Lau1DCubicSpline::calcDerivativesStandard() { // derivatives are determined such that the second derivative is continuous at internal knots // derivatives, k_i, are the solutions to a set of linear equations of the form: // a_i * k_i-1 + b_i * k+i + c_i * k_i+1 = d_i with a_0 = 0, c_n-1 = 0 // this is solved using the tridiagonal matrix algorithm as on en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm // first and last equations give boundary conditions // - for natural boundary, require f''(x) = 0 at end knot // - for 'not a knot' boundary, require f'''(x) continuous at second knot // - for clamped boundary, require predefined value of f'(x) at end knot // non-zero values of a_0 and c_n-1 would give cyclic boundary conditions a_[0] = 0.; c_[nKnots_-1] = 0.; // set left boundary condition if(leftBound_ == Lau1DCubicSpline::Natural) { b_[0] = 2./(x_[1]-x_[0]); c_[0] = 1./(x_[1]-x_[0]); d_[0] = 3.*(y_[1]-y_[0])/((x_[1]-x_[0])*(x_[1]-x_[0])); } else if(leftBound_ == Lau1DCubicSpline::NotAKnot) { // define the width, h, and the 'slope', delta, of the first cell Double_t h1(x_[1]-x_[0]), h2(x_[2]-x_[0]); Double_t delta1((y_[1]-y_[0])/h1), delta2((y_[2]-y_[1])/h2); // these coefficients can be determined by requiring f'''_0(x_1) = f'''_1(x_1) // the requirement f''_0(x_1) = f''_1(x_1) has been used to remove the dependence on k_2 b_[0] = h2; c_[0] = h1+h2; d_[0] = delta1*(2.*h2*h2 + 3.*h1*h2)/(h1+h2) + delta2*5.*h1*h1/(h1+h2); } else { //Clamped b_[0] = 1.; c_[0] = 0.; d_[0] = dydx0_; } // set right boundary condition if(rightBound_ == Lau1DCubicSpline::Natural) { a_[nKnots_-1] = 1./(x_[nKnots_-1]-x_[nKnots_-2]); b_[nKnots_-1] = 2./(x_[nKnots_-1]-x_[nKnots_-2]); d_[nKnots_-1] = 3.*(y_[nKnots_-1]-y_[nKnots_-2])/((x_[nKnots_-1]-x_[nKnots_-2])*(x_[nKnots_-1]-x_[nKnots_-2])); } else if(rightBound_ == Lau1DCubicSpline::NotAKnot) { // define the width, h, and the 'slope', delta, of the last cell Double_t hnm1(x_[nKnots_-1]-x_[nKnots_-2]), hnm2(x_[nKnots_-2]-x_[nKnots_-3]); Double_t deltanm1((y_[nKnots_-1]-y_[nKnots_-2])/hnm1), deltanm2((y_[nKnots_-2]-y_[nKnots_-3])/hnm2); // these coefficients can be determined by requiring f'''_n-3(x_n-2) = f'''_n-2(x_n-2) // the requirement f''_n-3(x_n-2) = f''_n-2(x_n-2) has been used to remove // the dependence on k_n-3 a_[nKnots_-1] = hnm2 + hnm1; b_[nKnots_-1] = hnm1; d_[nKnots_-1] = deltanm2*hnm1*hnm1/(hnm2+hnm1) + deltanm1*(2.*hnm2*hnm2 + 3.*hnm2*hnm1)/(hnm2+hnm1); } else { //Clamped a_[nKnots_-1] = 0.; b_[nKnots_-1] = 1.; d_[nKnots_-1] = dydxn_; } // the remaining equations ensure that f_i-1''(x_i) = f''_i(x_i) for all internal knots for(UInt_t i=1; i=0; --i) { dydx_[i] = d_[i] - c_[i]*dydx_[i+1]; } } void Lau1DCubicSpline::calcDerivativesAkima() { //derivatives are calculated according to the Akima method // J.ACM vol. 17 no. 4 pp 589-602 Double_t am1(0.), an(0.), anp1(0.); // a[i] is the slope of the segment from point i-1 to point i // // n.b. segment 0 is before point 0 and segment n is after point n-1 // internal segments are numbered 1 - n-1 for(UInt_t i=1; i -#include - -#include "Rtypes.h" -#include "TSystem.h" - -#include "Lau1DCubicSpline.hh" -#include "LauSplineDecayTimeEfficiency.hh" -#include "LauParameter.hh" -#include "LauParamFixed.hh" - -ClassImp(LauSplineDecayTimeEfficiency); - -LauSplineDecayTimeEfficiency::LauSplineDecayTimeEfficiency( std::unique_ptr effSpline ) : effSpline_{std::move(effSpline)} -{ - if ( not effSpline_ ) { - std::cerr<<"ERROR in LauSplineDecayTimeEfficiency::LauSplineDecayTimeEfficiency : supplied efficiency spline pointer is null"<Exit(EXIT_FAILURE); - } - - values_ = effSpline_->getYValues(); - - const std::size_t nPars { values_.size() }; - ownedParams_.reserve( nPars ); - - for( std::size_t index{0}; index < nPars; ++index ) { - ownedParams_.emplace_back( std::make_unique( Form( "dteff_knot_%lu", index ), values_[index], 0.0, 1.0, kTRUE ) ); - } - - params_.reserve( nPars ); - for ( auto& ptr : ownedParams_ ) { - params_.push_back( ptr.get() ); - } -} - -void LauSplineDecayTimeEfficiency::fixKnots() -{ - for ( auto& ptr : ownedParams_ ) { - ptr->fixed(true); - } -} - -void LauSplineDecayTimeEfficiency::floatKnots() -{ - for ( auto& ptr : ownedParams_ ) { - ptr->fixed(false); - } -} - -void LauSplineDecayTimeEfficiency::fixKnot( const std::size_t knotIndex, const bool fixed ) -{ - if ( knotIndex >= ownedParams_.size() ) { - std::cerr<<"ERROR in LauSplineDecayTimeEfficiency::fixKnot : supplied knot index is out of range"<fixed(fixed); -} - -Double_t LauSplineDecayTimeEfficiency::getEfficiency( const Double_t abscissa ) const -{ - return effSpline_->evaluate( abscissa ); -} - -void LauSplineDecayTimeEfficiency::initialise() -{ - LauParamFixed isFixed; - anyKnotFloating_ = not std::all_of( params_.begin(), params_.end(), isFixed ); - - this->updateParameterCache(); -} - -void LauSplineDecayTimeEfficiency::updateParameterCache() -{ - // Update the spline itself - effSpline_->updateYValues( params_ ); - - static auto assignValue = [](const LauAbsRValue* par){return par->unblindValue();}; - - // Update our local cache - std::transform( params_.begin(), params_.end(), values_.begin(), assignValue ); -} - -void LauSplineDecayTimeEfficiency::propagateParUpdates() -{ - // If none of the parameters are floating there's nothing to do - if ( not anyKnotFloating_ ) { - return; - } - - // Otherwise, determine which of the floating parameters have changed (if any) and act accordingly - - // TODO - these lambda's are repeated in lots of places, should we put them in some utility include? - static auto checkEquality = [](const LauAbsRValue* par, const Double_t cacheVal){return par->unblindValue() == cacheVal;}; - - anyKnotChanged_ = not std::equal(params_.begin(), params_.end(), values_.begin(), checkEquality); - - // Update the acceptance spline if any of the knot values have changed - if ( anyKnotChanged_ ) { - this->updateParameterCache(); - } -}