diff --git a/inc/LauSplineDecayTimeEfficiency.hh b/inc/LauSplineDecayTimeEfficiency.hh index 21393f8..86cecf8 100644 --- a/inc/LauSplineDecayTimeEfficiency.hh +++ b/inc/LauSplineDecayTimeEfficiency.hh @@ -1,446 +1,449 @@ /* 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 #include "Rtypes.h" #include "TFitResult.h" #include "TMath.h" #include "TSystem.h" #include "Lau1DCubicSpline.hh" #include "LauAbsDecayTimeEfficiency.hh" #include "LauParameter.hh" class TH1; //! Defines the allowed orders of spline polynomials enum class LauSplineOrder : std::size_t { Cubic = 3, //!< 3rd order Sixth = 6 //!< 6th order }; template class LauSplineDecayTimeEfficiency : public LauAbsDecayTimeEfficiency { public: //! 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; //! 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; //! Find the maximum value of the efficiency in the range /*! Works analytically for cubic splines and numerically for 6th order splines \return the maximum of the spline */ Double_t findMaximum() 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; //! 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_; } //! Fit our spline to a TH1 /*! Useful to obtain reasonable starting values and uncertainties for the knot parameters \param [in] hist the histogram to be fit to \param [in] printLevel the level of printout desired from fit \param [in] doWL If true do a weighted likelihood fit, else chi2 */ void fitToTH1( TH1* hist, const LauOutputLevel printLevel = LauOutputLevel::Quiet, const bool doWL = false ); 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 ) { constexpr Double_t maxVal { (Order == LauSplineOrder::Cubic) ? 1.0 : 2.0 }; ownedParams_.emplace_back( std::make_unique( Form( "%s_knot_%lu", prefix_.Data(), index ), values_[index], 0.0, maxVal, 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_; } } } template void LauSplineDecayTimeEfficiency::fitToTH1( TH1* hist, const LauOutputLevel printLevel, const bool doWL ) { std::cout << "INFO in LauSplineDecayTimeEfficiency::fitToTH1 : fitting " << prefix_ << " spline to supplied histogram \"" << hist->GetName() << "\"\n"; std::cout << " : to obtain initial estimates of parameter values and errors" << std::endl; TFitResultPtr fitResult { effSpline_->fitToTH1( hist, printLevel, doWL ) }; if ( fitResult->Status() != 0 || fitResult->CovMatrixStatus() < 2 ) { std::cerr << "WARNING in LauSplineDecayTimeEfficiency::fitToTH1 : fit failed with status " << fitResult->Status() << " and covariance matrix status " << fitResult->CovMatrixStatus() << "\n"; std::cerr << " : will not update spline parameters based on this fit result" << std::endl; return; } if ( fitResult->CovMatrixStatus() == 2 ) { std::cerr << "WARNING in LauSplineDecayTimeEfficiency::fitToTH1 : final covariance matrix not positive definite\n"; std::cerr << " : continuing, but double check your fit!" << std::endl; } const std::size_t nKnots { ownedParams_.size() }; for ( std::size_t knot{0}; knot < nKnots; ++knot ) { const Double_t knotVal { fitResult->Value(knot) }; const Double_t knotErr { fitResult->Error(knot) }; std::cout << " : Setting parameter " << ownedParams_[knot]->name() << " to " << knotVal << " +/- " << knotErr << std::endl; ownedParams_[knot]->value( knotVal ); ownedParams_[knot]->genValue( knotVal ); ownedParams_[knot]->initValue( knotVal ); ownedParams_[knot]->error( knotErr ); } this->updateParameterCache(); } template Double_t LauSplineDecayTimeEfficiency::findMaximum() const { if constexpr ( Order == LauSplineOrder::Cubic ) { return effSpline_->findMaximum(); } else { //Build a lambda to construct a TF1 auto sixthOrder = [this](Double_t x[], [[maybe_unused]] Double_t p[]) {return effSpline_->evaluate(x[0])*otherSpline_->getEfficiency(x[0]);}; - TF1 func("function", sixthOrder, effSpline_->getXValues().front(), effSpline_->getXValues().back(), 0); - return func.GetMaximum(func.GetXmin(), func.GetXmax()); + const Double_t tMin {effSpline_->getXValues().front()}; + const Double_t tMax {effSpline_->getXValues().back()*(1.0-std::numeric_limits::epsilon())}; + TF1 func("function", sixthOrder, tMin, tMax, 0); + return func.GetMaximum(tMin, tMax); } } #endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 17a2050..fc6af6c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,18 +1,19 @@ list(APPEND TEST_SOURCES TestCovariant TestCovariant2 TestNewKinematicsMethods TestFitSplineToTH1 TestFitDoubleSplineToTH1 TestSplineDTAdivision TestWriteSplineToJson TestReadSplineFromJson + TestSplineFindMax ) foreach( _test ${TEST_SOURCES}) add_executable(${_test} ${_test}.cc) target_link_libraries(${_test} PRIVATE Laura++) - #install(TARGETS ${_test} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) + install(TARGETS ${_test} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) endforeach() diff --git a/test/TestSplineFindMax.cc b/test/TestSplineFindMax.cc new file mode 100644 index 0000000..663f809 --- /dev/null +++ b/test/TestSplineFindMax.cc @@ -0,0 +1,77 @@ +#include +#include +#include +#include + +#include "TFile.h" +#include "TH1.h" +#include "TRegexp.h" +#include "TString.h" + +#include "Lau1DCubicSpline.hh" +#include "LauSplineDecayTimeEfficiency.hh" + +int main() +{ + std::unique_ptr DTA_file { TFile::Open("root://eoslhcb.cern.ch//eos/lhcb/wg/b2oc/B2D0pipi/Tuples22/DTAhists.root") }; + if ( ! DTA_file || DTA_file->IsZombie() ) { + std::cerr << "ERROR : cannot open histogram file, exiting" << std::endl; + return 1; + } + + const std::vector dtvals { 0.20, 0.50, 0.75, 1.00, 2.00, 3.00, 4.00, 6.00, 15.00 }; + + std::map> splines; + + for ( const auto&& obj : *(DTA_file->GetListOfKeys()) ) { + const TString histName { obj->GetName() }; + if ( histName.Contains("Sideband") || histName.Contains("Run1") ) { + continue; + } + + TH1* thid { dynamic_cast( DTA_file->Get(histName) ) }; + if ( ! thid ) { + std::cerr << "ERROR : cannot find histogram : \"" << histName << "\", skipping..." << std::endl; + continue; + } + thid->SetDirectory(0); + + Lau1DCubicSpline::SplineType splineType = Lau1DCubicSpline::SplineType::StandardSpline; + if ( ! histName.BeginsWith("signalMC_Kpi") ) { + TString denomName = Form( "signalMC%s", histName("_.+$").Data() ); + denomName.ReplaceAll("_KK_", "_Kpi_"); + denomName.ReplaceAll("_pipi_","_Kpi_"); + + TH1* denom { dynamic_cast( DTA_file->Get(denomName) ) }; + if(not denom){std::cerr << "ERROR : Could not find denominator!" << std::endl; continue;} + thid->Divide( denom ); + delete denom; + } else { + splineType = Lau1DCubicSpline::SplineType::AkimaSpline; + } + + std::vector effvals; + effvals.assign(dtvals.size(), 1.0); + for ( std::size_t i{0}; i < effvals.size(); ++i ) { + effvals[i] = thid->GetBinContent(thid->FindFixBin(dtvals[i])); + } + + splines[histName] = std::make_unique(dtvals,effvals,splineType); + + delete thid; + } + + DTA_file->Close(); + + LauSplineDecayTimeEfficiency denom_eff { "signalMC_Kpi", std::move(splines["signalMC_Kpi_Run2_DTAhist"]) }; + splines.erase("signalMC_Kpi_Run2_DTAhist"); + + std::cout << "Denominator maximum = " << denom_eff.findMaximum() << std::endl; + + for ( auto& [ name, spline ] : splines ) { + LauSplineDecayTimeEfficiency eff { name, denom_eff, std::move(spline) }; + std::cout << name << " maximum = " << eff.findMaximum() << std::endl; + } + + return 0; +}