diff --git a/inc/LauAbsDecayTimeIntegrator.hh b/inc/LauAbsDecayTimeIntegrator.hh index 29acc9e..0880576 100644 --- a/inc/LauAbsDecayTimeIntegrator.hh +++ b/inc/LauAbsDecayTimeIntegrator.hh @@ -1,123 +1,124 @@ /* 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 LauAbsDecayTimeIntegrator.hh \brief File containing declaration of LauAbsDecayTimeIntegrator class. */ #ifndef LAU_ABS_DECAYTIME_INTEGRATOR #define LAU_ABS_DECAYTIME_INTEGRATOR +#include #include #include "Rtypes.h" class LauParameter; /*! \struct LauDecayTimeNormTerms \brief Struct for gathering the normalisation terms for the decay time model */ struct LauDecayTimeNormTerms { //! Normalisation of the exponential term Double_t expTerm{0.0}; //! Normalisation of the cos term Double_t cosTerm{0.0}; //! Normalisation of the sin term Double_t sinTerm{0.0}; //! Normalisation of the cosh term Double_t coshTerm{0.0}; //! Normalisation of the sinh term Double_t sinhTerm{0.0}; //! Define operator+= so we can sum these LauDecayTimeNormTerms& operator+=( const LauDecayTimeNormTerms& rhs ) { expTerm += rhs.expTerm; cosTerm += rhs.cosTerm; sinTerm += rhs.sinTerm; coshTerm += rhs.coshTerm; sinhTerm += rhs.sinhTerm; return *this; } //! Define operator*= so we can scale by a constant LauDecayTimeNormTerms& operator*=( const Double_t scale ) { expTerm *= scale; cosTerm *= scale; sinTerm *= scale; coshTerm *= scale; sinhTerm *= scale; return *this; } }; /*! \class LauAbsDecayTimeIntegrator \brief Class for defining the abstract interface for performing the integration of the decay time model */ class LauAbsDecayTimeIntegrator { public: //! Cache information from data /*! \param [in] abscissaErrors the values of the per-event decay time error */ virtual void cacheInfo( const std::vector& abscissaErrors ) = 0; //! Propagate any updates to parameters and recalculate information as neeeded virtual void propagateParUpdates() = 0; //! Retrieve the normalisation terms for a given event /*! \param [in] iEvt the event index */ virtual LauDecayTimeNormTerms getNormTerms( const std::size_t iEvt ) const = 0; //! 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) */ virtual LauDecayTimeNormTerms calcNormTerms( const std::array& abscissaError = {1.0,1.0} ) const = 0; //! Constructor LauAbsDecayTimeIntegrator() = default; //! Destructor virtual ~LauAbsDecayTimeIntegrator() = default; //! Copy constructor (deleted) LauAbsDecayTimeIntegrator(const LauAbsDecayTimeIntegrator& other) = delete; //! Copy assignment operator (deleted) LauAbsDecayTimeIntegrator& operator=(const LauAbsDecayTimeIntegrator& other) = delete; //! Move constructor (deleted) LauAbsDecayTimeIntegrator(LauAbsDecayTimeIntegrator&& other) = default; //! Move assignment operator (deleted) LauAbsDecayTimeIntegrator& operator=(LauAbsDecayTimeIntegrator&& other) = default; private: ClassDef(LauAbsDecayTimeIntegrator, 0) }; #endif diff --git a/src/LauNonSmearedDecayTimeCalculator.cc b/src/LauNonSmearedDecayTimeCalculator.cc index ccf5be7..0350f6d 100644 --- a/src/LauNonSmearedDecayTimeCalculator.cc +++ b/src/LauNonSmearedDecayTimeCalculator.cc @@ -1,119 +1,120 @@ /* 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 LauNonSmearedDecayTimeCalculator.cc \brief File containing implementation of LauNonSmearedDecayTimeCalculator class. */ +#include #include #include #include "TMath.h" #include "TSystem.h" #include "LauNonSmearedDecayTimeCalculator.hh" ClassImp(LauNonSmearedDecayTimeCalculator); LauNonSmearedDecayTimeCalculator::LauNonSmearedDecayTimeCalculator( const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const LauDecayTimePhysicsModel& physModel ) : minAbscissa_{minAbscissaVal}, maxAbscissa_{maxAbscissaVal}, physModel_{physModel}, tauVal_{physModel_.getLifetimeValue()}, gammaVal_{physModel_.getGammaValue()}, deltaMVal_{physModel_.getDeltaMValue()}, deltaGammaVal_{physModel_.getDeltaGammaValue()}, physicsParamChanged_{physModel_.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 LauNonSmearedDecayTimeCalculator::LauNonSmearedDecayTimeCalculator : 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; } } void LauNonSmearedDecayTimeCalculator::cacheInfo( const std::vector& abscissas, [[maybe_unused]] const std::vector& abscissaErrors ) { abscissas_ = abscissas; terms_.clear(); terms_.resize( abscissas_.size() ); this->updateCache(); } void LauNonSmearedDecayTimeCalculator::updateCache() { const std::size_t nEvents { abscissas_.size() }; for ( std::size_t iEvt{0}; iEvt < nEvents; ++iEvt ) { terms_[iEvt] = this->calcTerms( abscissas_[iEvt], {1.0,1.0} ); } } void LauNonSmearedDecayTimeCalculator::propagateParUpdates() { if ( physicsParamChanged_ ) { this->updateCache(); } } LauDecayTimeTerms LauNonSmearedDecayTimeCalculator::calcTerms( const Double_t abscissa, [[maybe_unused]] const std::array& abscissaError ) { LauDecayTimeTerms terms; terms.expTerm = TMath::Exp( -abscissa * gammaVal_ ); // Calculate also the terms related to cosine and sine const LauDecayTime::FuncType type { physModel_.getFunctionType() }; if ( type == LauDecayTime::FuncType::ExpTrig ) { terms.coshTerm = terms.expTerm; terms.sinhTerm = 0.0; terms.cosTerm = TMath::Cos( deltaMVal_ * abscissa ) * terms.expTerm; terms.sinTerm = TMath::Sin( deltaMVal_ * abscissa ) * terms.expTerm; } // Calculate also the terms related to cosh, sinh, cosine, and sine else if ( type == LauDecayTime::FuncType::ExpHypTrig ) { terms.coshTerm = TMath::CosH( 0.5 * deltaGammaVal_ * abscissa ) * terms.expTerm; terms.sinhTerm = TMath::SinH( 0.5 * deltaGammaVal_ * abscissa ) * terms.expTerm; terms.cosTerm = TMath::Cos( deltaMVal_ * abscissa ) * terms.expTerm; terms.sinTerm = TMath::Sin( deltaMVal_ * abscissa ) * terms.expTerm; } return terms; } diff --git a/src/LauSmearedDecayTimeCalculator.cc b/src/LauSmearedDecayTimeCalculator.cc index 6f1ff1e..76f4acb 100644 --- a/src/LauSmearedDecayTimeCalculator.cc +++ b/src/LauSmearedDecayTimeCalculator.cc @@ -1,222 +1,223 @@ /* 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 LauSmearedDecayTimeCalculator.cc \brief File containing implementation of LauSmearedDecayTimeCalculator class. */ +#include #include #include #include #include "RooMath.h" #include "TMath.h" #include "TSystem.h" #include "LauConstants.hh" #include "LauSmearedDecayTimeCalculator.hh" ClassImp(LauSmearedDecayTimeCalculator); LauSmearedDecayTimeCalculator::LauSmearedDecayTimeCalculator( const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const LauDecayTimePhysicsModel& physModel, const LauDecayTimeResolution& resolModel ) : minAbscissa_{minAbscissaVal}, maxAbscissa_{maxAbscissaVal}, physModel_{physModel}, 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()}, 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 LauSmearedDecayTimeCalculator::LauSmearedDecayTimeCalculator : 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 LauSmearedDecayTimeCalculator::LauSmearedDecayTimeCalculator : 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; } } void LauSmearedDecayTimeCalculator::cacheInfo( const std::vector& abscissas, const std::vector& abscissaErrors ) { abscissas_ = abscissas; terms_.clear(); terms_.resize( abscissas_.size() ); abscissaErrors_.clear(); if ( abscissaErrors.empty() ) { if ( resolModel_.scaleWithPerEventError() ) { std::cerr << "ERROR in LauSmearedDecayTimeCalculator::cacheInfo : No per-event decay time errors provided" << std::endl; gSystem->Exit(EXIT_FAILURE); } abscissaErrors_.push_back( { 1.0, 1.0 } ); } else { if ( abscissaErrors.size() != abscissas_.size() ) { std::cerr << "ERROR in LauSmearedDecayTimeCalculator::cacheInfo : Sizes of decay time and decay time error vectors do not match" << std::endl; std::cerr << " : decay time: " << abscissas_.size() << std::endl; std::cerr << " : decay time error: " << abscissaErrors.size() << std::endl; gSystem->Exit(EXIT_FAILURE); } abscissaErrors_.reserve( abscissaErrors.size() ); for ( auto error : abscissaErrors ) { abscissaErrors_.push_back( { 1.0, error } ); } } this->updateCache(); } void LauSmearedDecayTimeCalculator::updateCache() { const std::size_t nEvents { abscissas_.size() }; for ( std::size_t iEvt{0}; iEvt < nEvents; ++iEvt ) { terms_[iEvt] = this->calcTerms( abscissas_[iEvt], abscissaErrors_[iEvt] ); } } void LauSmearedDecayTimeCalculator::propagateParUpdates() { if ( physicsParamChanged_ or resoParamChanged_ ) { this->updateCache(); } } LauDecayTimeTerms LauSmearedDecayTimeCalculator::calcTerms( const Double_t abscissa, const std::array& abscissaError ) { /* // TODO - this is how to calculate the Delta and DeltaExp cases - looks like they already include their normalisation // - to be tidied up and integrated into the rest of the framework in due course Double_t value(0.0); if (type == LauDecayTime::FuncType::Delta || type == LauDecayTime::FuncType::DeltaExp) { // Calculate the gaussian function(s) for (std::size_t iGauss(0); iGauss 1e-10) { Double_t scale = LauConstants::root2*sigma[iGauss]; Double_t scale2 = LauConstants::rootPiBy2*sigma[iGauss]; Double_t x = ( abscissa - mean[iGauss] ) / scale; Double_t exponent = -x*x; Double_t norm = scale2*(TMath::Erf((maxAbscissa_ - mean[iGauss])/scale) - TMath::Erf((minAbscissa_ - mean[iGauss])/scale)); value += frac[iGauss]*TMath::Exp(exponent)/norm; } } // TODO - calc exp term here if (type == LauDecayTime::FuncType::DeltaExp) { value *= fracPrompt; value += (1.0-fracPrompt)*expTerm_; } else { value = expTerm_; } } */ LauDecayTimeTerms terms; const LauDecayTime::FuncType type { physModel_.getFunctionType() }; for ( std::size_t iGauss{0}; iGauss < nGauss_; ++iGauss ) { const Double_t mean { meanVals_[iGauss] * abscissaError[scaleMeans_[iGauss]] }; const Double_t widthOverRoot2 { widthVals_[iGauss] * abscissaError[scaleWidths_[iGauss]] / LauConstants::root2 }; const Double_t x { ( abscissa - mean ) / ( 2.0 * widthOverRoot2 ) }; const std::complex z { gammaVal_ * widthOverRoot2 }; const Double_t expTerm { this->smearedGeneralTerm( z, x ).real() }; terms.expTerm += fractionVals_[iGauss] * expTerm; if ( type == LauDecayTime::FuncType::ExpTrig or type == LauDecayTime::FuncType::ExpHypTrig ) { const std::complex zTrig { gammaVal_ * widthOverRoot2, -deltaMVal_ * widthOverRoot2 }; const std::complex trigTerm { this->smearedGeneralTerm( zTrig, x ) }; const Double_t cosTerm { trigTerm.real() }; const Double_t sinTerm { trigTerm.imag() }; terms.cosTerm += fractionVals_[iGauss] * cosTerm; terms.sinTerm += fractionVals_[iGauss] * sinTerm; if ( type == LauDecayTime::FuncType::ExpTrig ) { terms.coshTerm += fractionVals_[iGauss] * expTerm; } else { const std::complex zH { (gammaVal_ - 0.5 * deltaGammaVal_) * widthOverRoot2 }; const std::complex zL { (gammaVal_ + 0.5 * deltaGammaVal_) * widthOverRoot2 }; const Double_t termH { this->smearedGeneralTerm( zH, x ).real() }; const Double_t termL { this->smearedGeneralTerm( zL, x ).real() }; const Double_t coshTerm { 0.5 * (termH + termL) }; const Double_t sinhTerm { 0.5 * (termH - termL) }; terms.coshTerm += fractionVals_[iGauss] * coshTerm; terms.sinhTerm += fractionVals_[iGauss] * sinhTerm; } } } return terms; } std::complex LauSmearedDecayTimeCalculator::smearedGeneralTerm( const std::complex& z, const Double_t x ) const { using namespace std::complex_literals; const Double_t x2 { x * x }; const std::complex arg1 { 1i * (z - x) }; const std::complex arg2 { -(x2) - (arg1*arg1) }; const std::complex conv { ( arg1.imag() < -5.0 ) ? 0.5 * std::exp(arg2) * RooMath::erfc( -1i * arg1 ) : 0.5 * std::exp(-(x2)) * RooMath::faddeeva(arg1) }; return conv; }