diff --git a/examples/GenFit3K.cc b/examples/GenFit3K.cc --- a/examples/GenFit3K.cc +++ b/examples/GenFit3K.cc @@ -138,10 +138,14 @@ res = sigModel->addResonance("NonReson", 0, LauAbsResonance::FlatNR); - // Reset the maximum signal DP ASq value - // This will be automatically adjusted to avoid bias or extreme - // inefficiency if you get the value wrong but best to set this by - // hand once you've found the right value through some trial and error. + // Set the maximum signal DP ASq value + // If you do not provide a value, one will be determined automatically, + // which should be close to the true maximum but is not guaranteed to + // be optimal. + // Any value, whether manually provided or automatically determined, + // will be automatically adjusted to avoid bias or extreme inefficiency + // but it is best to set this by hand once you've found the right value + // through some trial and error. sigModel->setASqMaxValue(15.8); // Create the fit model diff --git a/examples/GenFit3KS.cc b/examples/GenFit3KS.cc --- a/examples/GenFit3KS.cc +++ b/examples/GenFit3KS.cc @@ -115,10 +115,14 @@ /*reson =*/ sigModel->addResonance("f_2(2010)", 3, LauAbsResonance::RelBW); /*reson =*/ sigModel->addResonance("chi_c0", 3, LauAbsResonance::RelBW); - // Reset the maximum signal DP ASq value - // This will be automatically adjusted to avoid bias or extreme - // inefficiency if you get the value wrong but best to set this by - // hand once you've found the right value through some trial and error. + // Set the maximum signal DP ASq value + // If you do not provide a value, one will be determined automatically, + // which should be close to the true maximum but is not guaranteed to + // be optimal. + // Any value, whether manually provided or automatically determined, + // will be automatically adjusted to avoid bias or extreme inefficiency + // but it is best to set this by hand once you've found the right value + // through some trial and error. sigModel->setASqMaxValue(0.285); // Create the fit model diff --git a/examples/GenFit3pi.cc b/examples/GenFit3pi.cc --- a/examples/GenFit3pi.cc +++ b/examples/GenFit3pi.cc @@ -150,10 +150,14 @@ reson = sigModel->addResonance(nrName, 0, LauAbsResonance::BelleSymNRNoInter); reson->setResonanceParameter("alpha", 0.2); - // Reset the maximum signal DP ASq value - // This will be automatically adjusted to avoid bias or extreme - // inefficiency if you get the value wrong but best to set this by - // hand once you've found the right value through some trial and error. + // Set the maximum signal DP ASq value + // If you do not provide a value, one will be determined automatically, + // which should be close to the true maximum but is not guaranteed to + // be optimal. + // Any value, whether manually provided or automatically determined, + // will be automatically adjusted to avoid bias or extreme inefficiency + // but it is best to set this by hand once you've found the right value + // through some trial and error. sigModel->setASqMaxValue(0.35); // Create the fit model diff --git a/examples/GenFit3pi.py.in b/examples/GenFit3pi.py.in --- a/examples/GenFit3pi.py.in +++ b/examples/GenFit3pi.py.in @@ -131,10 +131,14 @@ reson = sigModel.addResonance('BelleNR_Swave', 0, ROOT.LauAbsResonance.BelleSymNRNoInter) reson.setResonanceParameter("alpha", 0.2) -# Reset the maximum signal DP ASq value -# This will be automatically adjusted to avoid bias or extreme -# inefficiency if you get the value wrong but best to set this by -# hand once you've found the right value through some trial and error. +# Set the maximum signal DP ASq value +# If you do not provide a value, one will be determined automatically, +# which should be close to the true maximum but is not guaranteed to +# be optimal. +# Any value, whether manually provided or automatically determined, +# will be automatically adjusted to avoid bias or extreme inefficiency +# but it is best to set this by hand once you've found the right value +# through some trial and error. sigModel.setASqMaxValue(0.35) # Create the fit model diff --git a/examples/GenFitBelleCPKpipi.cc b/examples/GenFitBelleCPKpipi.cc --- a/examples/GenFitBelleCPKpipi.cc +++ b/examples/GenFitBelleCPKpipi.cc @@ -130,12 +130,15 @@ negSigModel->setIntFileName("integ_neg.dat"); posSigModel->setIntFileName("integ_pos.dat"); - // Reset the maximum signal DP ASq value - // This will be automatically adjusted to avoid bias or extreme - // inefficiency if you get the value wrong but best to set this by - // hand once you've found the right value through some trial and - // error. - Double_t aSqMaxValue = 1.62; + // Set the maximum signal DP ASq value + // If you do not provide a value, one will be determined automatically, + // which should be close to the true maximum but is not guaranteed to + // be optimal. + // Any value, whether manually provided or automatically determined, + // will be automatically adjusted to avoid bias or extreme inefficiency + // but it is best to set this by hand once you've found the right value + // through some trial and error. + const Double_t aSqMaxValue{1.62}; negSigModel->setASqMaxValue(aSqMaxValue); posSigModel->setASqMaxValue(aSqMaxValue); diff --git a/examples/GenFitDs2KKpi.cc b/examples/GenFitDs2KKpi.cc --- a/examples/GenFitDs2KKpi.cc +++ b/examples/GenFitDs2KKpi.cc @@ -111,10 +111,14 @@ /*reson =*/ sigModel->addResonance("K*0(892)", 2, LauAbsResonance::RelBW); /*reson =*/ sigModel->addResonance("NonReson", 0, LauAbsResonance::FlatNR); - // Reset the maximum signal DP ASq value - // This will be automatically adjusted to avoid bias or extreme - // inefficiency if you get the value wrong but best to set this by - // hand once you've found the right value through some trial and error. + // Set the maximum signal DP ASq value + // If you do not provide a value, one will be determined automatically, + // which should be close to the true maximum but is not guaranteed to + // be optimal. + // Any value, whether manually provided or automatically determined, + // will be automatically adjusted to avoid bias or extreme inefficiency + // but it is best to set this by hand once you've found the right value + // through some trial and error. sigModel->setASqMaxValue(310.0); // Create the fit model diff --git a/examples/GenFitEFKLLM.cc b/examples/GenFitEFKLLM.cc --- a/examples/GenFitEFKLLM.cc +++ b/examples/GenFitEFKLLM.cc @@ -118,10 +118,14 @@ res = sigModel->addResonance("K*0_0(1430)", 1, LauAbsResonance::EFKLLM); res->setResonanceParameter("massFactor", -2.0); - // Reset the maximum signal DP ASq value - // This will be automatically adjusted to avoid bias or extreme - // inefficiency if you get the value wrong but best to set this by - // hand once you've found the right value through some trial and error. + // Set the maximum signal DP ASq value + // If you do not provide a value, one will be determined automatically, + // which should be close to the true maximum but is not guaranteed to + // be optimal. + // Any value, whether manually provided or automatically determined, + // will be automatically adjusted to avoid bias or extreme inefficiency + // but it is best to set this by hand once you've found the right value + // through some trial and error. sigModel->setASqMaxValue(14.5); // Create the fit model diff --git a/examples/SimFitTask.cc b/examples/SimFitTask.cc --- a/examples/SimFitTask.cc +++ b/examples/SimFitTask.cc @@ -137,11 +137,14 @@ res->fixMass(kFALSE); res = sigModel->addResonance("K*0_0(1430)", 1, LauAbsResonance::LASS); - // Reset the maximum signal DP ASq value - // This will be automatically adjusted to avoid bias or extreme - // inefficiency if you get the value wrong but best to set this by - // hand once you've found the right value through some trial and - // error. + // Set the maximum signal DP ASq value + // If you do not provide a value, one will be determined automatically, + // which should be close to the true maximum but is not guaranteed to + // be optimal. + // Any value, whether manually provided or automatically determined, + // will be automatically adjusted to avoid bias or extreme inefficiency + // but it is best to set this by hand once you've found the right value + // through some trial and error. sigModel->setASqMaxValue(1.25); TString integralsFileName("integ-"); diff --git a/inc/LauASqMaxFinder.hh b/inc/LauASqMaxFinder.hh new file mode 100644 --- /dev/null +++ b/inc/LauASqMaxFinder.hh @@ -0,0 +1,94 @@ + +/* +Copyright 2022 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 LauASqMaxFinder.hh + \brief File containing declaration of LauASqMaxFinder class. +*/ + +#ifndef LAU_ASQ_MAX_FINDER +#define LAU_ASQ_MAX_FINDER + +#include "LauFitObject.hh" + +#include +#include + +class LauIsobarDynamics; +class LauParameter; + +/*! \class LauASqMaxFinder + \brief Class for locating the |A|^2 max for a given model + + Currently supports finding the |A|^2 max for a LauIsobarDynamics object +*/ + +class LauASqMaxFinder : public LauFitObject { + + public: + //! Constructor + /*! + \param [in] iso the isobar dynamics for which we need to find the maximum |A|^2 value + */ + explicit LauASqMaxFinder(LauIsobarDynamics& iso); + + //! Run the minimisation to locate the maximum |A|^2 value + /*! + \return the maximum |A|^2 value + */ + Double_t find(); + + //! This function sets the parameter values from Minuit + /*! + This function has to be public since it is called from the global FCN. + It should not be called otherwise! + + \param [in] par an array storing the various parameter values + \param [in] npar the number of free parameters + */ + void setParsFromMinuit(Double_t* par, Int_t npar) override; + + //! Calculate the new value of the negative log likelihood + /*! + This function has to be public since it is called from the global FCN. + It should not be called otherwise! + */ + Double_t getTotNegLogLikelihood() override; + + private: + //! The isobar dynamics for which we need to find the maximum |A|^2 value + LauIsobarDynamics& iso_; + + //! Flag to control printing of warnings if the fit wanders outside the DP + Bool_t printMinimisationWarnings_{kTRUE}; + + //! The fit parameters (owned by us) + std::vector> ownedParams_; + + //! The fit parameters in a form to be passed to the minimiser + std::vector params_; + + ClassDefOverride(LauASqMaxFinder,0) +}; + +#endif diff --git a/inc/LauAbsFitter.hh b/inc/LauAbsFitter.hh --- a/inc/LauAbsFitter.hh +++ b/inc/LauAbsFitter.hh @@ -36,11 +36,11 @@ #ifndef LAU_ABS_FITTER #define LAU_ABS_FITTER -#include - #include "Rtypes.h" #include "TMatrixDfwd.h" +#include + class LauFitObject; class LauParameter; @@ -62,7 +62,7 @@ }; //! Destructor - virtual ~LauAbsFitter() {} + virtual ~LauAbsFitter() = default; //! Initialise the fitter, setting the information on the parameters /*! @@ -72,7 +72,7 @@ virtual void initialise( LauFitObject* fitObj, const std::vector& parameters ) = 0; //! Get the object that controls the calculation of the likelihood - virtual LauFitObject* getFitObject() =0; + virtual LauFitObject* getFitObject() = 0; //! Get the total number of fit parameters virtual UInt_t nParameters() const = 0; @@ -125,17 +125,22 @@ protected: //! Constructor - LauAbsFitter() {} + LauAbsFitter() = default; private: //! Copy constructor - private and not implemented - LauAbsFitter( const LauAbsFitter& ); + LauAbsFitter( const LauAbsFitter& ) = delete; + + //! Move constructor - private and not implemented + LauAbsFitter( LauAbsFitter&& ) = delete; //! Copy assignment operator - private and not implemented - LauAbsFitter& operator=( const LauAbsFitter& rhs ); + LauAbsFitter& operator=( const LauAbsFitter& rhs ) = delete; + + //! Move assignment operator - private and not implemented + LauAbsFitter& operator=( LauAbsFitter&& rhs ) = delete; ClassDef(LauAbsFitter,0); }; #endif - diff --git a/inc/LauFitObject.hh b/inc/LauFitObject.hh --- a/inc/LauFitObject.hh +++ b/inc/LauFitObject.hh @@ -26,10 +26,6 @@ \brief File containing declaration of LauFitObject class. */ -/*! \class LauFitObject - \brief The abstract interface for the objects that control the calculation of the likelihood. -*/ - #ifndef LAU_FIT_OBJECT #define LAU_FIT_OBJECT @@ -41,11 +37,15 @@ #include "LauAbsFitter.hh" +/*! \class LauFitObject + \brief The abstract interface for the objects that control the calculation of the likelihood. +*/ + class LauFitObject : public TObject { public: //! Destructor - virtual ~LauFitObject() {} + virtual ~LauFitObject() = default; //! Turn on or off the computation of asymmetric errors (e.g. MINOS routine in Minuit) /*! diff --git a/inc/LauFitter.hh b/inc/LauFitter.hh --- a/inc/LauFitter.hh +++ b/inc/LauFitter.hh @@ -26,25 +26,29 @@ \brief File containing declaration of LauFitter class. */ -/*! \class LauFitter - \brief Factory class for creating and providing access to the fitter. - - The fitter type can be set before first access to determine which fitter is used. -*/ - #ifndef LAU_FITTER #define LAU_FITTER +#include "LauPrint.hh" + #include "Rtypes.h" #include "TString.h" +#include + class LauAbsFitter; -class LauFitter { +/*! \class LauFitter + \brief Factory class for creating and providing access to the fitter. + + The fitter type and verbosity can be set before first access to determine which fitter is used. +*/ + +class LauFitter final { public: //! The types of fitter available - enum Type { + enum class Type { Minuit /*!< the Minuit fitter */ }; @@ -52,33 +56,63 @@ /*! \param [in] type the type of the fitter (default set to Minuit) */ - static void setFitterType( Type type ); + static void setFitterType( const Type type ); + + //! Set the verbosity level of the fitter + /*! + \param [in] level the level of verbosity of the fitter (default set to Standard) + */ + static void setFitterVerbosity( const LauOutputLevel level ); + + //! Set the maximum number of parameters for the fitter + /*! + \param [in] maxPars the maximum number of parameters for the fitter (default set to 100) + */ + static void setFitterMaxPars( const UInt_t maxPars ); //! Method that provides access to the singleton fitter /*! - \return a pointer to a singleton LauAbsFitter object + \return a reference to a singleton LauAbsFitter object + */ + static LauAbsFitter& fitter(); + + //! Destroy the current fitter + /*! + A new fitter will be created on the next call to LauFitter::fitter */ - static LauAbsFitter* fitter(); + static void destroyFitter(); private: //! Constructor - LauFitter() {} + LauFitter() = default; //! Destructor - virtual ~LauFitter() {} + ~LauFitter() = default; + + //! Copy constructor (deleted) + LauFitter( const LauFitter& ) = delete; - //! Copy constructor (not implemented) - LauFitter( const LauFitter& ); + //! Move constructor (deleted) + LauFitter( LauFitter&& ) = delete; - //! Copy assignment operator (not implemented) - LauFitter& operator=( const LauFitter& ); + //! Copy assignment operator (deleted) + LauFitter& operator=( const LauFitter& ) = delete; + + //! Move assignment operator (deleted) + LauFitter& operator=( LauFitter&& ) = delete; //! Pointer to the singleton fitter instance - static LauAbsFitter* theInstance_; + static std::unique_ptr theInstance_; //! The fitter type static Type fitterType_; + //! The fitter verbosity + static LauOutputLevel fitterVerbosity_; + + //! The maximum number of parameters for the fitter + static UInt_t fitterMaxPars_; + ClassDef(LauFitter,0); }; diff --git a/inc/LauIsobarDynamics.hh b/inc/LauIsobarDynamics.hh --- a/inc/LauIsobarDynamics.hh +++ b/inc/LauIsobarDynamics.hh @@ -208,9 +208,11 @@ //! Set the maximum value of A squared to be used in the accept/reject /*! + Disables the automatic determination of ASqMax + \param [in] value the new value */ - inline void setASqMaxValue(Double_t value) {aSqMaxSet_ = value;} + inline void setASqMaxValue(Double_t value) {aSqMaxSet_ = value; aSqMaxAuto_ = kFALSE;} //! Retrieve the maximum value of A squared to be used in the accept/reject /*! @@ -860,6 +862,7 @@ //! Force the symmetrisation of the integration in m13 <-> m23 for non-symmetric but flavour-conjugate final states Bool_t forceSymmetriseIntegration_; + //! The storage of the integration scheme std::vector dpPartialIntegralInfo_; @@ -956,6 +959,9 @@ //! The maximum value of A squared that has been seen so far while generating Double_t aSqMaxVar_; + //! Flag to generate aSqMaxSet_ once generate is called + Bool_t aSqMaxAuto_{kTRUE}; + //! The helicity flip flag for new amplitude components Bool_t flipHelicity_; diff --git a/inc/LauMinuit.hh b/inc/LauMinuit.hh --- a/inc/LauMinuit.hh +++ b/inc/LauMinuit.hh @@ -35,13 +35,14 @@ #ifndef LAU_MINUIT #define LAU_MINUIT -#include -#include +#include "LauAbsFitter.hh" +#include "LauPrint.hh" #include "Rtypes.h" #include "TMatrixD.h" -#include "LauAbsFitter.hh" +#include +#include class LauParameter; class TVirtualFitter; @@ -51,7 +52,7 @@ public: //! Destructor - virtual ~LauMinuit(); + virtual ~LauMinuit() = default; //! Initialise the fitter, setting the information on the parameters /*! @@ -118,37 +119,46 @@ friend class LauFitter; //! Constructor - explicit LauMinuit( Int_t maxPar = 100 ); + LauMinuit( const UInt_t maxPar = 100, const LauOutputLevel verbosity = LauOutputLevel::Standard ); //! Copy constructor - private and not implemented - LauMinuit( const LauMinuit& ); + LauMinuit( const LauMinuit& ) = delete; + + //! Move constructor - private and not implemented + LauMinuit( LauMinuit&& ) = delete; //! Copy assignment operator - private and not implemented - LauMinuit& operator=( const LauMinuit& rhs ); + LauMinuit& operator=( const LauMinuit& ) = delete; + + //! Move assignment operator - private and not implemented + LauMinuit& operator=( LauMinuit&& ) = delete; //! The interface to Minuit - TVirtualFitter* minuit_; + TVirtualFitter* minuit_{nullptr}; //! The maximum number of parameters - const UInt_t maxPar_; + const UInt_t maxPar_{100}; + + //! The verbosity level of the fitter + const LauOutputLevel outputLevel_{LauOutputLevel::Standard}; //! The fit parameters std::vector params_; //! The total number of parameters - UInt_t nParams_; + UInt_t nParams_{0}; //! The number of free parameters - UInt_t nFreeParams_; + UInt_t nFreeParams_{0}; //! Option to perform a two stage fit - Bool_t twoStageFit_; + Bool_t twoStageFit_{kFALSE}; //! Option to use asymmetric errors - Bool_t useAsymmFitErrors_; + Bool_t useAsymmFitErrors_{kFALSE}; //! The status of the fit - FitStatus fitStatus_; + FitStatus fitStatus_{-1,0.0,0.0}; //! The covariance matrix TMatrixD covMatrix_; @@ -157,4 +167,3 @@ }; #endif - diff --git a/inc/LauPrint.hh b/inc/LauPrint.hh --- a/inc/LauPrint.hh +++ b/inc/LauPrint.hh @@ -23,13 +23,7 @@ */ /*! \file LauPrint.hh - \brief File containing declaration of LauPrint class. -*/ - -/*! \class LauPrint - \brief Class to define various output print commands. - - Class to define various output print commands (for output tables etc..) + \brief File containing declaration of LauPrint class and LauOutputLevel enum. */ #ifndef LAU_PRINT @@ -40,6 +34,24 @@ #include +/*! \enum LauOutputLevel + \brief Enumeration to define verbosity level for various printouts +*/ + +enum class LauOutputLevel { + None = -1, //!< Zero printout + Quiet = 0, //!< Reduced printout + Standard = 1, //!< Normal level of printout + Verbose = 2, //!< Verbose printout + ExtraVerbose = 3 //!< Highly verbose printout +}; + +/*! \class LauPrint + \brief Class to define various output print commands. + + Class to define various output print commands (for output tables etc..) +*/ + class LauPrint { public: diff --git a/inc/Laura++_LinkDef.h b/inc/Laura++_LinkDef.h --- a/inc/Laura++_LinkDef.h +++ b/inc/Laura++_LinkDef.h @@ -51,6 +51,7 @@ #pragma link C++ class LauAbsResonance+; #pragma link C++ class LauAbsRValue+; #pragma link C++ class LauArgusPdf+; +#pragma link C++ class LauASqMaxFinder+; #pragma link C++ class LauAsymmCalc+; #pragma link C++ class LauBelleCPCoeffSet+; #pragma link C++ class LauBelleNR+; diff --git a/src/LauASqMaxFinder.cc b/src/LauASqMaxFinder.cc new file mode 100644 --- /dev/null +++ b/src/LauASqMaxFinder.cc @@ -0,0 +1,172 @@ + +/* +Copyright 2022 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 LauASqMaxFinder.cc + \brief File containing implementation of LauASqMaxFinder class. +*/ + +#include "LauASqMaxFinder.hh" + +#include "LauAbsResonance.hh" +#include "LauFitter.hh" +#include "LauIsobarDynamics.hh" +#include "LauKinematics.hh" +#include "LauParameter.hh" + +#include "TSystem.h" + +#include +#include +#include + +ClassImp( LauASqMaxFinder ) + +LauASqMaxFinder::LauASqMaxFinder(LauIsobarDynamics& iso) : + iso_{iso} +{ + const LauKinematics* kinematics { iso_.getKinematics() }; + const Double_t m13SqMin{kinematics->getm13SqMin()}; + const Double_t m13SqMax{kinematics->getm13SqMax()}; + const Double_t m23SqMin{kinematics->getm23SqMin()}; + const Double_t m23SqMax{kinematics->getm23SqMax()}; + + ownedParams_.resize(2); + ownedParams_[0] = std::make_unique("m13Sq", 0.5*(m13SqMin+m13SqMax), m13SqMin, m13SqMax, kFALSE ); + ownedParams_[1] = std::make_unique("m23Sq", 0.5*(m23SqMin+m23SqMax), m23SqMin, m23SqMax, kFALSE ); + + params_.resize(2); + params_[0] = ownedParams_[0].get(); + params_[1] = ownedParams_[1].get(); +} + +void LauASqMaxFinder::setParsFromMinuit(Double_t* par, Int_t npar) +{ + if ( static_cast(npar) != params_.size() ) { + std::cerr << "ERROR in LauASqMaxFinder::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n"; + std::cerr << " Expected: " << params_.size() << ".\n" << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + + for ( Int_t i{0}; i < npar; ++i ) { + params_[i]->value( par[i] ); + } + + return; +} + +Double_t LauASqMaxFinder::getTotNegLogLikelihood() +{ + const Double_t m13Sq { params_[0]->unblindValue() }; + const Double_t m23Sq { params_[1]->unblindValue() }; + + if ( ! iso_.getKinematics()->withinDPLimits(m13Sq, m23Sq) ) { + if ( printMinimisationWarnings_ ) { + std::cerr << "WARNING in LauASqMaxFinder::getTotNegLogLikelihood : Fit has wandered outside the phase space boundary\n"; + std::cerr << " : Returning 0.0 value (further warnings of this type will be suppressed)" << std::endl; + printMinimisationWarnings_ = kFALSE; + } + return 0.0; + } + + iso_.calcLikelihoodInfo( m13Sq, m23Sq ); + return -iso_.getEvtIntensity(); +} + +Double_t LauASqMaxFinder::find() +{ + Double_t ASqMax{-1.0}; + + const LauKinematics* kinematics { iso_.getKinematics() }; + LauKinematics localKin { + kinematics->getm1(), + kinematics->getm2(), + kinematics->getm3(), + kinematics->getmParent() + }; + + const Double_t m13SqMin { localKin.getm13Min() }; + const Double_t m13SqMax { localKin.getm13Max() }; + const Double_t m23SqMin { localKin.getm23Min() }; + const Double_t m23SqMax { localKin.getm23Max() }; + + LauFitter::setFitterMaxPars( 2 ); + LauFitter::setFitterVerbosity( LauOutputLevel::None ); + + for ( UInt_t iRes{0}; iRes < iso_.getnTotAmp(); ++iRes ) { + // Find the center of the resonance and set that as the start of the search + const LauAbsResonance* res { const_cast(iso_).getResonance(iRes) }; + const Int_t bachelor { res->getPairInt() }; + const Double_t mass { res->getMass() }; + const Double_t width { res->getWidth() }; + + const std::array massVals { mass - width, mass, mass + width }; + const std::array cosHelVals { -1.0, 0.0, 1.0 }; + for ( Double_t massVal : massVals ) { + for ( Double_t cosHel : cosHelVals ) { + + switch ( bachelor ) { + case 0 : + localKin.updateKinematics( 0.5*(m13SqMax+m13SqMin), 0.5*(m23SqMax+m23SqMin) ); + break; + case 1 : + localKin.updateKinematicsFrom23( massVal, cosHel ); + break; + case 2 : + localKin.updateKinematicsFrom13( massVal, cosHel ); + break; + case 3 : + localKin.updateKinematicsFrom12( massVal, cosHel ); + break; + + default : + break; + } + + params_[0]->initValue( localKin.getm13Sq() ); + params_[1]->initValue( localKin.getm23Sq() ); + + params_[0]->error(0.1); + params_[1]->error(0.1); + + LauFitter::fitter().initialise( this, params_ ); + + const LauAbsFitter::FitStatus fitResult { LauFitter::fitter().minimise() }; + + if ( fitResult.status != 0 and -fitResult.NLL > ASqMax ) { + ASqMax = -fitResult.NLL; + } + + } + } + } + + LauFitter::destroyFitter(); + + if ( ASqMax < 0.0 ) { + std::cerr << "WARNING in LauASqMaxFinder::find : Could not find ASqMax\n"; + std::cerr << " : Returning -1.0" << std::endl; + } + + return ASqMax; +} diff --git a/src/LauAbsFitModel.cc b/src/LauAbsFitModel.cc --- a/src/LauAbsFitModel.cc +++ b/src/LauAbsFitModel.cc @@ -566,36 +566,36 @@ this->checkInitFitParams(); // Initialise the fitter - LauFitter::fitter()->useAsymmFitErrors( this->useAsymmFitErrors() ); - LauFitter::fitter()->twoStageFit( this->twoStageFit() ); - LauFitter::fitter()->initialise( this, fitVars_ ); + LauFitter::fitter().useAsymmFitErrors( this->useAsymmFitErrors() ); + LauFitter::fitter().twoStageFit( this->twoStageFit() ); + LauFitter::fitter().initialise( this, fitVars_ ); - this->startNewFit( LauFitter::fitter()->nParameters(), LauFitter::fitter()->nFreeParameters() ); + this->startNewFit( LauFitter::fitter().nParameters(), LauFitter::fitter().nFreeParameters() ); // Now ready for minimisation step std::cout << "\nINFO in LauAbsFitModel::fitExpt : Start minimisation...\n"; - LauAbsFitter::FitStatus fitResult = LauFitter::fitter()->minimise(); + LauAbsFitter::FitStatus fitResult = LauFitter::fitter().minimise(); // If we're doing a two stage fit we can now release (i.e. float) // the 2nd stage parameters and re-fit if (this->twoStageFit()) { if ( fitResult.status != 3 ) { std::cerr << "WARNING in LauAbsFitModel:fitExpt : Not running second stage fit since first stage failed." << std::endl; - LauFitter::fitter()->releaseSecondStageParameters(); + LauFitter::fitter().releaseSecondStageParameters(); } else { - LauFitter::fitter()->releaseSecondStageParameters(); - this->startNewFit( LauFitter::fitter()->nParameters(), LauFitter::fitter()->nFreeParameters() ); - fitResult = LauFitter::fitter()->minimise(); + LauFitter::fitter().releaseSecondStageParameters(); + this->startNewFit( LauFitter::fitter().nParameters(), LauFitter::fitter().nFreeParameters() ); + fitResult = LauFitter::fitter().minimise(); } } - const TMatrixD& covMat = LauFitter::fitter()->covarianceMatrix(); + const TMatrixD& covMat = LauFitter::fitter().covarianceMatrix(); this->storeFitStatus( fitResult, covMat ); // Store the final fit results and errors into protected internal vectors that // all sub-classes can use within their own finalFitResults implementation // used below (e.g. putting them into an ntuple in a root file) - LauFitter::fitter()->updateParameters(); + LauFitter::fitter().updateParameters(); } void LauAbsFitModel::calculateSPlotData() diff --git a/src/LauFitter.cc b/src/LauFitter.cc --- a/src/LauFitter.cc +++ b/src/LauFitter.cc @@ -26,21 +26,25 @@ \brief File containing implementation of LauFitter methods. */ -#include - #include "LauFitter.hh" -#include "LauAbsFitter.hh" + #include "LauMinuit.hh" +#include +#include +#include + -LauAbsFitter* LauFitter::theInstance_ = 0; -LauFitter::Type LauFitter::fitterType_ = LauFitter::Minuit; +std::unique_ptr LauFitter::theInstance_ = nullptr; +LauFitter::Type LauFitter::fitterType_ = LauFitter::Type::Minuit; +LauOutputLevel LauFitter::fitterVerbosity_ = LauOutputLevel::Standard; +UInt_t LauFitter::fitterMaxPars_ = 100; ClassImp(LauFitter) -void LauFitter::setFitterType( Type type ) +void LauFitter::setFitterType( const Type type ) { - if ( theInstance_ != 0 ) { + if ( theInstance_ != nullptr ) { std::cerr << "ERROR in LauFitter::setFitterType : The fitter has already been created, cannot change the type now." << std::endl; return; } @@ -48,17 +52,48 @@ fitterType_ = type; } -LauAbsFitter* LauFitter::fitter() +void LauFitter::setFitterVerbosity( const LauOutputLevel level ) +{ + if ( theInstance_ != nullptr ) { + std::cerr << "ERROR in LauFitter::setFitterVerbosity : The fitter has already been created, cannot change the verbosity now." << std::endl; + return; + } + + fitterVerbosity_ = level; +} + +void LauFitter::setFitterMaxPars( const UInt_t maxPars ) +{ + if ( theInstance_ != nullptr ) { + std::cerr << "ERROR in LauFitter::setFitterMaxPars : The fitter has already been created, cannot change the maximum number of parameters now." << std::endl; + return; + } + + fitterMaxPars_ = maxPars; +} + +LauAbsFitter& LauFitter::fitter() { - // Returns a pointer to a singleton LauAbsFitter object. + // Returns a reference to a singleton LauAbsFitter object. // Creates the object the first time it is called. - if ( theInstance_ == 0 ) { - if ( fitterType_ == Minuit ) { - theInstance_ = new LauMinuit(); + if ( theInstance_ == nullptr ) { + if ( fitterType_ == Type::Minuit ) { + // NB cannot use std::make_unique here since the LauMinuit constructor is private + theInstance_.reset( new LauMinuit( fitterMaxPars_, fitterVerbosity_ ) ); } } - return theInstance_; + return *theInstance_; } +void LauFitter::destroyFitter() +{ + // destroy the current fitter + theInstance_.reset(); + + // restore the default settings + fitterType_ = LauFitter::Type::Minuit; + fitterVerbosity_ = LauOutputLevel::Standard; + fitterMaxPars_ = 100; +} diff --git a/src/LauIsobarDynamics.cc b/src/LauIsobarDynamics.cc --- a/src/LauIsobarDynamics.cc +++ b/src/LauIsobarDynamics.cc @@ -58,6 +58,7 @@ #include "LauResonanceInfo.hh" #include "LauResonanceMaker.hh" #include "LauRhoOmegaMix.hh" +#include "LauASqMaxFinder.hh" ClassImp(LauIsobarDynamics) @@ -2287,10 +2288,22 @@ // We need to make sure to calculate everything for every resonance integralsToBeCalculated_.clear(); - for ( UInt_t i(0); i < nAmp_+nIncohAmp_; ++i ) { + for ( UInt_t i{0}; i < nAmp_+nIncohAmp_; ++i ) { integralsToBeCalculated_.insert(i); } + if ( aSqMaxAuto_ ) { + std::cout<<"INFO in LauIsobarDynamics::generate : Starting auto-location of |A|^2 maximum"< "< -#include - -#include "TMatrixD.h" -#include "TVirtualFitter.h" +#include "LauMinuit.hh" #include "LauFitObject.hh" #include "LauFitter.hh" -#include "LauMinuit.hh" #include "LauParameter.hh" #include "LauParamFixed.hh" +#include "TMatrixD.h" +#include "TVirtualFitter.h" + +#include +#include +#include + // It's necessary to define an external function that specifies the address of the function // that Minuit needs to minimise. Minuit doesn't know about any classes - therefore // use gMinuit->SetFCN(external_function), gMinuit->SetObjectFit(this). @@ -49,31 +51,31 @@ ClassImp(LauMinuit) -LauMinuit::LauMinuit( Int_t maxPar ) : LauAbsFitter(), - minuit_(0), - maxPar_(maxPar), - nParams_(0), - nFreeParams_(0), - twoStageFit_(kFALSE), - useAsymmFitErrors_(kFALSE), - fitStatus_({-1,0.0,0.0}) +LauMinuit::LauMinuit( const UInt_t maxPar, const LauOutputLevel verbosity ) : LauAbsFitter(), + maxPar_{maxPar}, + outputLevel_{verbosity} { TVirtualFitter::SetDefaultFitter( "Minuit" ); - minuit_ = TVirtualFitter::Fitter( 0, maxPar_ ); -} + minuit_ = TVirtualFitter::Fitter( nullptr, maxPar_ ); -LauMinuit::~LauMinuit() -{ + // Set the printout level + std::array argL { static_cast(outputLevel_) }; + minuit_->ExecuteCommand("SET PRINT", argL.data(), argL.size()); + if ( outputLevel_ == LauOutputLevel::None ) { + minuit_->ExecuteCommand("SET NOWARNINGS", argL.data(), 0); + } } void LauMinuit::initialise( LauFitObject* fitObj, const std::vector& parameters ) { // Check whether we're going to use asymmetric errors - if (useAsymmFitErrors_ == kTRUE) { - std::cout << "INFO in LauMinuit::fit : We are going to calculate the asymmetric fit errors." << std::endl; - std::cout << " : This will, in general, significantly increase the CPU time required for fitting." << std::endl; - } else { - std::cout << "INFO in LauMinuit::fit : Only parabolic errors will be calculated." << std::endl; + if ( outputLevel_ > LauOutputLevel::Quiet ) { + if (useAsymmFitErrors_ == kTRUE) { + std::cout << "INFO in LauMinuit::fit : We are going to calculate the asymmetric fit errors." << std::endl; + std::cout << " : This will, in general, significantly increase the CPU time required for fitting." << std::endl; + } else { + std::cout << "INFO in LauMinuit::fit : Only parabolic errors will be calculated." << std::endl; + } } // Store the parameters @@ -87,8 +89,10 @@ minuit_->Clear(); nParams_ = params_.size(); - std::cout << "INFO in LauMinuit::initialise : Setting fit parameters" << std::endl; - std::cout << " : Total number of parameters = " << nParams_ << std::endl; + if ( outputLevel_ > LauOutputLevel::Quiet ) { + std::cout << "INFO in LauMinuit::initialise : Setting fit parameters" << std::endl; + std::cout << " : Total number of parameters = " << nParams_ << std::endl; + } // Define the default relative error const Double_t defaultError(0.01); @@ -114,12 +118,16 @@ } Bool_t fixVar = params_[i]->fixed(); - std::cout << " : Setting parameter " << i << " called " << name << " to have initial value " << initVal << ", error " << initErr << " and range " << minVal << " to " << maxVal << std::endl; + if ( outputLevel_ > LauOutputLevel::Quiet ) { + std::cout << " : Setting parameter " << i << " called " << name << " to have initial value " << initVal << ", error " << initErr << " and range " << minVal << " to " << maxVal << std::endl; + } minuit_->SetParameter(i, name, initVal, initErr, minVal, maxVal); // Fix parameter if required if (fixVar == kTRUE) { - std::cout << " : Fixing parameter " << i << std::endl; + if ( outputLevel_ > LauOutputLevel::Quiet ) { + std::cout << " : Fixing parameter " << i << std::endl; + } minuit_->FixParameter(i); } } @@ -131,25 +139,26 @@ // for maximum likelihood fit. Very important command, otherwise all // extracted errors will be too big, and pull distributions will be too narrow! // TODO - The alternative to this is to make FCN = -2log(L) rather than -log(L) - Double_t argL[2]; - argL[0] = 0.5; - fitStatus_.status = minuit_->ExecuteCommand("SET ERR", argL, 1); + + std::array argL { 0.5 }; + fitStatus_.status = minuit_->ExecuteCommand("SET ERR", argL.data(), argL.size()); //argL[0] = 0; - //fitStatus_.status = minuit_->ExecuteCommand("SET STRATEGY", argL, 1); + //fitStatus_.status = minuit_->ExecuteCommand("SET STRATEGY", argL.data(), argL.size()); } LauFitObject* LauMinuit::getFitObject() { - return (minuit_!=0) ? dynamic_cast( minuit_->GetObjectFit() ) : 0; + return (minuit_!=nullptr) ? dynamic_cast( minuit_->GetObjectFit() ) : nullptr; } const LauAbsFitter::FitStatus& LauMinuit::minimise() { - Double_t arglist[2]; - arglist[0] = 1000*nParams_; // maximum iterations - arglist[1] = 0.05; // tolerance -> min EDM = 0.001*tolerance (0.05) - fitStatus_.status = minuit_->ExecuteCommand("MIGRAD", arglist, 2); + std::array arglist { + 1000.0*nParams_, // maximum iterations + 0.05 // tolerance -> min EDM = 0.001*tolerance (0.05) + }; + fitStatus_.status = minuit_->ExecuteCommand("MIGRAD", arglist.data(), arglist.size()); // Dummy variables - need to feed them to the function // used for getting NLL, EDM and error matrix status @@ -158,30 +167,38 @@ if (fitStatus_.status != 0) { - std::cerr << "ERROR in LauMinuit::minimise : Error in minimising loglike." << std::endl; + if ( outputLevel_ > LauOutputLevel::None ) { + std::cerr << "ERROR in LauMinuit::minimise : Error in minimising loglike." << std::endl; + } } else { // Check that the error matrix is ok fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); - std::cout << "INFO in LauMinuit::minimise : Error matrix status after MIGRAD is: " << fitStatus_.status << std::endl; + if ( outputLevel_ > LauOutputLevel::Quiet ) { + std::cout << "INFO in LauMinuit::minimise : Error matrix status after MIGRAD is: " << fitStatus_.status << std::endl; + } // 0= not calculated at all // 1= approximation only, not accurate // 2= full matrix, but forced positive-definite // 3= full accurate covariance matrix // Fit result was OK. Now get the more precise errors. - fitStatus_.status = minuit_->ExecuteCommand("HESSE", arglist, 1); + fitStatus_.status = minuit_->ExecuteCommand("HESSE", arglist.data(), 1); if (fitStatus_.status != 0) { - std::cerr << "ERROR in LauMinuit::minimise : Error in HESSE routine." << std::endl; + if ( outputLevel_ > LauOutputLevel::None ) { + std::cerr << "ERROR in LauMinuit::minimise : Error in HESSE routine." << std::endl; + } } else { // Check that the error matrix is ok fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); - std::cout << "INFO in LauMinuit::minimise : Error matrix status after HESSE is: " << fitStatus_.status << std::endl; + if ( outputLevel_ > LauOutputLevel::Quiet ) { + std::cout << "INFO in LauMinuit::minimise : Error matrix status after HESSE is: " << fitStatus_.status << std::endl; + } // 0= not calculated at all // 1= approximation only, not accurate // 2= full matrix, but forced positive-definite @@ -192,7 +209,7 @@ if (useAsymmFitErrors_ == kTRUE) { LauFitObject* fitObj = this->getFitObject(); fitObj->withinAsymErrorCalc( kTRUE ); - fitStatus_.status = minuit_->ExecuteCommand("MINOS", arglist, 1); + fitStatus_.status = minuit_->ExecuteCommand("MINOS", arglist.data(), 1); fitObj->withinAsymErrorCalc( kFALSE ); if (fitStatus_.status != 0) { std::cerr << "ERROR in LauMinuit::minimise : Error in MINOS routine." << std::endl; @@ -203,13 +220,17 @@ // Print results fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); - std::cout << "INFO in LauMinuit::minimise : Final error matrix status is: " << fitStatus_.status << std::endl; + if ( outputLevel_ > LauOutputLevel::None ) { + std::cout << "INFO in LauMinuit::minimise : Final error matrix status is: " << fitStatus_.status << std::endl; + } // 0= not calculated at all // 1= approximation only, not accurate // 2= full matrix, but forced positive-definite // 3= full accurate covariance matrix - minuit_->PrintResults(3, fitStatus_.NLL); + if ( outputLevel_ > LauOutputLevel::None ) { + minuit_->PrintResults(3, fitStatus_.NLL); + } // Retrieve the covariance matrix from the fitter // For some reason the array returned is as if the matrix is of dimension nParams_ x nParams_ @@ -226,9 +247,8 @@ void LauMinuit::fixSecondStageParameters() { - for (UInt_t i = 0; i < nParams_; ++i) { - Bool_t secondStage = params_[i]->secondStage(); - if (secondStage == kTRUE) { + for (UInt_t i{0}; i < nParams_; ++i) { + if ( params_[i]->secondStage() ) { params_[i]->fixed(kTRUE); minuit_->FixParameter(i); } @@ -239,9 +259,8 @@ void LauMinuit::releaseSecondStageParameters() { - for (UInt_t i = 0; i < nParams_; ++i) { - Bool_t secondStage = params_[i]->secondStage(); - if (secondStage == kTRUE) { + for (UInt_t i{0}; i < nParams_; ++i) { + if ( params_[i]->secondStage() ) { params_[i]->fixed(kFALSE); minuit_->ReleaseParameter(i); } @@ -252,13 +271,13 @@ void LauMinuit::updateParameters() { - for (UInt_t i = 0; i < nParams_; ++i) { + for (UInt_t i{0}; i < nParams_; ++i) { // Get the value and errors from MINUIT - Double_t value = minuit_->GetParameter(i); - Double_t error(0.0); - Double_t negError(0.0); - Double_t posError(0.0); - Double_t globalcc(0.0); + Double_t value { minuit_->GetParameter(i) }; + Double_t error{0.0}; + Double_t negError{0.0}; + Double_t posError{0.0}; + Double_t globalcc{0.0}; minuit_->GetErrors(i, posError, negError, error, globalcc); params_[i]->valueAndErrors(value, error, negError, posError); params_[i]->globalCorrelationCoeff(globalcc); @@ -271,7 +290,7 @@ // Routine that specifies the negative log-likelihood function for the fit. // Used by the MINUIT minimising code. - LauFitObject* theModel = LauFitter::fitter()->getFitObject(); + LauFitObject* theModel = LauFitter::fitter().getFitObject(); // Set the internal parameters for this model using parameters from Minuit (pars): theModel->setParsFromMinuit( par, npar ); diff --git a/src/LauSimFitCoordinator.cc b/src/LauSimFitCoordinator.cc --- a/src/LauSimFitCoordinator.cc +++ b/src/LauSimFitCoordinator.cc @@ -645,36 +645,36 @@ this->checkInitFitParams(); // Initialise the fitter - LauFitter::fitter()->useAsymmFitErrors( this->useAsymmFitErrors() ); - LauFitter::fitter()->twoStageFit( this->twoStageFit() ); - LauFitter::fitter()->initialise( this, params_ ); + LauFitter::fitter().useAsymmFitErrors( this->useAsymmFitErrors() ); + LauFitter::fitter().twoStageFit( this->twoStageFit() ); + LauFitter::fitter().initialise( this, params_ ); - this->startNewFit( LauFitter::fitter()->nParameters(), LauFitter::fitter()->nFreeParameters() ); + this->startNewFit( LauFitter::fitter().nParameters(), LauFitter::fitter().nFreeParameters() ); // Now ready for minimisation step std::cout << "\nINFO in LauSimFitCoordinator::fitExpt : Start minimisation...\n"; - LauAbsFitter::FitStatus fitResult = LauFitter::fitter()->minimise(); + LauAbsFitter::FitStatus fitResult = LauFitter::fitter().minimise(); // If we're doing a two stage fit we can now release (i.e. float) // the 2nd stage parameters and re-fit if (this->twoStageFit()) { if ( fitResult.status != 3 ) { std::cerr << "ERROR in LauSimFitCoordinator:fitExpt : Not running second stage fit since first stage failed." << std::endl; - LauFitter::fitter()->releaseSecondStageParameters(); + LauFitter::fitter().releaseSecondStageParameters(); } else { - LauFitter::fitter()->releaseSecondStageParameters(); - this->startNewFit( LauFitter::fitter()->nParameters(), LauFitter::fitter()->nFreeParameters() ); - fitResult = LauFitter::fitter()->minimise(); + LauFitter::fitter().releaseSecondStageParameters(); + this->startNewFit( LauFitter::fitter().nParameters(), LauFitter::fitter().nFreeParameters() ); + fitResult = LauFitter::fitter().minimise(); } } - const TMatrixD& covMat = LauFitter::fitter()->covarianceMatrix(); + const TMatrixD& covMat = LauFitter::fitter().covarianceMatrix(); this->storeFitStatus( fitResult, covMat ); // Store the final fit results and errors into protected internal vectors that // all sub-classes can use within their own finalFitResults implementation // used below (e.g. putting them into an ntuple in a root file) - LauFitter::fitter()->updateParameters(); + LauFitter::fitter().updateParameters(); } void LauSimFitCoordinator::setParsFromMinuit(Double_t* par, Int_t npar)