diff --git a/inc/LauDecayTimePhysicsModel.hh b/inc/LauDecayTimePhysicsModel.hh index e127f1e..55d2f75 100644 --- a/inc/LauDecayTimePhysicsModel.hh +++ b/inc/LauDecayTimePhysicsModel.hh @@ -1,329 +1,329 @@ /* 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 LauDecayTimePhysicsModel.hh \brief File containing declaration of LauDecayTimePhysicsModel class. */ #ifndef LAU_DECAYTIME_PHYSICSMODEL #define LAU_DECAYTIME_PHYSICSMODEL #include "LauDecayTime.hh" #include "LauJsonTools.hh" #include "Rtypes.h" #include "TString.h" #include #include #include #include #include class LauAbsRValue; /*! \class LauDecayTimePhysicsModel \brief Class for defining the physics and associated parameters for the decay time model */ class LauDecayTimePhysicsModel final { public: //! Constructor /*! \param type the type of the physics model \param parameters the parameters of the physics model */ LauDecayTimePhysicsModel( const LauDecayTime::FuncType type, std::vector> parameters ); //! Retrieve the type of the physics function /*! \return the function type */ LauDecayTime::FuncType getFunctionType() const { return type_; } //! Retrieve the parameters of the physics model so that they can be loaded into a fit /*! \return the parameters of the physics model */ std::vector getParameters() { return params_; } //! Retrieve the up-to-date value of the lifetime /*! \return the lifetime value */ const Double_t& getLifetimeValue() const { return tauVal_; } //! Retrieve the up-to-date value of the inverse lifetime /*! \return the inverse lifetime value */ const Double_t& getGammaValue() const { return gammaVal_; } //! Retrieve the up-to-date value of the mass difference /*! \return the mass difference value */ const Double_t& getDeltaMValue() const { return deltaMVal_; } //! Retrieve the up-to-date value of the width difference /*! \return the width difference value */ const Double_t& getDeltaGammaValue() const { return deltaGammaVal_; } //! Retrieve the up-to-date value of the prompt fraction /*! \return the prompt fraction value */ const Double_t& getFracPromptValue() const { return fracPromptVal_; } //! Retrieve whether any of the parameter values have changed in the last fit iteration /*! \return the any param changed flag */ const bool& anythingChanged() const { return anythingChanged_; } //! Retrieve whether the lifetime value has changed in the last fit iteration /*! \return the lifetime changed flag */ const bool& lifetimeChanged() const { return tauChanged_; } //! Retrieve the up-to-date value of the mass difference /*! \return the mass difference value */ const bool& deltaMChanged() const { return deltaMChanged_; } //! Retrieve the up-to-date value of the width difference /*! \return the width difference value */ const bool& deltaGammaChanged() const { return deltaGammaChanged_; } //! Retrieve the up-to-date value of the prompt fraction /*! \return the prompt fraction value */ const bool& fracPromptChanged() const { return fracPromptChanged_; } //! Initialise the parameter cache /*! Must be called prior to starting fitting or generation */ void initialise(); //! Propagate any updates to parameters and recalculate information as neeeded /*! Should be called at each fit iteration */ void propagateParUpdates(); //! Construct a collection of physics model objects based on values read from a JSON value /*! The JSON value should be an Object, which must contain the following elements: - "parameters", an Array, which contains an Object for each parameter of the models to be constructed - "signalModel", an Object, which contains: + "functionType", a String to specify the function type (must match a state of LauDecayTime::FuncType) + "parameters", an Array of Strings that specifies the parameters from the top-level "parameters" Array to be used to build the signal model - "backgroundModels", an Array, which contains Objects for each background component, each of which contain: + "name", a String specifying the background component name + "model", an Object to define the model (definition as per the "signalModel" entry above) \param [in] j the JSON value to deserialise \return the collection of newly constructed models */ static std::map> readFromJson( const nlohmann::json& j ); //! Construct a collection of physics model objects based on values read from a JSON file /*! The JSON value structure is as defined in readFromJson(nlohmann::json) \param [in] fileName the name of the file from which the JSON should be read \param [in] elementName the optional name of the JSON element that contains the definitions (defaults to using the root record) \return the collection of newly constructed models */ static std::map> readFromJson( const TString& fileName, const TString& elementName = "" ); //! Write a collection of physics model objects to a JSON file /*! \tparam PhysModelPtr a pointer-like type that points to LauDecayTimePhysicsModel objects \param [in] models the collection of models to be written out \param [in] fileName the name of the file to which the JSON should be written \param [in] elementName the (optional) name of the JSON element to which the coefficients should be written \param [in] append if true, append the spline to the existing JSON within the file - if using this option it is then mandatory to provide elementName \param [in] indent the indentation level to use in the output in number of spaces (defaults to 4) */ template static void writeToJson( const std::map& models, const TString& fileName, const TString& elementName = "", const bool append = false, const int indent = 4 ); //! Retrieve the specified parameter /*! \param [in] parName the parameter to retrieve */ const LauAbsRValue* findParameter(const TString& parName) const; //! Retrieve the specified parameter /*! \param [in] parName the parameter to retrieve */ LauAbsRValue* findParameter(const TString& parName); private: //! Update the cached parameter values void updateParameterCache(); - //! Read the signal model part of the JSON serialisation + //! Read the physics model part of the JSON serialisation /*! \param [in] j the JSON value from which to read the model \param [in] parameters the collection of parameters, some of which should be used to construct the model \return the newly constructed signal model */ - static std::unique_ptr readSignalModelFromJson( const nlohmann::json& j, std::vector>& parameters ); + static std::unique_ptr readPhysicsModelFromJson( const nlohmann::json& j, std::vector>& parameters ); //! Read the background models part of the JSON serialisation /*! \param [in] j the JSON value from which to read the model \param [in] parameters the collection of parameters, which should be used to construct the models \return the newly constructed background models */ static std::map> readBackgroundModelsFromJson( const nlohmann::json& j, std::vector>& parameters ); //! Which type of decay time function is this? const LauDecayTime::FuncType type_; //! Store of all parameters of the physics model std::vector> paramsOwned_; //! Store of all parameters of the physics model std::vector params_; //! Lifetime parameter LauAbsRValue* tau_{nullptr}; //! Mass difference parameter LauAbsRValue* deltaM_{nullptr}; //! Width difference parameter LauAbsRValue* deltaGamma_{nullptr}; //! Parameter for the fraction of prompt events in DeltaExp LauAbsRValue* fracPrompt_{nullptr}; //! Cached value of lifetime Double_t tauVal_{0.0}; //! Cached value of 1/lifetime Double_t gammaVal_{0.0}; //! Cached value of mass difference Double_t deltaMVal_{0.0}; //! Cached value of width difference Double_t deltaGammaVal_{0.0}; //! Cached value of prompt fraction Double_t fracPromptVal_{0.0}; //! Is the lifetime floating in the fit? bool tauFloating_{false}; //! Is the mass difference floating in the fit? bool deltaMFloating_{false}; //! Is the width difference floating in the fit? bool deltaGammaFloating_{false}; //! Is the prompt fraction floating in the fit? bool fracPromptFloating_{false}; //! Are any of the physics parameters floating in the fit? bool anythingFloating_{false}; //! Has the lifetime parameter changed in the last fit iteration? bool tauChanged_{false}; //! Has the mass difference parameter changed in the last fit iteration? bool deltaMChanged_{false}; //! Has the width difference parameter changed in the last fit iteration? bool deltaGammaChanged_{false}; //! Has the prompt fraction parameter changed in the last fit iteration? bool fracPromptChanged_{false}; //! Have any of the physics parameters changed in the last fit iteration? bool anythingChanged_{false}; ClassDef(LauDecayTimePhysicsModel, 0) }; template void LauDecayTimePhysicsModel::writeToJson( const std::map& models, const TString& fileName, const TString& elementName, const bool append, const int indent ) { using nlohmann::json; // Check for a signal model auto iter { models.find("signal") }; if ( iter == models.end() ) { std::cerr << "ERROR in LauDecayTimePhysicsModel::writeToJson : cannot locate the signal model in the collection, aborting" << std::endl; return; } json j; json& parametersJson = j["parameters"] = json::array(); json& signalModelJson = j["signalModel"] = json::object(); json& backgroundModelsJson = j["backgroundModels"] = json::array(); for ( const auto& [ name, model ] : models ) { if ( name == "signal" ) { signalModelJson["functionType"] = model->getFunctionType(); json& sigPars = signalModelJson["parameters"] = json::array(); for ( auto par : model->getParameters() ) { sigPars.push_back( par->name() ); parametersJson.push_back( *par ); } } else { json obj = json::object(); obj["name"] = name; json& bkgndModel = obj["model"] = json::object(); bkgndModel["functionType"] = model->getFunctionType(); json& bkgndPars = bkgndModel["parameters"] = json::array(); for ( auto par : model->getParameters() ) { bkgndPars.push_back( par->name() ); parametersJson.push_back( *par ); } backgroundModelsJson.push_back( obj ); } } const bool writeOK { LauJsonTools::writeJsonFile( j, fileName.Data(), elementName.Data(), append, indent ) }; if ( ! writeOK ) { std::cerr << "ERROR in LauDecayTimePhysicsModel::writeToJson : could not successfully write to file \"" << fileName << std::endl; } } #endif diff --git a/src/LauDecayTimePhysicsModel.cc b/src/LauDecayTimePhysicsModel.cc index 6c4c664..ae71b2a 100644 --- a/src/LauDecayTimePhysicsModel.cc +++ b/src/LauDecayTimePhysicsModel.cc @@ -1,352 +1,327 @@ /* 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 LauDecayTimePhysicsModel.cc \brief File containing implementation of LauDecayTimePhysicsModel class. */ #include #include #include #include "TSystem.h" #include "LauAbsRValue.hh" #include "LauDecayTime.hh" #include "LauDecayTimePhysicsModel.hh" ClassImp(LauDecayTimePhysicsModel); LauDecayTimePhysicsModel::LauDecayTimePhysicsModel( const LauDecayTime::FuncType type, std::vector> parameters ) : type_{type}, paramsOwned_{std::move(parameters)} { params_.resize( paramsOwned_.size() ); std::transform( paramsOwned_.begin(), paramsOwned_.end(), params_.begin(), [](std::unique_ptr& par){return par.get();} ); bool foundParams{true}; switch ( type_ ) { case LauDecayTime::FuncType::Hist : std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : Hist type needs no physics model" << std::endl; gSystem->Exit(EXIT_FAILURE); break; case LauDecayTime::FuncType::Delta : if ( ! params_.empty() ) { std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : Delta type model requires no parameters" << std::endl; gSystem->Exit(EXIT_FAILURE); } break; case LauDecayTime::FuncType::Exp : tau_ = this->findParameter("tau"); foundParams &= (tau_ != nullptr); if ( params_.size() != 1 || (!foundParams)) { std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : Exp type model requires the following parameters:\n"; std::cerr << " - the lifetime of the exponential decay: \"tau\"" << std::endl; gSystem->Exit(EXIT_FAILURE); } break; case LauDecayTime::FuncType::DeltaExp : tau_ = this->findParameter("tau"); fracPrompt_ = this->findParameter("frac_prompt"); foundParams &= (tau_ != nullptr); foundParams &= (fracPrompt_ != nullptr); if ( params_.size() != 2 || (!foundParams)) { std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : DeltaExp type model requires the following parameters:\n"; std::cerr << " - the lifetime of the exponential decay: \"tau\"\n"; std::cerr << " - the fraction of the prompt part: \"frac_prompt\"" << std::endl; gSystem->Exit(EXIT_FAILURE); } break; case LauDecayTime::FuncType::ExpTrig : tau_ = this->findParameter("tau"); deltaM_ = this->findParameter("deltaM"); foundParams &= (tau_ != nullptr); foundParams &= (deltaM_ != nullptr); if ( params_.size() != 2 || (!foundParams)) { std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : ExpTrig type model requires the following parameters:\n"; std::cerr << " - the lifetime of the exponential decay: \"tau\"\n"; std::cerr << " - the mass difference: \"deltaM\"" << std::endl; gSystem->Exit(EXIT_FAILURE); } break; case LauDecayTime::FuncType::ExpHypTrig : tau_ = this->findParameter("tau"); deltaM_ = this->findParameter("deltaM"); deltaGamma_ = this->findParameter("deltaGamma"); foundParams &= (tau_ != nullptr); foundParams &= (deltaM_ != nullptr); foundParams &= (deltaGamma_ != nullptr); if ( params_.size() != 3 || (!foundParams)) { std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : ExpHypTrig type model requires the following parameters:\n"; std::cerr << " - the lifetime of the exponential decay: \"tau\"\n"; std::cerr << " - the mass difference: \"deltaM\"\n"; std::cerr << " - the width difference: \"deltaGamma\"" << std::endl; gSystem->Exit(EXIT_FAILURE); } break; } } LauAbsRValue* LauDecayTimePhysicsModel::findParameter(const TString& parName) { for ( LauAbsRValue* par : params_ ) { if ( par && par->name().Contains(parName) ) { return par; } } std::cerr << "ERROR in LauDecayTimePhysicsModel::findParameter : Parameter \"" << parName << "\" not found." << std::endl; return nullptr; } const LauAbsRValue* LauDecayTimePhysicsModel::findParameter(const TString& parName) const { for ( const LauAbsRValue* par : params_ ) { if ( par && par->name().Contains(parName) ) { return par; } } std::cerr << "ERROR in LauDecayTimePhysicsModel::findParameter : Parameter \"" << parName << "\" not found." << std::endl; return nullptr; } void LauDecayTimePhysicsModel::initialise() { tauFloating_ = tauChanged_ = ( tau_ != nullptr ); deltaMFloating_ = deltaMChanged_ = ( deltaM_ != nullptr ); deltaGammaFloating_ = deltaGammaChanged_ = ( deltaGamma_ != nullptr ); fracPromptFloating_ = fracPromptChanged_ = ( fracPrompt_ != nullptr ); this->updateParameterCache(); tauFloating_ = ( tau_ != nullptr ) and not tau_->fixed(); deltaMFloating_ = ( deltaM_ != nullptr ) and not deltaM_->fixed(); deltaGammaFloating_ = ( deltaGamma_ != nullptr ) and not deltaGamma_->fixed(); fracPromptFloating_ = ( fracPrompt_ != nullptr ) and not fracPrompt_->fixed(); anythingFloating_ = ( tauFloating_ or deltaMFloating_ or deltaGammaFloating_ or fracPromptFloating_ ); std::cout << "INFO in LauDecayTimePhysicsModel::initialise : tau floating set to: " << (tauFloating_ ? "True" : "False") << std::endl; std::cout << "INFO in LauDecayTimePhysicsModel::initialise : deltaM floating set to: " << (deltaMFloating_ ? "True" : "False") << std::endl; std::cout << "INFO in LauDecayTimePhysicsModel::initialise : deltaGamma floating set to: " << (deltaGammaFloating_ ? "True" : "False") << std::endl; std::cout << "INFO in LauDecayTimePhysicsModel::initialise : fracPrompt floating set to: " << (fracPromptFloating_ ? "True" : "False") << std::endl; } void LauDecayTimePhysicsModel::updateParameterCache() { // Get the updated values of all parameters if ( tauChanged_ ) { tauVal_ = tau_->unblindValue(); gammaVal_ = 1.0 / tauVal_; } if ( deltaMChanged_ ) { deltaMVal_ = deltaM_->unblindValue(); } if ( deltaGammaChanged_ ) { deltaGammaVal_ = deltaGamma_->unblindValue(); } if ( fracPromptChanged_ ) { fracPromptVal_ = fracPrompt_->unblindValue(); } } void LauDecayTimePhysicsModel::propagateParUpdates() { // If none of the parameters are floating there's nothing to do if ( not anythingFloating_ ) { return; } // Otherwise, determine whether any of the floating parameters have changed and act accordingly static auto checkEquality = [](const LauAbsRValue* par, const Double_t cacheVal){return par->unblindValue() == cacheVal;}; tauChanged_ = tauFloating_ and not checkEquality(tau_, tauVal_); deltaMChanged_ = deltaMFloating_ and not checkEquality(deltaM_, deltaMVal_); deltaGammaChanged_ = deltaGammaFloating_ and not checkEquality(deltaGamma_, deltaGammaVal_); fracPromptChanged_ = fracPromptFloating_ and not checkEquality(fracPrompt_, fracPromptVal_); anythingChanged_ = ( tauChanged_ or deltaMChanged_ or deltaGammaChanged_ or fracPromptChanged_ ); if ( anythingChanged_ ) { this->updateParameterCache(); } } std::map> LauDecayTimePhysicsModel::readFromJson( const TString& fileName, const TString& elementName) { // NB deliberately not using uniform initialisation here because of this issue: // https://json.nlohmann.me/home/faq/#brace-initialization-yields-arrays const nlohmann::json j = LauJsonTools::readJsonFile( fileName.Data(), elementName.Data(), LauJsonTools::JsonType::Object ); if ( j.is_null() ) { if ( elementName != "" ) { std::cerr << "ERROR in LauDecayTimePhysicsModel::readFromJson : unable to retrieve JSON object from element \"" << elementName << "\" in file \"" << fileName << "\"" << std::endl; } else { std::cerr << "ERROR in LauDecayTimePhysicsModel::readFromJson : unable to retrieve JSON object from root element of file \"" << fileName << "\"" << std::endl; } return {}; } return readFromJson( j ); } std::map> LauDecayTimePhysicsModel::readFromJson( const nlohmann::json& j ) { using JsonType = LauJsonTools::JsonType; // the object should have the following elements: // - an array of parameter objects // - an object for the signal model settings // - an array of background model settings const std::vector mandatoryElements { std::make_pair("parameters", JsonType::Array), std::make_pair("signalModel", JsonType::Object), std::make_pair("backgroundModels", JsonType::Array) }; if ( ! LauJsonTools::checkObjectElements( j, mandatoryElements ) ) { std::cerr << "ERROR in LauDecayTimePhysicsModel::readFromJson : JSON object does not contain required elements: \"parameters\" (array), \"signalModel\" (object) and \"backgroundModels\" (array)" << std::endl; return {}; } auto parameters { LauAbsRValue::readFromJson( j.at("parameters") ) }; - auto sigModel { readSignalModelFromJson( j.at("signalModel"), parameters ) }; + auto sigModel { readPhysicsModelFromJson( j.at("signalModel"), parameters ) }; if ( ! sigModel ) { return {}; } auto backgroundModels { readBackgroundModelsFromJson( j.at("backgroundModels"), parameters ) }; if ( backgroundModels.size() != j.at("backgroundModels").size() ) { return {}; } backgroundModels.insert( std::make_pair( "signal", std::move(sigModel) ) ); return backgroundModels; } -std::unique_ptr LauDecayTimePhysicsModel::readSignalModelFromJson( const nlohmann::json& j, std::vector>& parameters ) +std::unique_ptr LauDecayTimePhysicsModel::readPhysicsModelFromJson( const nlohmann::json& j, std::vector>& parameters ) { using JsonType = LauJsonTools::JsonType; // the object should have the following elements: // - a string representing the LauDecayTime::FuncType function type // - an array of strings specifying the parameters we need const std::vector mandatoryElements { std::make_pair("functionType", JsonType::String), std::make_pair("parameters", JsonType::Array) }; if ( ! LauJsonTools::checkObjectElements( j, mandatoryElements ) ) { - std::cerr << "ERROR in LauDecayTimePhysicsModel::readSignalModelFromJson : JSON object from \"signalModel\" element does not contain required elements: \"functionType\" (string) and \"parameters\" (array)" << std::endl; + std::cerr << "ERROR in LauDecayTimePhysicsModel::readPhysicsModelFromJson : JSON object from \"signalModel\" element does not contain required elements: \"functionType\" (string) and \"parameters\" (array)" << std::endl; return {}; } auto funcType { LauJsonTools::getValue( j, "functionType" ) }; auto parNames { LauJsonTools::getValue>( j, "parameters" ) }; std::vector> params; params.reserve( parNames.size() ); for ( const auto& name : parNames ) { auto iter { std::find_if( parameters.begin(), parameters.end(), [&name]( std::unique_ptr& par ){ return par->name() == name; } ) }; if ( iter == parameters.end() ) { - const TString errMsg {"ERROR in LauDecayTimePhysicsModel::readSignalModelFromJson : Cannot locate parameter \"" + name + "\""}; + const TString errMsg {"ERROR in LauDecayTimePhysicsModel::readPhysicsModelFromJson : Cannot locate parameter \"" + name + "\""}; throw LauJsonTools::InvalidJson{errMsg.Data()}; } params.emplace_back( std::move(*iter) ); parameters.erase(iter); } return std::make_unique( funcType, std::move(params) ); } std::map> LauDecayTimePhysicsModel::readBackgroundModelsFromJson( const nlohmann::json& j, std::vector>& parameters ) { using LauJsonTools::JsonType; using LauJsonTools::checkObjectElements; using LauJsonTools::getValue; // each object in the array should have the following elements: // - a string for the name of the background component // - an object defining the model, which itself has elements: // = a string representing the LauDecayTime::FuncType function type // = an array of strings specifying the parameters we need std::vector mandatoryElements { std::make_pair("name", JsonType::String), std::make_pair("model", JsonType::Object) }; for ( const auto& comp : j ) { if ( ! checkObjectElements( comp, mandatoryElements ) ) { std::cerr << "ERROR in LauDecayTimePhysicsModel::readBackgroundModelsFromJson : JSON object in \"backgroundModels\" array does not contain required elements: \"name\" (string) and \"model\" (object)" << std::endl; return {}; } } mandatoryElements = { std::make_pair("functionType", JsonType::String), std::make_pair("parameters", JsonType::Array) }; std::map> models; for ( const auto& comp : j ) { auto name { getValue( comp, "name" ) }; - - const auto& model { comp.at("model") }; - if ( ! checkObjectElements( model, mandatoryElements ) ) { - std::cerr << "ERROR in LauDecayTimePhysicsModel::readBackgroundModelsFromJson : JSON object defining model in \"backgroundModels\" array does not contain required elements: \"functionType\" (string) and \"parameters\" (array)" << std::endl; - return {}; - } - - auto funcType { getValue( model, "functionType" ) }; - auto parNames { getValue>( model, "parameters" ) }; - - std::vector> params; - params.reserve( parNames.size() ); - - for ( const auto& parName : parNames ) { - - auto iter { std::find_if( parameters.begin(), parameters.end(), [&parName]( std::unique_ptr& par ){ return par->name() == parName; } ) }; - if ( iter == parameters.end() ) { - const TString errMsg {"ERROR in LauDecayTimePhysicsModel::readBackgroundModelsFromJson : Cannot locate parameter \"" + parName + "\""}; - throw LauJsonTools::InvalidJson{errMsg.Data()}; - } - - params.emplace_back( std::move(*iter) ); - parameters.erase(iter); - } - - models[name] = std::make_unique( funcType, std::move(params) ); + models[name] = readPhysicsModelFromJson( comp.at("model"), parameters ); } return models; }