diff --git a/inc/LauNonSmearedBinnedEfficiencyDecayTimeIntegrator.hh b/inc/LauNonSmearedBinnedEfficiencyDecayTimeIntegrator.hh index cdbfd08..0e904e1 100644 --- a/inc/LauNonSmearedBinnedEfficiencyDecayTimeIntegrator.hh +++ b/inc/LauNonSmearedBinnedEfficiencyDecayTimeIntegrator.hh @@ -1,85 +1,85 @@ /* 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 LauNonSmearedBinnedEfficiencyDecayTimeIntegrator.hh \brief File containing declaration of LauNonSmearedBinnedEfficiencyDecayTimeIntegrator class. */ #ifndef LAU_NONSMEARED_BINNEDEFFICIENCY_DECAYTIME_INTEGRATOR #define LAU_NONSMEARED_BINNEDEFFICIENCY_DECAYTIME_INTEGRATOR #include #include "Rtypes.h" #include "LauAbsDecayTimeIntegrator.hh" #include "LauBinnedDecayTimeEfficiency.hh" #include "LauDecayTimePhysicsModel.hh" #include "LauNonSmearedUniformEfficiencyDecayTimeIntegrator.hh" #include "LauUniformDecayTimeEfficiency.hh" class LauParameter; /*! \class LauNonSmearedBinnedEfficiencyDecayTimeIntegrator - \brief Class for defining the decay time integrator for the case of no resolution and uniform efficiency + \brief Class for defining the decay time integrator for the case of no resolution and binned efficiency */ class LauNonSmearedBinnedEfficiencyDecayTimeIntegrator : public LauAbsDecayTimeIntegrator { public: //! Constructor LauNonSmearedBinnedEfficiencyDecayTimeIntegrator( const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const LauDecayTimePhysicsModel& physModel, const LauBinnedDecayTimeEfficiency& effModel ); //! Cache information from data /*! \param [in] abscissaErrors the values of the per-event decay time error */ void cacheInfo( const std::vector& abscissaErrors ) override; //! Propagate any updates to parameters and recalculate information as neeeded void propagateParUpdates() override; //! Retrieve the normalisation terms (optionally for a given event) /*! \param [in] iEvt the event index */ LauDecayTimeNormTerms getNormTerms( const std::size_t iEvt ) const override; //! Calculate the normalisation terms for a given value of the decay time error /*! \param [in] abscissaError the decay time error value (actually {1.0,err} to allow branchless control of scaling) */ LauDecayTimeNormTerms calcNormTerms( const std::array& abscissaError = {1.0,1.0} ) const override; private: //! The collection of uniform efficiency objects which we pass to the integrators std::vector efficiencies_; //! The collection of uniform efficiency integrator objects to which we delegate std::vector integrators_; ClassDefOverride(LauNonSmearedBinnedEfficiencyDecayTimeIntegrator, 0) }; #endif diff --git a/inc/LauNonSmearedSplineEfficiencyDecayTimeIntegrator.hh b/inc/LauNonSmearedSplineEfficiencyDecayTimeIntegrator.hh index 570aacc..03ba648 100644 --- a/inc/LauNonSmearedSplineEfficiencyDecayTimeIntegrator.hh +++ b/inc/LauNonSmearedSplineEfficiencyDecayTimeIntegrator.hh @@ -1,392 +1,394 @@ /* 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 LauNonSmearedSplineEfficiencyDecayTimeIntegrator.hh \brief File containing declaration of LauNonSmearedSplineEfficiencyDecayTimeIntegrator class. */ #ifndef LAU_NONSMEARED_SPLINEEFFICIENCY_DECAYTIME_INTEGRATOR #define LAU_NONSMEARED_SPLINEEFFICIENCY_DECAYTIME_INTEGRATOR #include #include #include #include "Rtypes.h" #include "TSystem.h" #include "LauAbsDecayTimeIntegrator.hh" #include "LauDecayTime.hh" #include "LauDecayTimePhysicsModel.hh" #include "LauSplineDecayTimeEfficiency.hh" class LauParameter; /*! \class LauNonSmearedSplineEfficiencyDecayTimeIntegrator - \brief Class for defining the decay time integrator for the case of no resolution and uniform efficiency + \brief Class for defining the decay time integrator for the case of no resolution and spline efficiency + + The implementation uses the formalism from arXiv:1407.0748 */ template class LauNonSmearedSplineEfficiencyDecayTimeIntegrator : public LauAbsDecayTimeIntegrator { public: //! Constructor LauNonSmearedSplineEfficiencyDecayTimeIntegrator( const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const LauDecayTimePhysicsModel& physModel, const LauSplineDecayTimeEfficiency& effModel ); //! Cache information from data /*! \param [in] abscissaErrors the values of the per-event decay time error */ void cacheInfo( const std::vector& abscissaErrors ) override; //! Propagate any updates to parameters and recalculate information as neeeded void propagateParUpdates() override; //! Retrieve the normalisation terms (optionally for a given event) /*! \param [in] iEvt the event index */ LauDecayTimeNormTerms getNormTerms( [[maybe_unused]] const std::size_t iEvt ) const override { return normTerms_; } //! Calculate the normalisation terms for a given value of the decay time error /*! \param [in] abscissaError the decay time error value (actually {1.0,err} to allow branchless control of scaling) */ LauDecayTimeNormTerms calcNormTerms( const std::array& abscissaError = {1.0,1.0} ) const override; private: //! The number of coefficients of each spline segment static constexpr std::size_t nCoeffs { static_cast>(Order) + 1 }; //! Define the RealArray typedef using RealArray = std::array; //! Define the Real2DArray typedef using Real2DArray = std::array; //! Define the ComplexArray typedef using ComplexArray = std::array,nCoeffs>; //! Update the cache of the spline coefficients void updateCoeffCache(); //! Update the cache of the integrals for each power of the polynomial /*! \param [in] forceUpdate used to force an update the first time around */ void updateIkCache( const bool forceUpdate = false ); //! Sum over the polynomial orders of the coefficients multiplied by the supplied integrals /*! \param [in] cacheIk integrals for each polynomial degree */ std::complex sumIkTerms( const std::vector& cacheIk ) const; //! Calculate and fill the supplied integral cache given the value of u /*! \param [out] cacheIk cache of integrals for each polynomial degree \param [in] u = Gamma - i * DeltaM in general */ void calcIntegrals( std::vector& cacheIk, const std::complex& u ); //! Calculate a particular term in the integrals /*! \param [in] k the power of the monomial \param [in] minAbs the lower end of the integration range \param [in] maxAbs the upper end of the integration range \param [in] u = Gamma - i * DeltaM in general */ std::complex calcIk( const std::size_t k, const Double_t minAbs, const Double_t maxAbs, const std::complex& u ) const; //! The minimum value of the decay time const Double_t minAbscissa_; //! The maximum value of the decay time const Double_t maxAbscissa_; //! The physics model const LauDecayTimePhysicsModel& physModel_; //! The efficiency model const LauSplineDecayTimeEfficiency& effModel_; //! Cached value of lifetime const Double_t& tauVal_; //! Cached value of 1/lifetime const Double_t& gammaVal_; //! Cached value of mass difference const Double_t& deltaMVal_; //! Cached value of width difference const Double_t& deltaGammaVal_; //! Have any of the physics parameters changed in the latest fit iteration const bool& physicsParamChanged_; //! Has the lifetime changed in the latest fit iteration const bool& tauChanged_; //! Has the mass difference changed in the latest fit iteration const bool& deltaMChanged_; //! Has the width difference changed in the latest fit iteration const bool& deltaGammaChanged_; //! Have any of the efficiency parameters changed in the latest fit iteration const bool& effParamChanged_; //! Factorial values RealArray factorial_; //! Binomial coefficient values Real2DArray binomial_; //! Values of the spline coefficients for each spline segment std::vector coeffs_; //! Caches of the integrals for the exp term std::vector expTermIkVals_; //! Caches of the integrals for the cos and sin terms std::vector trigTermIkVals_; //! Caches of the integrals for the B_H part of the cosh and sinh terms std::vector hypHTermIkVals_; //! Caches of the integrals for the B_L part of the cosh and sinh terms std::vector hypLTermIkVals_; //! The final integrals LauDecayTimeNormTerms normTerms_; ClassDefOverride(LauNonSmearedSplineEfficiencyDecayTimeIntegrator, 0) }; templateClassImp(LauNonSmearedSplineEfficiencyDecayTimeIntegrator); template LauNonSmearedSplineEfficiencyDecayTimeIntegrator::LauNonSmearedSplineEfficiencyDecayTimeIntegrator( const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const LauDecayTimePhysicsModel& physModel, const LauSplineDecayTimeEfficiency& effModel ) : minAbscissa_{minAbscissaVal}, maxAbscissa_{maxAbscissaVal}, physModel_{physModel}, effModel_{effModel}, tauVal_{physModel_.getLifetimeValue()}, gammaVal_{physModel_.getGammaValue()}, deltaMVal_{physModel_.getDeltaMValue()}, deltaGammaVal_{physModel_.getDeltaGammaValue()}, physicsParamChanged_{physModel_.anythingChanged()}, tauChanged_{physModel_.lifetimeChanged()}, deltaMChanged_{physModel_.deltaMChanged()}, deltaGammaChanged_{physModel_.deltaGammaChanged()}, effParamChanged_{effModel_.anythingChanged()} { switch ( physModel_.getFunctionType() ) { case LauDecayTime::FuncType::Hist : case LauDecayTime::FuncType::Delta : case LauDecayTime::FuncType::DeltaExp : // These do not make sense for us: // - Hist has no need of an integrator // - Delta and DeltaExp require a resolution function std::cerr << "ERROR in LauNonSmearedSplineEfficiencyDecayTimeIntegrator::LauNonSmearedSplineEfficiencyDecayTimeIntegrator : Unsupported function type in the physics model" << std::endl; gSystem->Exit(EXIT_FAILURE); break; case LauDecayTime::FuncType::Exp : case LauDecayTime::FuncType::ExpTrig : case LauDecayTime::FuncType::ExpHypTrig : // All fine, we can deal with these break; } // TODO - we should check that the range of the spline matches the range we've been given // populate the binomial and factorial arrays for ( std::size_t i{0}; i < nCoeffs; ++i ) { factorial_[i] = TMath::Factorial(i); for ( std::size_t j{0}; j < nCoeffs; ++j ) { if ( j <= i ) { binomial_[i][j] = TMath::Binomial(i,j); } else { binomial_[i][j] = 0.0; } } } } template void LauNonSmearedSplineEfficiencyDecayTimeIntegrator::cacheInfo( [[maybe_unused]] const std::vector& abscissaErrors ) { const std::size_t nSplineSegments { effModel_.nSegments() }; coeffs_.clear(); coeffs_.resize(nSplineSegments); expTermIkVals_.clear(); expTermIkVals_.resize(nSplineSegments); trigTermIkVals_.clear(); trigTermIkVals_.resize(nSplineSegments); hypHTermIkVals_.clear(); hypHTermIkVals_.resize(nSplineSegments); hypLTermIkVals_.clear(); hypLTermIkVals_.resize(nSplineSegments); this->updateCoeffCache(); this->updateIkCache(true); normTerms_ = this->calcNormTerms(); } template void LauNonSmearedSplineEfficiencyDecayTimeIntegrator::updateCoeffCache() { const std::size_t nSplineSegments { effModel_.nSegments() }; for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { coeffs_[iSeg] = effModel_.getCoefficients(iSeg); } } template void LauNonSmearedSplineEfficiencyDecayTimeIntegrator::updateIkCache( const bool forceUpdate ) { if ( tauChanged_ or forceUpdate ) { const std::complex uExp { gammaVal_ }; this->calcIntegrals( expTermIkVals_, uExp ); } const LauDecayTime::FuncType funcType { physModel_.getFunctionType() }; if ( funcType == LauDecayTime::FuncType::ExpTrig or funcType == LauDecayTime::FuncType::ExpHypTrig ) { if ( tauChanged_ or deltaMChanged_ or forceUpdate ) { const std::complex uTrig { gammaVal_, -deltaMVal_ }; this->calcIntegrals( trigTermIkVals_, uTrig ); } if ( funcType == LauDecayTime::FuncType::ExpHypTrig ) { if ( tauChanged_ or deltaGammaChanged_ or forceUpdate ) { const std::complex uH { gammaVal_ - 0.5 * deltaGammaVal_ }; const std::complex uL { gammaVal_ + 0.5 * deltaGammaVal_ }; this->calcIntegrals( hypHTermIkVals_, uH ); this->calcIntegrals( hypLTermIkVals_, uL ); } } } } template LauDecayTimeNormTerms LauNonSmearedSplineEfficiencyDecayTimeIntegrator::calcNormTerms( [[maybe_unused]] const std::array& abscissaError ) const { LauDecayTimeNormTerms normTerms; const Double_t expIntegral { this->sumIkTerms(expTermIkVals_).real() }; normTerms.expTerm = expIntegral; const LauDecayTime::FuncType funcType { physModel_.getFunctionType() }; if ( funcType == LauDecayTime::FuncType::ExpTrig or funcType == LauDecayTime::FuncType::ExpHypTrig ) { const std::complex trigIntegral { this->sumIkTerms(trigTermIkVals_) }; normTerms.cosTerm = trigIntegral.real(); normTerms.sinTerm = trigIntegral.imag(); if ( funcType == LauDecayTime::FuncType::ExpTrig ) { normTerms.coshTerm = normTerms.expTerm; } else { const Double_t integralH { this->sumIkTerms(hypHTermIkVals_).real() }; const Double_t integralL { this->sumIkTerms(hypLTermIkVals_).real() }; const Double_t coshIntegral { 0.5 * (integralH + integralL) }; const Double_t sinhIntegral { 0.5 * (integralH - integralL) }; normTerms.coshTerm = coshIntegral; normTerms.sinhTerm = sinhIntegral; } } return normTerms; } template std::complex LauNonSmearedSplineEfficiencyDecayTimeIntegrator::sumIkTerms( const std::vector& cacheIk ) const { std::complex integral { 0.0 }; const std::size_t nSplineSegments { effModel_.nSegments() }; for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { - // eqn 30 in https://arxiv.org/pdf/1407.0748.pdf, using I_k from Appendix B.1 (with the typo corrected) + // eqn 30 in arXiv:1407.0748 using I_k from Appendix B.1 (with the typo corrected) for ( std::size_t k{0}; k < nCoeffs; ++k ) { integral += ( cacheIk[iSeg][k] * coeffs_[iSeg][k] ); } } return integral; } template void LauNonSmearedSplineEfficiencyDecayTimeIntegrator::propagateParUpdates() { // if any of the efficiency parameters changed in this iteration // we need to get new coefficients for each spline segment if ( effParamChanged_ ) { this->updateCoeffCache(); } // if any of the physics parameters changed in this iteration // we need to calculate new values of the Ik integrals if ( physicsParamChanged_ ) { this->updateIkCache(); } if ( effParamChanged_ or physicsParamChanged_ ) { normTerms_ = this->calcNormTerms(); } } template void LauNonSmearedSplineEfficiencyDecayTimeIntegrator::calcIntegrals( std::vector& cacheIk, const std::complex& u ) { const std::vector& tVals { effModel_.knotPositions() }; const std::size_t nSplineSegments { effModel_.nSegments() }; for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { const Double_t minAbs { tVals[iSeg] }; const Double_t maxAbs { tVals[iSeg+1] }; for(std::size_t k{0}; k < nCoeffs; ++k) { cacheIk[iSeg][k] = this->calcIk(k, minAbs, maxAbs, u); } } } template std::complex LauNonSmearedSplineEfficiencyDecayTimeIntegrator::calcIk( const std::size_t k, const Double_t minAbs, const Double_t maxAbs, const std::complex& u ) const { - // calculate I_k using https://arxiv.org/pdf/1407.0748.pdf Appendix B.1 (with n -> n+1 typo corrected - see below) + // calculate I_k using arXiv:1407.0748 Appendix B.1 (with n -> n+1 typo corrected - see below) // u = Gamma - i * DeltaM in general auto G = [&u,this](const std::size_t n){return -factorial_[n]/std::pow(u,n+1);}; //power of n+1 used rather than n, this is due to typo in the paper auto H = [&u](const std::size_t n, const Double_t t){return std::pow(t,n)*std::exp(-u*t);}; std::complex valIk { 0.0, 0.0 }; for (std::size_t j{0}; j <= k; ++j) {valIk += binomial_[k][j]*G(j)*( H( k-j, maxAbs ) - H( k-j, minAbs ) );} return valIk; } #endif diff --git a/inc/LauSmearedBinnedEfficiencyDecayTimeIntegrator.hh b/inc/LauSmearedBinnedEfficiencyDecayTimeIntegrator.hh index 5d648d7..18e84ab 100644 --- a/inc/LauSmearedBinnedEfficiencyDecayTimeIntegrator.hh +++ b/inc/LauSmearedBinnedEfficiencyDecayTimeIntegrator.hh @@ -1,85 +1,85 @@ /* 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 LauSmearedBinnedEfficiencyDecayTimeIntegrator.hh \brief File containing declaration of LauSmearedBinnedEfficiencyDecayTimeIntegrator class. */ #ifndef LAU_SMEARED_BINNEDEFFICIENCY_DECAYTIME_INTEGRATOR #define LAU_SMEARED_BINNEDEFFICIENCY_DECAYTIME_INTEGRATOR #include #include "Rtypes.h" #include "LauAbsDecayTimeIntegrator.hh" #include "LauBinnedDecayTimeEfficiency.hh" #include "LauDecayTimePhysicsModel.hh" #include "LauSmearedUniformEfficiencyDecayTimeIntegrator.hh" #include "LauUniformDecayTimeEfficiency.hh" class LauParameter; /*! \class LauSmearedBinnedEfficiencyDecayTimeIntegrator - \brief Class for defining the decay time integrator for the case of no resolution and uniform efficiency + \brief Class for defining the decay time integrator for the case of resolution and binned efficiency */ class LauSmearedBinnedEfficiencyDecayTimeIntegrator : public LauAbsDecayTimeIntegrator { public: //! Constructor LauSmearedBinnedEfficiencyDecayTimeIntegrator( const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const LauDecayTimePhysicsModel& physModel, const LauBinnedDecayTimeEfficiency& effModel, const LauDecayTimeResolution& resolModel ); //! Cache information from data /*! \param [in] abscissaErrors the values of the per-event decay time error */ void cacheInfo( const std::vector& abscissaErrors ) override; //! Propagate any updates to parameters and recalculate information as neeeded void propagateParUpdates() override; //! Retrieve the normalisation terms (optionally for a given event) /*! \param [in] iEvt the event index */ LauDecayTimeNormTerms getNormTerms( const std::size_t iEvt ) const override; //! Calculate the normalisation terms for a given value of the decay time error /*! \param [in] abscissaError the decay time error value (actually {1.0,err} to allow branchless control of scaling) */ LauDecayTimeNormTerms calcNormTerms( const std::array& abscissaError ) const override; private: //! The collection of uniform efficiency objects which we pass to the integrators std::vector efficiencies_; //! The collection of uniform efficiency integrator objects to which we delegate std::vector integrators_; ClassDefOverride(LauSmearedBinnedEfficiencyDecayTimeIntegrator, 0) }; #endif diff --git a/inc/LauSmearedSplineEfficiencyDecayTimeIntegrator.hh b/inc/LauSmearedSplineEfficiencyDecayTimeIntegrator.hh index 1b57fee..bd74019 100644 --- a/inc/LauSmearedSplineEfficiencyDecayTimeIntegrator.hh +++ b/inc/LauSmearedSplineEfficiencyDecayTimeIntegrator.hh @@ -1,866 +1,874 @@ /* 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 LauSmearedSplineEfficiencyDecayTimeIntegrator.hh \brief File containing declaration of LauSmearedSplineEfficiencyDecayTimeIntegrator class. */ #ifndef LAU_SMEARED_SPLINEEFFICIENCY_DECAYTIME_INTEGRATOR #define LAU_SMEARED_SPLINEEFFICIENCY_DECAYTIME_INTEGRATOR #include #include #include #include "RooMath.h" #include "Rtypes.h" #include "TSystem.h" #include "LauAbsDecayTimeIntegrator.hh" #include "LauConstants.hh" #include "LauDecayTime.hh" #include "LauDecayTimePhysicsModel.hh" #include "LauDecayTimeResolution.hh" #include "LauSplineDecayTimeEfficiency.hh" class LauParameter; /*! \class LauSmearedSplineEfficiencyDecayTimeIntegrator - \brief Class for defining the decay time integrator for the case of no resolution and uniform efficiency + \brief Class for defining the decay time integrator for the case of resolution and spline efficiency + + The implementation uses the formalism from arXiv:1407.0748 */ template class LauSmearedSplineEfficiencyDecayTimeIntegrator : public LauAbsDecayTimeIntegrator { public: //! Constructor LauSmearedSplineEfficiencyDecayTimeIntegrator( const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const LauDecayTimePhysicsModel& physModel, const LauSplineDecayTimeEfficiency& effModel, const LauDecayTimeResolution& resolModel ); //! Cache information from data /*! \param [in] abscissaErrors the values of the per-event decay time error */ void cacheInfo( const std::vector& abscissaErrors ) override; //! Propagate any updates to parameters and recalculate information as neeeded void propagateParUpdates() override; //! Retrieve the normalisation terms (optionally for a given event) /*! \param [in] iEvt the event index */ LauDecayTimeNormTerms getNormTerms( const std::size_t iEvt ) const override { return normTerms_[iEvt]; } //! Calculate the normalisation terms for a given value of the decay time error /*! \param [in] abscissaError the decay time error value (actually {1.0,err} to allow branchless control of scaling) */ LauDecayTimeNormTerms calcNormTerms( const std::array& abscissaError ) const override; private: //! The number of coefficients of each spline segment static constexpr std::size_t nCoeffs { static_cast>(Order) + 1 }; //! Define the RealArray typedef using RealArray = std::array; //! Define the Real2DArray typedef using Real2DArray = std::array; //! Define the ComplexArray typedef using ComplexArray = std::array,nCoeffs>; //! Update the cache of the spline coefficients void updateCoeffCache(); //! Update the cache of the powers of the mean and width void updateMeanAndWidthPowersCache(); //! Update the cache of the K and M vectors for each power of the polynomial /*! \param [in] forceUpdate used to force an update the first time around or if a resolution parameter changes */ void updateKvecMvecCache( const bool forceUpdate = false ); //! Calculate the normalisation terms for a given event /*! \param [in] iEvt the event index */ LauDecayTimeNormTerms calcNormTerms( const std::size_t iEvt ) const; //! Sum over the polynomial orders of the coefficients multiplied by the supplied integral components /*! \param [in] coeffs the spline coefficients \param [in] kVec the cached K-vector \param [in] mVec the cached M-vector \param [in] widthPowers the cached powers of the Gaussian width \param [in] meanPowers the cached powers of the Gaussian mean */ std::complex sumIkTerms( const RealArray& coeffs, const ComplexArray& kVec, const ComplexArray& mVec, const RealArray& widthPowers, const RealArray& meanPowers ) const; //! Sum each integral term for the given event, spline segment and Gaussian /*! \param [in] iEvt the index of the event \param [in] iSeg the index of the spline segment \param [in] iGauss the index of the Gaussian */ LauDecayTimeNormTerms sumIntegrals( const std::size_t iEvt, const std::size_t iSeg, const std::size_t iGauss ) const; //! Calculate the mean powers vector for the given Gaussian /*! \param [in] abscissaError the decay time error value (actually {1.0,err} to allow branchless control of scaling) \param [in] iGauss the index of the Gaussian */ RealArray generateMeanPowers(const std::array& abscissaError, const std::size_t iGauss) const; //! Calculate the width powers vector for the given Gaussian /*! \param [in] abscissaError the decay time error value (actually {1.0,err} to allow branchless control of scaling) \param [in] iGauss the index of the Gaussian */ RealArray generateWidthPowers(const std::array& abscissaError, const std::size_t iGauss) const; //! Calculate the K vector for the given value of z /*! \param [in] z = (Gamma - i Delta m)*sigma/sqrt(2) in general */ ComplexArray generateKvector(const std::complex& z) const; //! Calculate the M vector for the given value of z in the given range /*! \param [in] minAbs the lower end of the integration range \param [in] maxAbs the upper end of the integration range \param [in] z = (Gamma - i Delta m)*sigma/sqrt(2) in general \param [in] widthOverRoot2 the Gaussian width / sqrt(2) \param [in] mean the Gaussian mean */ ComplexArray generateMvector(const Double_t minAbs, const Double_t maxAbs, const std::complex& z, const Double_t widthOverRoot2, const Double_t mean) const; //! Calculate the M vector for the given value of z at one of the limits /*! \param [in] abs the value of the limit \param [in] z = (Gamma - i Delta m)*sigma/sqrt(2) in general \param [in] widthOverRoot2 the Gaussian width / sqrt(2) \param [in] mean the Gaussian mean */ ComplexArray generateMvectorLimit(const Double_t abs, const std::complex& z, const Double_t widthOverRoot2, const Double_t mean) const; //! The minimum value of the decay time const Double_t minAbscissa_; //! The maximum value of the decay time const Double_t maxAbscissa_; //! The physics model const LauDecayTimePhysicsModel& physModel_; //! The efficiency model const LauSplineDecayTimeEfficiency& effModel_; //! The resolution model const LauDecayTimeResolution& resolModel_; //! Cached value of lifetime const Double_t& tauVal_; //! Cached value of 1/lifetime const Double_t& gammaVal_; //! Cached value of mass difference const Double_t& deltaMVal_; //! Cached value of width difference const Double_t& deltaGammaVal_; //! Cached value of prompt fraction const Double_t& fracPromptVal_; //! Cached value of the number of Gaussians in the resolution model const std::size_t& nGauss_; //! Cached value of the Gaussian fractions const std::vector& fractionVals_; //! Cached value of the Gaussian means const std::vector& meanVals_; //! Cached value of the Gaussian widths const std::vector& widthVals_; //! Cached value of the mean scaling flags const std::vector& scaleMeans_; //! Cached value of the width scaling flags const std::vector& scaleWidths_; //! Have any of the physics parameters changed in the latest fit iteration const bool& physicsParamChanged_; //! Has the lifetime changed in the latest fit iteration const bool& tauChanged_; //! Has the mass difference changed in the latest fit iteration const bool& deltaMChanged_; //! Has the width difference changed in the latest fit iteration const bool& deltaGammaChanged_; //! Has the width difference changed in the latest fit iteration const bool& fracPromptChanged_; //! Have any of the efficiency parameters changed in the latest fit iteration const bool& effParamChanged_; //! Have any of the resolution parameters changed in the latest fit iteration const bool& resoParamChanged_; //! Factorial values RealArray factorial_; //! Binomial coefficient values Real2DArray binomial_; //! Values of the spline coefficients for each spline segment std::vector coeffs_; //! 0th to Order powers of the Gaussian means /*! outer vector has nEvents entries inner vector has nGauss_ entries */ std::vector> meanPowerVals_; //! 1st to Order+1 powers of the Gaussian widths / sqrt(2) /*! outer vector has nEvents entries inner vector has nGauss_ entries */ std::vector> widthPowerVals_; //! K-vector for the exponential terms /*! outer vector has nEvents entries inner vector has nGauss_ entries */ std::vector> expTermKvecVals_; //! K-vector for the cosine and sine terms /*! outer vector has nEvents entries inner vector has nGauss_ entries */ std::vector> trigTermKvecVals_; //! K-vector for the B_H part of the cosh and sinh terms /*! outer vector has nEvents entries inner vector has nGauss_ entries */ std::vector> hypHTermKvecVals_; //! K-vector for the B_L part of the cosh and sinh terms /*! outer vector has nEvents entries inner vector has nGauss_ entries */ std::vector> hypLTermKvecVals_; //! M-vector for the exponential terms /*! outer vector has nEvents entries middle vector has nSplineSegments entries inner vector has nGauss_ entries */ std::vector>> expTermMvecVals_; //! M-vector for the cosine and sine terms /*! outer vector has nEvents entries middle vector has nSplineSegments entries inner vector has nGauss_ entries */ std::vector>> trigTermMvecVals_; //! M-vector for the B_H part of the cosh and sinh terms /*! outer vector has nEvents entries middle vector has nSplineSegments entries inner vector has nGauss_ entries */ std::vector>> hypHTermMvecVals_; //! M-vector for the B_L part of the cosh and sinh terms /*! outer vector has nEvents entries middle vector has nSplineSegments entries inner vector has nGauss_ entries */ std::vector>> hypLTermMvecVals_; //! Cache of the decay time errors for each event /*! The inner array contains { 1.0, error } to allow un-branched control of whether or not to scale */ std::vector> abscissaErrors_; //! Cache of the normalisation terms for each event std::vector normTerms_; ClassDefOverride(LauSmearedSplineEfficiencyDecayTimeIntegrator, 0) }; templateClassImp(LauSmearedSplineEfficiencyDecayTimeIntegrator); template LauSmearedSplineEfficiencyDecayTimeIntegrator::LauSmearedSplineEfficiencyDecayTimeIntegrator( const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const LauDecayTimePhysicsModel& physModel, const LauSplineDecayTimeEfficiency& effModel, const LauDecayTimeResolution& resolModel ) : minAbscissa_{minAbscissaVal}, maxAbscissa_{maxAbscissaVal}, physModel_{physModel}, effModel_{effModel}, resolModel_{resolModel}, tauVal_{physModel_.getLifetimeValue()}, gammaVal_{physModel_.getGammaValue()}, deltaMVal_{physModel_.getDeltaMValue()}, deltaGammaVal_{physModel_.getDeltaGammaValue()}, fracPromptVal_{physModel_.getFracPromptValue()}, nGauss_{resolModel_.nGauss()}, fractionVals_{resolModel_.getFractionValues()}, meanVals_{resolModel_.getMeanValues()}, widthVals_{resolModel_.getWidthValues()}, scaleMeans_{resolModel_.scaleMeans()}, scaleWidths_{resolModel_.scaleWidths()}, physicsParamChanged_{physModel_.anythingChanged()}, tauChanged_{physModel_.lifetimeChanged()}, deltaMChanged_{physModel_.deltaMChanged()}, deltaGammaChanged_{physModel_.deltaGammaChanged()}, fracPromptChanged_{physModel_.fracPromptChanged()}, effParamChanged_{effModel_.anythingChanged()}, resoParamChanged_{resolModel_.anythingChanged()} { switch ( physModel_.getFunctionType() ) { case LauDecayTime::FuncType::Hist : // This does not make sense for us: // - Hist has no need of an integrator std::cerr << "ERROR in LauSmearedSplineEfficiencyDecayTimeIntegrator::LauSmearedSplineEfficiencyDecayTimeIntegrator : Unsupported function type in the physics model" << std::endl; gSystem->Exit(EXIT_FAILURE); break; case LauDecayTime::FuncType::Delta : case LauDecayTime::FuncType::DeltaExp : // TODO These are not yet implemented std::cerr << "ERROR in LauSmearedSplineEfficiencyDecayTimeIntegrator::LauSmearedSplineEfficiencyDecayTimeIntegrator : Function type not yet supported" << std::endl; gSystem->Exit(EXIT_FAILURE); break; case LauDecayTime::FuncType::Exp : case LauDecayTime::FuncType::ExpTrig : case LauDecayTime::FuncType::ExpHypTrig : // All fine, we can deal with these break; } // TODO - we should check that the range of the spline matches the range we've been given // populate the binomial and factorial arrays for ( std::size_t i{0}; i < nCoeffs; ++i ) { factorial_[i] = TMath::Factorial(i); for ( std::size_t j{0}; j < nCoeffs; ++j ) { if ( j <= i ) { binomial_[i][j] = TMath::Binomial(i,j); } else { binomial_[i][j] = 0.0; } } } } template void LauSmearedSplineEfficiencyDecayTimeIntegrator::cacheInfo( const std::vector& abscissaErrors ) { abscissaErrors_.clear(); if ( abscissaErrors.empty() ) { if ( resolModel_.scaleWithPerEventError() ) { std::cerr << "ERROR in LauSmearedSplineEfficiencyDecayTimeIntegrator::cacheInfo : No per-event decay time errors provided" << std::endl; gSystem->Exit(EXIT_FAILURE); } abscissaErrors_.push_back( { 1.0, 1.0 } ); } else { abscissaErrors_.reserve( abscissaErrors.size() ); for ( auto error : abscissaErrors ) { abscissaErrors_.push_back( { 1.0, error } ); } } const std::size_t nSplineSegments { effModel_.nSegments() }; coeffs_.clear(); coeffs_.resize(nSplineSegments); const std::size_t nEvents { abscissaErrors_.size() }; meanPowerVals_.clear(); meanPowerVals_.resize(nEvents); for ( auto& innerVec : meanPowerVals_ ) { innerVec.resize(nGauss_); } widthPowerVals_.clear(); widthPowerVals_.resize(nEvents); for ( auto& innerVec : widthPowerVals_ ) { innerVec.resize(nGauss_); } expTermKvecVals_.clear(); expTermKvecVals_.resize(nEvents); for ( auto& innerVec : expTermKvecVals_ ) { innerVec.resize(nGauss_); } trigTermKvecVals_.clear(); trigTermKvecVals_.resize(nEvents); for ( auto& innerVec : trigTermKvecVals_ ) { innerVec.resize(nGauss_); } hypHTermKvecVals_.clear(); hypHTermKvecVals_.resize(nEvents); for ( auto& innerVec : hypHTermKvecVals_ ) { innerVec.resize(nGauss_); } hypLTermKvecVals_.clear(); hypLTermKvecVals_.resize(nEvents); for ( auto& innerVec : hypLTermKvecVals_ ) { innerVec.resize(nGauss_); } // TODO - reorder iSeg and iGauss in M-vectors? what's better for vectorisation? expTermMvecVals_.clear(); expTermMvecVals_.resize(nEvents); for ( auto& middleVec : expTermMvecVals_) { middleVec.resize(nSplineSegments); for ( auto& innerVec : middleVec ) { innerVec.resize(nGauss_); } } trigTermMvecVals_.clear(); trigTermMvecVals_.resize(nEvents); for ( auto& middleVec : trigTermMvecVals_) { middleVec.resize(nSplineSegments); for ( auto& innerVec : middleVec ) { innerVec.resize(nGauss_); } } hypHTermMvecVals_.clear(); hypHTermMvecVals_.resize(nEvents); for ( auto& middleVec : hypHTermMvecVals_) { middleVec.resize(nSplineSegments); for ( auto& innerVec : middleVec ) { innerVec.resize(nGauss_); } } hypLTermMvecVals_.clear(); hypLTermMvecVals_.resize(nEvents); for ( auto& middleVec : hypLTermMvecVals_) { middleVec.resize(nSplineSegments); for ( auto& innerVec : middleVec ) { innerVec.resize(nGauss_); } } this->updateCoeffCache(); this->updateMeanAndWidthPowersCache(); this->updateKvecMvecCache(true); normTerms_.clear(); normTerms_.resize( nEvents ); for ( std::size_t iEvt{0}; iEvt < nEvents; ++iEvt ) { normTerms_[iEvt] = this->calcNormTerms( iEvt ); } } template void LauSmearedSplineEfficiencyDecayTimeIntegrator::propagateParUpdates() { // if any of the efficiency parameters changed in this iteration // we need to get new coefficients for each spline segment if ( effParamChanged_ ) { this->updateCoeffCache(); } // if any of the resolution parameters changed in this iteration // we need to update the powers of the means and widths if ( resoParamChanged_ ) { this->updateMeanAndWidthPowersCache(); } // if any of the physics or resolution parameters changed in this iteration // we need to calculate new values of the K- and M- vectors if ( physicsParamChanged_ or resoParamChanged_ ) { this->updateKvecMvecCache(resoParamChanged_); } if ( effParamChanged_ or physicsParamChanged_ or resoParamChanged_ ) { const std::size_t nEvents { abscissaErrors_.size() }; for ( std::size_t iEvt{0}; iEvt < nEvents; ++iEvt ) { normTerms_[iEvt] = this->calcNormTerms( iEvt ); } } } template void LauSmearedSplineEfficiencyDecayTimeIntegrator::updateCoeffCache() { const std::size_t nSplineSegments { effModel_.nSegments() }; for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { coeffs_[iSeg] = effModel_.getCoefficients(iSeg); } } template void LauSmearedSplineEfficiencyDecayTimeIntegrator::updateMeanAndWidthPowersCache() { // Calculate powers of mean and width/sqrt(2) needed by all terms in the smeared spline normalisation const std::size_t nEvents { abscissaErrors_.size() }; for ( std::size_t iEvt{0}; iEvt < nEvents; ++iEvt ) { for ( std::size_t iGauss{0}; iGauss < nGauss_; ++iGauss ) { meanPowerVals_[iEvt][iGauss] = this->generateMeanPowers( abscissaErrors_[iEvt], iGauss ); widthPowerVals_[iEvt][iGauss] = this->generateWidthPowers( abscissaErrors_[iEvt], iGauss ); } } } template typename LauSmearedSplineEfficiencyDecayTimeIntegrator::RealArray LauSmearedSplineEfficiencyDecayTimeIntegrator::generateMeanPowers(const std::array& abscissaError, const std::size_t iGauss) const { const Double_t mean { meanVals_[iGauss] * abscissaError[scaleMeans_[iGauss]] }; const Double_t meanSq { mean * mean }; RealArray meanPowers; meanPowers[0] = 1.0; meanPowers[1] = mean; meanPowers[2] = meanSq; meanPowers[3] = mean*meanSq; if constexpr ( Order == LauSplineOrder::Sixth ) { meanPowers[4] = meanPowers[2]*meanPowers[2]; meanPowers[5] = meanPowers[2]*meanPowers[3]; meanPowers[6] = meanPowers[3]*meanPowers[3]; } return meanPowers; } template typename LauSmearedSplineEfficiencyDecayTimeIntegrator::RealArray LauSmearedSplineEfficiencyDecayTimeIntegrator::generateWidthPowers(const std::array& abscissaError, const std::size_t iGauss) const { const Double_t widthOverRoot2 { widthVals_[iGauss] * abscissaError[scaleWidths_[iGauss]] / LauConstants::root2 }; const Double_t widthOverRoot2Sq { widthOverRoot2 * widthOverRoot2 }; RealArray widthPowers; widthPowers[0] = widthOverRoot2; widthPowers[1] = widthOverRoot2Sq; widthPowers[2] = widthOverRoot2*widthOverRoot2Sq; widthPowers[3] = widthOverRoot2Sq*widthOverRoot2Sq; if constexpr ( Order == LauSplineOrder::Sixth ) { widthPowers[4] = widthPowers[1]*widthPowers[2]; widthPowers[5] = widthPowers[2]*widthPowers[2]; widthPowers[6] = widthPowers[2]*widthPowers[3]; } return widthPowers; } template typename LauSmearedSplineEfficiencyDecayTimeIntegrator::ComplexArray LauSmearedSplineEfficiencyDecayTimeIntegrator::generateKvector(const std::complex& z) const { + // Expressions for K from section 4.1 in arXiv:1407.0748 + const std::complex zr { 1.0/z }; const std::complex zr2 { zr*zr }; ComplexArray K; K[0] = 0.5*zr; K[1] = 0.5*zr2; K[2] = zr*(1.0+zr2); K[3] = 3.0*zr2*(1.0+zr2); if constexpr ( Order == LauSplineOrder::Sixth ) { K[4] = 6.0*zr*(2.0*zr2*(1.0+zr2)+1.0); K[5] = 30.0*zr2*(2.0*zr2*(1.0+zr2)+1.0); K[6] = 60.0*zr*(3.0*zr2*(2.0*zr2*(1.0+zr2)+1.0)+1.0); } return K; } template typename LauSmearedSplineEfficiencyDecayTimeIntegrator::ComplexArray LauSmearedSplineEfficiencyDecayTimeIntegrator::generateMvector(const Double_t minAbs, const Double_t maxAbs, const std::complex& z, const Double_t widthOverRoot2, const Double_t mean) const { + // Expressions for M = J(x2) - J(x1) from section 4.1 in arXiv:1407.0748 + ComplexArray minM { this->generateMvectorLimit( minAbs, z, widthOverRoot2, mean ) }; ComplexArray maxM { this->generateMvectorLimit( maxAbs, z, widthOverRoot2, mean ) }; ComplexArray M; for ( std::size_t i{0}; i < nCoeffs; ++i ) { M[i] = maxM[i] - minM[i]; } return M; } template typename LauSmearedSplineEfficiencyDecayTimeIntegrator::ComplexArray LauSmearedSplineEfficiencyDecayTimeIntegrator::generateMvectorLimit(const Double_t abs, const std::complex& z, const Double_t widthOverRoot2, const Double_t mean) const { + // Expressions for J from section 4.1 in arXiv:1407.0748 + using namespace std::complex_literals; const Double_t x { (abs - mean) / (2.0 * widthOverRoot2) }; const Double_t x2 { x * x }; const Double_t ex2 { TMath::Exp(-(x2)) }; const std::complex arg { 1i * (z - x) }; //fad = the faddeeva term times the ex2 value (done in different ways depending on the domain) std::complex fad; if ( arg.imag() < -5.0 ) { fad = std::exp(-(x2) - (arg*arg)) * RooMath::erfc(-1i * arg); } else { fad = ex2 * RooMath::faddeeva(arg); } const Double_t invRootPi { 1.0 / LauConstants::rootPi }; ComplexArray M; M[0] = RooMath::erf(x) - fad; M[1] = -2.0 * (invRootPi*ex2 + x*fad); M[2] = -2.0 * (2.0*x*invRootPi*ex2 + (2.0*x2 - 1.0)*fad); M[3] = -4.0 * ((2.0*x2 - 1.0)*invRootPi*ex2 + x*(2.0*x2-3.0)*fad); if constexpr ( Order == LauSplineOrder::Sixth ) { const Double_t x3 { x * x2 }; M[4] = -4.0 * ((2.0*x2 - 3.0)*2.0*x*ex2*invRootPi + (3.0+4.0*x2*(x2-3.0))*fad); M[5] = -8.0 * ((3.0 + (x2 - 3.0)*4.0*x2)*ex2*invRootPi + (15.0 + (x2-5.0)*4.0*x2)*x*fad); M[6] = -8.0 * ((15.0 + (x2 - 5.0)*4.0*x2)*2.0*x*ex2*invRootPi + (8.0*x3*x3 - 60.0*x2*x2 + 90.0*x2 - 15.0)*fad); } return M; } template void LauSmearedSplineEfficiencyDecayTimeIntegrator::updateKvecMvecCache( const bool forceUpdate ) { const std::size_t nEvents { abscissaErrors_.size() }; const std::size_t nSplineSegments { effModel_.nSegments() }; const std::vector& tVals { effModel_.knotPositions() }; std::complex z { 0.0, 0.0 }; if ( tauChanged_ or forceUpdate ) { for ( std::size_t iEvt{0}; iEvt < nEvents; ++iEvt ) { for ( std::size_t iGauss{0}; iGauss < nGauss_; ++iGauss ) { const Double_t widthOverRoot2 { widthPowerVals_[iEvt][iGauss][0] }; const Double_t mean { meanPowerVals_[iEvt][iGauss][1] }; z = gammaVal_ * widthOverRoot2; expTermKvecVals_[iEvt][iGauss] = this->generateKvector(z); for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { const Double_t minAbs { tVals[iSeg] }; const Double_t maxAbs { tVals[iSeg+1] }; expTermMvecVals_[iEvt][iSeg][iGauss] = this->generateMvector(minAbs, maxAbs, z, widthOverRoot2, mean); } } } } const LauDecayTime::FuncType funcType { physModel_.getFunctionType() }; if ( funcType == LauDecayTime::FuncType::ExpTrig or funcType == LauDecayTime::FuncType::ExpHypTrig ) { if ( tauChanged_ or deltaMChanged_ or forceUpdate ) { for ( std::size_t iEvt{0}; iEvt < nEvents; ++iEvt ) { for ( std::size_t iGauss{0}; iGauss < nGauss_; ++iGauss ) { const Double_t widthOverRoot2 { widthPowerVals_[iEvt][iGauss][0] }; const Double_t mean { meanPowerVals_[iEvt][iGauss][1] }; z.real( gammaVal_ * widthOverRoot2 ); z.imag( -deltaMVal_ * widthOverRoot2 ); trigTermKvecVals_[iEvt][iGauss] = this->generateKvector(z); for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { const Double_t minAbs { tVals[iSeg] }; const Double_t maxAbs { tVals[iSeg+1] }; trigTermMvecVals_[iEvt][iSeg][iGauss] = this->generateMvector(minAbs, maxAbs, z, widthOverRoot2, mean); } } } } if ( funcType == LauDecayTime::FuncType::ExpHypTrig ) { if ( tauChanged_ or deltaGammaChanged_ or forceUpdate ) { for ( std::size_t iEvt{0}; iEvt < nEvents; ++iEvt ) { for ( std::size_t iGauss{0}; iGauss < nGauss_; ++iGauss ) { const Double_t widthOverRoot2 { widthPowerVals_[iEvt][iGauss][0] }; const Double_t mean { meanPowerVals_[iEvt][iGauss][1] }; z.real( ( gammaVal_ - 0.5 * deltaGammaVal_ ) * widthOverRoot2 );; z.imag( 0.0 ); hypHTermKvecVals_[iEvt][iGauss] = this->generateKvector(z); for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { const Double_t minAbs { tVals[iSeg] }; const Double_t maxAbs { tVals[iSeg+1] }; hypHTermMvecVals_[iEvt][iSeg][iGauss] = this->generateMvector(minAbs, maxAbs, z, widthOverRoot2, mean); } z.real( ( gammaVal_ + 0.5 * deltaGammaVal_ ) * widthOverRoot2 ); z.imag( 0.0 ); hypLTermKvecVals_[iEvt][iGauss] = this->generateKvector(z); for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { const Double_t minAbs { tVals[iSeg] }; const Double_t maxAbs { tVals[iSeg+1] }; hypLTermMvecVals_[iEvt][iSeg][iGauss] = this->generateMvector(minAbs, maxAbs, z, widthOverRoot2, mean); } } } } } } } template LauDecayTimeNormTerms LauSmearedSplineEfficiencyDecayTimeIntegrator::calcNormTerms( const std::array& abscissaError ) const { LauDecayTimeNormTerms normTerms; const LauDecayTime::FuncType funcType { physModel_.getFunctionType() }; const std::size_t nSplineSegments { effModel_.nSegments() }; const std::vector& tVals { effModel_.knotPositions() }; for ( std::size_t iGauss{0}; iGauss < nGauss_; ++iGauss ) { const RealArray meanPowers { this->generateMeanPowers( abscissaError, iGauss ) }; const RealArray widthPowers { this->generateWidthPowers( abscissaError, iGauss ) }; const Double_t mean { meanPowers[1] }; const Double_t widthOverRoot2 { widthPowers[0] }; const Double_t fraction { fractionVals_[iGauss] }; const std::complex zExp { gammaVal_ * widthOverRoot2 }; const std::complex zTrig { gammaVal_ * widthOverRoot2, -deltaMVal_ * widthOverRoot2 }; const std::complex zHypH { ( gammaVal_ - 0.5 * deltaGammaVal_ ) * widthOverRoot2 }; const std::complex zHypL { ( gammaVal_ + 0.5 * deltaGammaVal_ ) * widthOverRoot2 }; ComplexArray expTermKvec { this->generateKvector(zExp) }; ComplexArray trigTermKvec; ComplexArray hypHTermKvec; ComplexArray hypLTermKvec; if ( funcType == LauDecayTime::FuncType::ExpTrig or funcType == LauDecayTime::FuncType::ExpHypTrig ) { trigTermKvec = this->generateKvector(zTrig); if ( funcType == LauDecayTime::FuncType::ExpHypTrig ) { hypHTermKvec = this->generateKvector(zHypH); hypLTermKvec = this->generateKvector(zHypL); } } for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { const Double_t minAbs { tVals[iSeg] }; const Double_t maxAbs { tVals[iSeg+1] }; ComplexArray expTermMvec { this->generateMvector(minAbs, maxAbs, zExp, widthOverRoot2, mean) }; ComplexArray trigTermMvec; ComplexArray hypHTermMvec; ComplexArray hypLTermMvec; if ( funcType == LauDecayTime::FuncType::ExpTrig or funcType == LauDecayTime::FuncType::ExpHypTrig ) { trigTermMvec = this->generateMvector(minAbs, maxAbs, zTrig, widthOverRoot2, mean); if ( funcType == LauDecayTime::FuncType::ExpHypTrig ) { hypHTermMvec = this->generateMvector(minAbs, maxAbs, zHypH, widthOverRoot2, mean); hypLTermMvec = this->generateMvector(minAbs, maxAbs, zHypL, widthOverRoot2, mean); } } const Double_t expIntegral { this->sumIkTerms(coeffs_[iSeg], expTermKvec, expTermMvec, widthPowers, meanPowers).real() }; normTerms.expTerm += fraction * expIntegral; if ( funcType == LauDecayTime::FuncType::ExpTrig or funcType == LauDecayTime::FuncType::ExpHypTrig ) { const std::complex trigIntegral { this->sumIkTerms(coeffs_[iSeg], trigTermKvec, trigTermMvec, widthPowers, meanPowers) }; normTerms.cosTerm += fraction * trigIntegral.real(); normTerms.sinTerm += fraction * trigIntegral.imag(); if ( funcType == LauDecayTime::FuncType::ExpTrig ) { normTerms.coshTerm += fraction * expIntegral; } else { const Double_t integralH { this->sumIkTerms(coeffs_[iSeg], hypHTermKvec, hypHTermMvec, widthPowers, meanPowers).real() }; const Double_t integralL { this->sumIkTerms(coeffs_[iSeg], hypLTermKvec, hypLTermMvec, widthPowers, meanPowers).real() }; const Double_t coshIntegral { 0.5 * (integralH + integralL) }; const Double_t sinhIntegral { 0.5 * (integralH - integralL) }; normTerms.coshTerm += fraction * coshIntegral; normTerms.sinhTerm += fraction * sinhIntegral; } } } } return normTerms; } template LauDecayTimeNormTerms LauSmearedSplineEfficiencyDecayTimeIntegrator::calcNormTerms( const std::size_t iEvt ) const { LauDecayTimeNormTerms normTerms; // Sum the integral for each spline segment and each Gaussian const std::size_t nSplineSegments { effModel_.nSegments() }; for ( std::size_t iSeg{0}; iSeg < nSplineSegments; ++iSeg ) { for ( std::size_t iGauss{0}; iGauss < nGauss_; ++iGauss ) { normTerms += this->sumIntegrals( iEvt, iSeg, iGauss ); } } return normTerms; } template LauDecayTimeNormTerms LauSmearedSplineEfficiencyDecayTimeIntegrator::sumIntegrals( const std::size_t iEvt, const std::size_t iSeg, const std::size_t iGauss ) const { LauDecayTimeNormTerms normTerms; const Double_t fraction { fractionVals_[iGauss] }; const Double_t expIntegral { this->sumIkTerms(coeffs_[iSeg], expTermKvecVals_[iEvt][iGauss], expTermMvecVals_[iEvt][iSeg][iGauss], widthPowerVals_[iEvt][iGauss], meanPowerVals_[iEvt][iGauss]).real() }; normTerms.expTerm = fraction * expIntegral; const LauDecayTime::FuncType funcType { physModel_.getFunctionType() }; if ( funcType == LauDecayTime::FuncType::ExpTrig or funcType == LauDecayTime::FuncType::ExpHypTrig ) { const std::complex trigIntegral { this->sumIkTerms(coeffs_[iSeg], trigTermKvecVals_[iEvt][iGauss], trigTermMvecVals_[iEvt][iSeg][iGauss], widthPowerVals_[iEvt][iGauss], meanPowerVals_[iEvt][iGauss]) }; normTerms.cosTerm = fraction * trigIntegral.real(); normTerms.sinTerm = fraction * trigIntegral.imag(); if ( funcType == LauDecayTime::FuncType::ExpTrig ) { normTerms.coshTerm = normTerms.expTerm; } else { const Double_t integralH { this->sumIkTerms(coeffs_[iSeg], hypHTermKvecVals_[iEvt][iGauss], hypHTermMvecVals_[iEvt][iSeg][iGauss], widthPowerVals_[iEvt][iGauss], meanPowerVals_[iEvt][iGauss]).real() }; const Double_t integralL { this->sumIkTerms(coeffs_[iSeg], hypLTermKvecVals_[iEvt][iGauss], hypLTermMvecVals_[iEvt][iSeg][iGauss], widthPowerVals_[iEvt][iGauss], meanPowerVals_[iEvt][iGauss]).real() }; const Double_t coshIntegral { 0.5 * (integralH + integralL) }; const Double_t sinhIntegral { 0.5 * (integralH - integralL) }; normTerms.coshTerm = fraction * coshIntegral; normTerms.sinhTerm = fraction * sinhIntegral; } } return normTerms; } template std::complex LauSmearedSplineEfficiencyDecayTimeIntegrator::sumIkTerms( const RealArray& coeffs, const ComplexArray& kVec, const ComplexArray& mVec, const RealArray& widthPowers, const RealArray& meanPowers ) const { - // Triple sum to get N (eqn 31 and 29 in https://arxiv.org/pdf/1407.0748.pdf) + // Triple sum to get N - eqns 31 and 29 in arXiv:1407.0748 std::complex integral { 0.0 }; for ( std::size_t k{0}; k < nCoeffs; ++k ) { for ( std::size_t n{0}; n <= k; ++n ) { for ( std::size_t i{0}; i <= n; ++i ) { //The binomial coefficient terms const Double_t b { binomial_[n][i]*binomial_[k][n] }; integral += widthPowers[n]*coeffs[k]*meanPowers[k-n]*kVec[i]*mVec[n-i]*b; }}} return integral; } #endif diff --git a/inc/LauSmearedUniformEfficiencyDecayTimeIntegrator.hh b/inc/LauSmearedUniformEfficiencyDecayTimeIntegrator.hh index 7f00151..e26f6f2 100644 --- a/inc/LauSmearedUniformEfficiencyDecayTimeIntegrator.hh +++ b/inc/LauSmearedUniformEfficiencyDecayTimeIntegrator.hh @@ -1,155 +1,155 @@ /* 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 LauSmearedUniformEfficiencyDecayTimeIntegrator.hh \brief File containing declaration of LauSmearedUniformEfficiencyDecayTimeIntegrator class. */ #ifndef LAU_SMEARED_UNIFORMEFFICIENCY_DECAYTIME_INTEGRATOR #define LAU_SMEARED_UNIFORMEFFICIENCY_DECAYTIME_INTEGRATOR #include #include #include "Rtypes.h" #include "LauAbsDecayTimeIntegrator.hh" #include "LauDecayTimePhysicsModel.hh" #include "LauDecayTimeResolution.hh" #include "LauUniformDecayTimeEfficiency.hh" class LauParameter; /*! \class LauSmearedUniformEfficiencyDecayTimeIntegrator - \brief Class for defining the decay time integrator for the case of no resolution and uniform efficiency + \brief Class for defining the decay time integrator for the case of resolution and uniform efficiency */ class LauSmearedUniformEfficiencyDecayTimeIntegrator : public LauAbsDecayTimeIntegrator { public: //! Constructor LauSmearedUniformEfficiencyDecayTimeIntegrator( const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const LauDecayTimePhysicsModel& physModel, const LauUniformDecayTimeEfficiency& effModel, const LauDecayTimeResolution& resolModel ); //! Cache information from data /*! \param [in] abscissaErrors the values of the per-event decay time error */ void cacheInfo( const std::vector& abscissaErrors ) override; //! Propagate any updates to parameters and recalculate information as neeeded void propagateParUpdates() override; //! Retrieve the normalisation terms (optionally for a given event) /*! \param [in] iEvt the event index */ LauDecayTimeNormTerms getNormTerms( const std::size_t iEvt ) const override { return normTerms_[iEvt]; } //! Calculate the normalisation terms for a given value of the decay time error /*! \param [in] abscissaError the decay time error value (actually {1.0,err} to allow branchless control of scaling) */ LauDecayTimeNormTerms calcNormTerms( const std::array& abscissaError ) const override; private: //! Calculate the integrals for one of the Gaussians for a given value of the decay time error /*! \param [in] abscissaError the decay time error value (actually {1.0,err} to allow branchless control of scaling) \param [in] iGauss the index of the Gaussian */ LauDecayTimeNormTerms calcIntegrals( const std::array& abscissaError, const std::size_t iGauss ) const; //! Calculate the integral for a given value of z, Gaussian width and mean /*! \param [in] z = (Gamma - i Delta m)*sigma/sqrt(2) in general \param [in] sigmaOverRoot2 = sigma/sqrt(2) \param [in] mu = Gaussian mean */ std::complex smearedGeneralIntegral(const std::complex& z, const Double_t sigmaOverRoot2, const Double_t mu) const; //! The minimum value of the decay time const Double_t minAbscissa_; //! The maximum value of the decay time const Double_t maxAbscissa_; //! The physics model const LauDecayTimePhysicsModel& physModel_; //! The efficiency model const LauUniformDecayTimeEfficiency& effModel_; //! The resolution model const LauDecayTimeResolution& resolModel_; //! Cached value of lifetime const Double_t& tauVal_; //! Cached value of 1/lifetime const Double_t& gammaVal_; //! Cached value of mass difference const Double_t& deltaMVal_; //! Cached value of width difference const Double_t& deltaGammaVal_; //! Cached value of prompt fraction const Double_t& fracPromptVal_; //! Cached value of the number of Gaussians in the resolution model const std::size_t& nGauss_; //! Cached value of the Gaussian fractions const std::vector& fractionVals_; //! Cached value of the Gaussian means const std::vector& meanVals_; //! Cached value of the Gaussian widths const std::vector& widthVals_; //! Cached value of the mean scaling flags const std::vector& scaleMeans_; //! Cached value of the width scaling flags const std::vector& scaleWidths_; //! Have any of the physics parameters changed in the latest fit iteration const bool& physicsParamChanged_; //! Have any of the resolution parameters changed in the latest fit iteration const bool& resoParamChanged_; //! Cache of the decay time errors for each event /*! The inner array contains { 1.0, error } to allow un-branched control of whether or not to scale */ std::vector> abscissaErrors_; //! Cache of the normalisation terms for each event std::vector normTerms_; ClassDefOverride(LauSmearedUniformEfficiencyDecayTimeIntegrator, 0) }; #endif