diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc index cf7c363..f5cd131 100644 --- a/examples/Test_Dpipi.cc +++ b/examples/Test_Dpipi.cc @@ -1,444 +1,444 @@ #include #include #include #include #include "TFile.h" #include "TH2.h" #include "TRandom.h" #include "TString.h" #include "TSystem.h" #include "TF1.h" #include "TCanvas.h" #include "Lau1DHistPdf.hh" #include "Lau1DCubicSpline.hh" #include "LauBinnedDecayTimeEfficiency.hh" #include "LauBkgndDPModel.hh" #include "LauDaughters.hh" #include "LauDecayTime.hh" #include "LauDecayTimePdf.hh" #include "LauDecayTimePhysicsModel.hh" #include "LauDecayTimeResolution.hh" #include "LauEffModel.hh" #include "LauFlavTag.hh" #include "LauIsobarDynamics.hh" #include "LauMagPhaseCoeffSet.hh" #include "LauRandom.hh" #include "LauRealImagCoeffSet.hh" #include "LauSplineDecayTimeEfficiency.hh" #include "LauTimeDepFitModel.hh" #include "LauVetoes.hh" #include "Test_Dpipi_ProgOpts.hh" int main(const int argc, const char ** argv) { const TestDpipi_ProgramSettings settings{argc,argv}; if ( settings.helpRequested ) { return EXIT_SUCCESS; } if ( ! settings.parsedOK ) { return EXIT_FAILURE; } LauRandom::setSeed(settings.RNGseed); LauDaughters* daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0"); LauDaughters* daughtersB0 = new LauDaughters("B0", "pi+", "pi-", "D0_bar"); // efficiency LauVetoes* vetoes = new LauVetoes(); LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes); LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes); // background types /* std::map bkgndInfo { {"Bkgnd1",LauFlavTag::BkgndType::Combinatorial}, {"Bkgnd2",LauFlavTag::BkgndType::SignalLike} }; */ // setup flavour tagging const Bool_t useAveDeltaCalibVals { settings.useAveDeltaCalibVals }; const Bool_t useEtaPrime { kFALSE }; LauFlavTag* flavTag = new LauFlavTag(useAveDeltaCalibVals,useEtaPrime); //LauFlavTag* flavTag = new LauFlavTag(useAveDeltaCalibVals,useEtaPrime,BkgndTypes); if (settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS) { flavTag->setDecayFlvVarName("decayFlv"); } flavTag->setTrueTagVarName("trueTag"); TFile* etaFile = TFile::Open("ft-eta-hist.root"); TH1* etaHist = dynamic_cast(etaFile->Get("ft_eta_hist")); // Crude check as to whether we're doing perfect vs realistic mis-tag // - in the former case all entries should be in the first bin // If the tagging is perfect then don't interpolate the eta histogram // and also make it perfectly efficient, otherwise do interpolate and // make it 50% efficient const Double_t meanEta { etaHist->GetMean() }; const Bool_t interpolateEtaHist { meanEta > etaHist->GetBinWidth(1) }; Lau1DHistPdf* etaHistPDF = new Lau1DHistPdf( "eta", etaHist, 0.0, 0.5, interpolateEtaHist, kFALSE ); const Double_t tagEffVal { (interpolateEtaHist) ? 0.5 : 1.0 }; std::pair tagEff {tagEffVal, useAveDeltaCalibVals ? 0.0 : tagEffVal}; // use a null calibration for the time being, so p0 = and p1 = 1 std::pair calib0 {meanEta, useAveDeltaCalibVals ? 0.0 : meanEta}; std::pair calib1 {1.0, useAveDeltaCalibVals ? 0.0 : 1.0}; flavTag->addTagger("OSTagger", "tagVal_OS", "mistagVal_OS", etaHistPDF, tagEff, calib0, calib1); flavTag->floatAllCalibPars(); // flavour tagging setup for backgrounds /* std::map BkgndEtas; TFile* bkgndEtaFile = TFile::Open("ft-bkgnd-eta-hists.root"); for(auto& [name, _] : BkgndTypes) { TH1* bkgndEtaHist = dynamic_cast(bkgndEtaFile->Get( name+"_eta_hist" )); Lau1DHistPdf* bkgndEtaHistPDF = new Lau1DHistPdf("eta",bkgndEtaHist,0.0,0.5,kTRUE,kFALSE); BkgndEtas.emplace( std::make_pair(name, bkgndEtaHistPDF) ); } for(auto& [name,hist] : BkgndEtas) {flavTag->setBkgndParams(name, "IFT", hist, tagEff );} */ // signal dynamics LauIsobarDynamics* sigModelB0bar = new LauIsobarDynamics(daughtersB0bar, effModelB0bar); sigModelB0bar->setIntFileName("integ_B0bar.dat"); sigModelB0bar->addResonance("D*+_2", 2, LauAbsResonance::RelBW); sigModelB0bar->addResonance("D*+_0", 2, LauAbsResonance::RelBW); sigModelB0bar->addResonance("rho0(770)", 3, LauAbsResonance::RelBW); sigModelB0bar->addResonance("f_0(980)", 3, LauAbsResonance::RelBW); sigModelB0bar->addResonance("f_2(1270)", 3, LauAbsResonance::RelBW); LauIsobarDynamics* sigModelB0 = new LauIsobarDynamics(daughtersB0, effModelB0); sigModelB0->setIntFileName("integ_B0.dat"); sigModelB0->addResonance("D*-_2", 1, LauAbsResonance::RelBW); sigModelB0->addResonance("D*-_0", 1, LauAbsResonance::RelBW); sigModelB0->addResonance("rho0(770)", 3, LauAbsResonance::RelBW); sigModelB0->addResonance("f_0(980)", 3, LauAbsResonance::RelBW); sigModelB0->addResonance("f_2(1270)", 3, LauAbsResonance::RelBW); // fit model LauTimeDepFitModel* fitModel = new LauTimeDepFitModel(sigModelB0bar,sigModelB0,flavTag); std::vector coeffset; coeffset.push_back( new LauRealImagCoeffSet("D*+_2", 1.00, 0.00, kTRUE, kTRUE) ); coeffset.push_back( new LauRealImagCoeffSet("D*+_0", 0.53*TMath::Cos( 3.00), 0.53*TMath::Sin( 3.00), kFALSE, kFALSE) ); coeffset.push_back( new LauRealImagCoeffSet("rho0(770)", 1.22*TMath::Cos( 2.25), 1.22*TMath::Sin( 2.25), kFALSE, kFALSE) ); coeffset.push_back( new LauRealImagCoeffSet("f_0(980)", 0.19*TMath::Cos(-2.48), 0.19*TMath::Sin(-2.48), kFALSE, kFALSE) ); coeffset.push_back( new LauRealImagCoeffSet("f_2(1270)", 0.75*TMath::Cos( 2.97), 0.75*TMath::Sin( 2.97), kFALSE, kFALSE) ); for (std::vector::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) { fitModel->setAmpCoeffSet(*iter); } // background DP models /* LauBkgndDPModel* Bkgnd1Model = new LauBkgndDPModel(daughtersB0, vetoes); LauBkgndDPModel* Bkgnd1Modelbar = new LauBkgndDPModel(daughtersB0bar, vetoes); LauBkgndDPModel* Bkgnd2Model = new LauBkgndDPModel(daughtersB0, vetoes); LauBkgndDPModel* Bkgnd2Modelbar = new LauBkgndDPModel(daughtersB0bar, vetoes); fitModel->setBkgndDPModels( "Bkgnd1", Bkgnd1Model, Bkgnd1Modelbar ); fitModel->setBkgndDPModels( "Bkgnd2", Bkgnd2Model, Bkgnd2Modelbar ); */ // decay type and mixing parameter const Bool_t fixPhiMix{ settings.fixPhiMix || settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS }; const Bool_t useSinCos{ settings.useSinCos }; fitModel->setCPEigenvalue( settings.dType ); fitModel->setPhiMix( 2.0*LauConstants::beta, fixPhiMix, useSinCos ); // production asymmetries fitModel->setAsymmetries( 0.0, kTRUE ); /* for(auto& [name, _] : BkgndTypes) { fitModel->setBkgndAsymmetries( name, 0.0, kTRUE ); } */ // decay time PDFs const Double_t minDt(0.0); const Double_t maxDt(15.0); const Double_t minDtErr(0.0); const Double_t maxDtErr(0.15); LauParameter * tau = new LauParameter("dt_tau", 1.519, 0.5, 5.0, settings.fixLifetime); LauParameter * freq = new LauParameter("dt_deltaM", 0.5064, 0.0, 1.0, settings.fixDeltaM); std::vector dtPhysPars { tau, freq }; auto dtPhysModel = std::make_unique( LauDecayTime::FuncType::ExpTrig, dtPhysPars ); const std::vector scale { settings.perEventTimeErr && true, settings.perEventTimeErr && true, }; const std::size_t nGauss{scale.size()}; LauParameter * mean0 = new LauParameter("dt_mean_0", scale[0] ? -2.01290e-03 : -2.25084e-03, -0.01, 0.01, kTRUE ); LauParameter * mean1 = new LauParameter("dt_mean_1", scale[1] ? -2.01290e-03 : -5.04275e-03, -0.01, 0.01, kTRUE ); LauParameter * sigma0 = new LauParameter("dt_sigma_0", scale[0] ? 9.95145e-01 : 3.03923e-02, 0.0, 2.0, kTRUE ); LauParameter * sigma1 = new LauParameter("dt_sigma_1", scale[1] ? 1.81715e+00 : 6.22376e-02, 0.0, 2.5, kTRUE ); LauParameter * frac1 = new LauParameter("dt_frac_1", scale[0] && scale[1] ? 1.-9.35940e-01 : 1.-7.69603e-01, 0.0, 1.0, kTRUE); std::vector dtResoPars { mean0, mean1, sigma0, sigma1, frac1 }; auto dtResoModel = std::make_unique( nGauss, dtResoPars, scale ); // Decay time error histogram // (always set this so that it gets generated properly, whether we're using it in the PDF or not) TFile* dtErrFile = TFile::Open("dte-hist.root"); TH1* dtErrHist = dynamic_cast(dtErrFile->Get("dte_hist")); LauDecayTimePdf * dtPdf{nullptr}; if ( settings.timeResolution ) { if ( settings.perEventTimeErr ) { dtPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, "decayTimeErr", minDtErr, maxDtErr, std::move(dtPhysModel), std::move(dtResoModel), dtErrHist ); } else { dtPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(dtPhysModel), std::move(dtResoModel) ); } } else { dtPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(dtPhysModel) ); } // Decay time acceptance histogram TFile* dtaFile = TFile::Open("dta-hist.root"); TH1* dtaHist = dynamic_cast(dtaFile->Get("dta_hist")); // Create the spline knot positions and // starting Y values, to be fit to dtaHist const std::vector dtvals { 0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 7.0, 10.0, 15.0}; const std::vector effvals {0.000, 0.010, 0.022, 0.035, 0.042, 0.050, 0.051, 0.052, 0.055}; const std::vector corvals { 1.10, 1.08, 1.06, 1.04, 1.02, 1.00, 0.98, 0.96, 0.94}; switch(settings.timeEffModel) { case LauDecayTime::EfficiencyMethod::Spline: { fitModel->setASqMaxValue(0.06); - auto dtEffSpline = std::make_unique(dtvals,effvals,Lau1DCubicSpline::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); + auto dtEffSpline = std::make_unique( dtvals, effvals, Lau1DCubicSpline::SplineType::AkimaSpline ); auto dtaModel = std::make_unique>( "dteff_QFS", std::move(dtEffSpline) ); dtaModel->fitToTH1(dtaHist); // Set which knots to float and which to fix (at least 1 knot must be fixed, not the first one) // Knots should only be floating if requested AND the B lifetime is fixed! if ( settings.fixSplineKnots or not settings.fixLifetime ) { dtaModel->fixKnots(); } else { dtaModel->floatKnots(); dtaModel->fixKnot( 0, true ); dtaModel->fixKnot( 3, true ); } if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) { // For the QFS mode we just use the cubic model as it is dtPdf->setSplineEfficiency( std::move(dtaModel) ); } else { // For the CP modes we modify it using a corrective spline - auto dtCorrSpline = std::make_unique(dtvals,corvals,Lau1DCubicSpline::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); + auto dtCorrSpline = std::make_unique( dtvals, corvals, Lau1DCubicSpline::SplineType::AkimaSpline ); auto dtaCPModel = std::make_unique>( "dteff_CP", *dtaModel, std::move( dtCorrSpline ) ); dtPdf->setSplineEfficiency( std::move(dtaCPModel) ); } break; } case LauDecayTime::EfficiencyMethod::Binned: { fitModel->setASqMaxValue(0.06); auto dtaBinnedModel = std::make_unique( *dtaHist ); dtPdf->setBinnedEfficiency( std::move(dtaBinnedModel) ); break; } case LauDecayTime::EfficiencyMethod::Uniform: { fitModel->setASqMaxValue(4.45); break; } } fitModel->setSignalDtPdf( dtPdf ); //Background dt part /* TFile* background_dt = TFile::Open("Lifetimes_PV_WL.root"); TH1* bkgnd1DtHist = dynamic_cast( background_dt->Get("Bkgnd1") ); TH1* bkgnd2DtHist = dynamic_cast( background_dt->Get("Bkgnd2") ); TH1* bkgnd1DtErrHist = dynamic_cast( background_dt->Get("Bkgnd1_Err") ); TH1* bkgnd2DtErrHist = dynamic_cast( background_dt->Get("Bkgnd2_Err") ); LauDecayTimePdf* bkgnd1DtPdf{nullptr}; LauDecayTimePdf* bkgnd2DtPdf{nullptr}; if ( settings.timeResolution and settings.perEventTimeErr ) { bkgnd1DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd1DtHist, bkgnd1DtErrHist ); bkgnd2DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd2DtHist, bkgnd2DtErrHist ); } else { bkgnd1DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd1DtHist ); bkgnd2DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd2DtHist ); } fitModel->setBkgndDtPdf("Bkgnd1",bkgnd1DtPdf); fitModel->setBkgndDtPdf("Bkgnd2",bkgnd2DtPdf); */ // set the signal yield Double_t nSigEvents{0}; switch (settings.dType) { case LauTimeDepFitModel::CPEigenvalue::CPEven : nSigEvents = 15000; break; case LauTimeDepFitModel::CPEigenvalue::CPOdd : nSigEvents = 5000; break; case LauTimeDepFitModel::CPEigenvalue::QFS : nSigEvents = 50000; break; } // set the background yields /* const Double_t nBkgnd1(100), nBkgnd2(200); LauParameter* Bkgnd1Yield = new LauParameter("Bkgnd1",nBkgnd1,-1.0*nBkgnd1,2.0*nBkgnd1,kFALSE); LauParameter* Bkgnd2Yield = new LauParameter("Bkgnd2",nBkgnd2,-1.0*nBkgnd2,2.0*nBkgnd2,kFALSE); fitModel->setNBkgndEvents(Bkgnd1Yield); fitModel->setNBkgndEvents(Bkgnd2Yield); */ std::cout<<"nSigEvents = "<setNSigEvents(nSigPar); // set the number of experiments if (settings.command == Command::Generate) { fitModel->setNExpts(settings.nExptGen, settings.firstExptGen); } else { fitModel->setNExpts(settings.nExptFit, settings.firstExptFit); } fitModel->useAsymmFitErrors(kFALSE); fitModel->useRandomInitFitPars(kFALSE); fitModel->writeLatexTable(kFALSE); const Bool_t haveBkgnds = ( fitModel->nBkgndClasses() > 0 ); fitModel->doPoissonSmearing(haveBkgnds); fitModel->doEMLFit(haveBkgnds); TString dTypeStr; switch (settings.dType) { case LauTimeDepFitModel::CPEigenvalue::CPEven : dTypeStr = "CPEven"; break; case LauTimeDepFitModel::CPEigenvalue::CPOdd : dTypeStr = "CPOdd"; break; case LauTimeDepFitModel::CPEigenvalue::QFS : dTypeStr = "QFS"; break; } TString dataFile(""); TString treeName("fitTree"); TString rootFileName(""); TString tableFileName(""); TString fitToyFileName(""); TString splotFileName(""); dataFile = "TEST-Dpipi"; dataFile += "_"+dTypeStr; switch(settings.timeEffModel) { case LauDecayTime::EfficiencyMethod::Spline: dataFile += "_Spline"; break; case LauDecayTime::EfficiencyMethod::Binned: dataFile += "_Binned"; break; case LauDecayTime::EfficiencyMethod::Uniform: dataFile += "_Uniform"; break; } if (settings.timeResolution) { if (settings.perEventTimeErr) { dataFile += "_DTRperevt"; } else { dataFile += "_DTRavg"; } } else { dataFile += "_DTRoff"; } dataFile += "_expts"; dataFile += settings.firstExptGen; dataFile += "-"; dataFile += settings.firstExptGen+settings.nExptGen-1; dataFile += ".root"; if (settings.command == Command::Generate) { rootFileName = "dummy.root"; tableFileName = "genResults"; } else { rootFileName = "fit"; rootFileName += settings.iFit; rootFileName += "_Results_"; rootFileName += dTypeStr; rootFileName += "_expts"; rootFileName += settings.firstExptFit; rootFileName += "-"; rootFileName += settings.firstExptFit+settings.nExptFit-1; rootFileName += ".root"; tableFileName = "fit"; tableFileName += settings.iFit; tableFileName += "_Results_"; tableFileName += dTypeStr; tableFileName += "_expts"; tableFileName += settings.firstExptFit; tableFileName += "-"; tableFileName += settings.firstExptFit+settings.nExptFit-1; fitToyFileName = "fit"; fitToyFileName += settings.iFit; fitToyFileName += "_ToyMC_"; fitToyFileName += dTypeStr; fitToyFileName += "_expts"; fitToyFileName += settings.firstExptFit; fitToyFileName += "-"; fitToyFileName += settings.firstExptFit+settings.nExptFit-1; fitToyFileName += ".root"; splotFileName = "fit"; splotFileName += settings.iFit; splotFileName += "_sPlot_"; splotFileName += dTypeStr; splotFileName += "_expts"; splotFileName += settings.firstExptFit; splotFileName += "-"; splotFileName += settings.firstExptFit+settings.nExptFit-1; splotFileName += ".root"; } // Generate toy from the fitted parameters //fitModel->compareFitData(1, fitToyFileName); // Write out per-event likelihoods and sWeights //fitModel->writeSPlotData(splotFileName, "splot", kFALSE); // Execute the generation/fit switch (settings.command) { case Command::Generate : fitModel->run( "gen", dataFile, treeName, rootFileName, tableFileName ); break; case Command::Fit : fitModel->run( "fit", dataFile, treeName, rootFileName, tableFileName ); break; case Command::SimFit : fitModel->runTask( dataFile, treeName, rootFileName, tableFileName, "localhost", settings.port ); break; } return EXIT_SUCCESS; } diff --git a/inc/Lau1DCubicSpline.hh b/inc/Lau1DCubicSpline.hh index 0cc2f6e..83f1bd1 100644 --- a/inc/Lau1DCubicSpline.hh +++ b/inc/Lau1DCubicSpline.hh @@ -1,261 +1,261 @@ /* Copyright 2015 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 Lau1DCubicSpline.hh \brief File containing declaration of Lau1DCubicSpline class. */ /*! \class Lau1DCubicSpline \brief Class for defining a 1D cubic spline based on a set of knots. Class for defining a 1D cubic spline based on a set of knots. Interpolation between the knots is performed either by one of two types of cubic spline (standard or Akima) or by linear interpolation. The splines are defined by a piecewise cubic function which, between knots i and i+1, has the form f_i(x) = (1 - t)*y_i + t*y_i+1 + t*(1 - t)(a*(1 - t) + b*t) where t = (x - x_i)/(x_i+1 - x_i), a = k_i *(x_i+1 - x_i) - (y_i+1 - y_i), b = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i) and k_i is (by construction) the first derivative at knot i. f(x) and f'(x) are continuous at the internal knots by construction. For the standard splines, f''(x) is required to be continuous at all internal knots placing n-2 constraints on the n parameters, k_i. The final two constraints are set by the boundary conditions. At each boundary, the function may be: (i) Clamped : f'(x) = C at the last knot (ii) Natural : f''(x) = 0 at the last knot (iii) Not a knot : f'''(x) continuous at the second last knot The algorithms used in these splines can be found on: http://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm For the Akima splines, the values of k_i are determined from the slopes of the four nearest segments (a_i-1, a_i, a_i+1 and a_i+2) as k_i = ( | a_i+2 - a_i+1 | a_i + | a_i - a_i-1 | a_i+1 ) / ( | a_i+2 - a_i+1 | + | a_i - a_i-1 | ) and as k_i = ( a_i + a_i+1 ) / 2 in the special case a_i-1 == a_i != a_i+1 == a_i+2. Boundary conditions are specified by the relations a_2 - a_1 = a_1 - a_0 = a_0 - a_-1 and a_n-1 - a_n-2 = a_n - a_n-1 = a_n+1 - a_n The algorithms used in these splines can be found in: J.ACM vol. 17 no. 4 pp 589-602 */ #ifndef LAU_1DCUBICSPLINE #define LAU_1DCUBICSPLINE #include #include #include "Rtypes.h" #include "TF1.h" #include "TFitResultPtr.h" #include "LauPrint.hh" class TH1; class LauAbsRValue; class LauParameter; class Lau1DCubicSpline final { public: //! Define the allowed interpolation types - enum LauSplineType { + enum class SplineType { StandardSpline, /*!< standard cubic splines with f''(x) continuous at all internal knots */ AkimaSpline, /*!< Akima cubic splines with f'(x) at each knot defined locally by the positions of only five knots */ LinearInterpolation /*! Linear interpolation between each pair of knots */ }; //! Define the allowed boundary condition types /*! These are only supported by standard splines */ - enum LauSplineBoundaryType { + enum class BoundaryType { Clamped, /*!< clamped boundary - f'(x) = C */ Natural, /*!< natural boundary - f''(x) = 0 */ NotAKnot /*!< 'not a knot' boundary - f'''(x) continuous at second last knot */ }; //! Constructor /*! \param [in] xs the x-values of the knots \param [in] ys the y-values of the knots \param [in] type the type of spline (Standard, Akima, Linear) \param [in] leftBound the left-hand boundary condition \param [in] rightBound the right-hand boundary condition \param [in] dydx0 the gradient at the left-hand boundary of a clamped spline \param [in] dydxn the gradient at the right-hand boundary of a clamped spline */ Lau1DCubicSpline(const std::vector& xs, const std::vector& ys, - LauSplineType type = Lau1DCubicSpline::StandardSpline, - LauSplineBoundaryType leftBound = Lau1DCubicSpline::NotAKnot, - LauSplineBoundaryType rightBound = Lau1DCubicSpline::NotAKnot, - Double_t dydx0 = 0.0, Double_t dydxn = 0.0); + const SplineType type = SplineType::StandardSpline, + const BoundaryType leftBound = BoundaryType::NotAKnot, + const BoundaryType rightBound = BoundaryType::NotAKnot, + const Double_t dydx0 = 0.0, const Double_t dydxn = 0.0); //! Evaluate the function at given point /*! \param [in] x the x-coordinate \return the value of the spline at x */ - Double_t evaluate(Double_t x) const; + Double_t evaluate(const Double_t x) const; //! Update the y-values of the knots /*! \param [in] ys the y-values of the knots */ void updateYValues(const std::vector& ys); //! Update the y-values of the knots /*! \param [in] ys the y-values of the knots */ void updateYValues(const std::vector& ys); //! Update the y-values of the knots /*! \param [in] ys the y-values of the knots */ void updateYValues(const std::vector& ys); //! Update the type of interpolation to perform /*! \param [in] type the type of interpolation */ - void updateType(LauSplineType type); + void updateType(const SplineType type); //! Update the boundary conditions for the spline /*! \param [in] leftBound the left-hand boundary condition \param [in] rightBound the right-hand boundary condition \param [in] dydx0 the gradient at the left-hand boundary of a clamped spline \param [in] dydxn the gradient at the right-hand boundary of a clamped spline */ - void updateBoundaryConditions(LauSplineBoundaryType leftBound, - LauSplineBoundaryType rightBound, - Double_t dydx0 = 0.0, - Double_t dydxn = 0.0); + void updateBoundaryConditions(const BoundaryType leftBound, + const BoundaryType rightBound, + const Double_t dydx0 = 0.0, + const Double_t dydxn = 0.0); //! Get the number of knots std::size_t getnKnots() const {return nKnots_;} //! Get y values const std::vector& getYValues() const {return y_;} //! Get x values const std::vector& getXValues() const {return x_;} //! Get the coefficients of a given spline segment in the form a + bx + cx^2 + dx^3 /*! \param [in] segIndex refers to the index of the knot at the left end of the segment (segIndex = 0 gets the coefficients of the the segment between x_0 and x_1) \param [in] normalise if true, the coefficients returned will be normalised by the integral of the complete spline (defaults to false) \return coefficients {a, b, c, d} */ std::array getCoefficients(const std::size_t segIndex, const bool normalise = false) const; //! Get the integral over all the spline segments Double_t integral() const; //! Make a TF1 showing the spline with its current knot values /*! \param [in] normalise whether or not you want the spline normalised \return 1D function object */ TF1* makeTF1(const bool normalise = false) const; //! Fit the the normalisation of the spline to a TH1 /*! \param [in] hist the histogram to be fitted \param [in] printLevel the level of printout desired from fit \return a TF1 fit to the histogram */ - TF1* normaliseToTH1(TH1* hist, LauOutputLevel printLevel = LauOutputLevel::Standard) const; + TF1* normaliseToTH1(TH1* hist, const LauOutputLevel printLevel = LauOutputLevel::Standard) const; //! Fit the spline to a TH1 /*! \param [in] hist the histogram to be fitted \param [in] printLevel the level of printout desired from fit \param [in] doWL If true do a weighted likelihood fit, else chi2 \return a shared-ownership smart pointer to the fit results */ - TFitResultPtr fitToTH1(TH1* hist, LauOutputLevel printLevel = LauOutputLevel::Standard, const bool doWL = false); + TFitResultPtr fitToTH1(TH1* hist, const LauOutputLevel printLevel = LauOutputLevel::Standard, const bool doWL = false); private: //! Initialise the class void init(); //! Calculate the first derivative at each knot void calcDerivatives(); //! Calculate the first derivatives according to the standard method void calcDerivativesStandard(); //! Calculate the first derivatives according to the Akima method void calcDerivativesAkima(); //! The number of knots in the spline const std::size_t nKnots_; //! The x-value at each knot - std::vector x_; + const std::vector x_; //! The y-value at each knot std::vector y_; //! The first derivative at each knot std::vector dydx_; //! The 'a' coefficients used to determine the derivatives std::vector a_; //! The 'b' coefficients used to determine the derivatives std::vector b_; //! The 'c' coefficients used to determine the derivatives std::vector c_; //! The 'd' coefficients used to determine the derivatives std::vector d_; //! The type of interpolation to be performed - LauSplineType type_; + SplineType type_; //! The left-hand boundary condition on the spline - LauSplineBoundaryType leftBound_; + BoundaryType leftBound_; //! The right-hand boundary condition on the spline - LauSplineBoundaryType rightBound_; + BoundaryType rightBound_; //! The gradient at the left boundary for a clamped spline Double_t dydx0_; //! The gradient at the right boundary for a clamped spline Double_t dydxn_; ClassDef(Lau1DCubicSpline, 0) // Class for defining a 1D cubic spline }; #endif diff --git a/inc/LauAbsModIndPartWave.hh b/inc/LauAbsModIndPartWave.hh index ead918f..113df1e 100644 --- a/inc/LauAbsModIndPartWave.hh +++ b/inc/LauAbsModIndPartWave.hh @@ -1,249 +1,247 @@ /* Copyright 2014 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 LauAbsModIndPartWave.hh \brief File containing declaration of LauAbsModIndPartWave class. */ /*! \class LauAbsModIndPartWave \brief Abstract base class for defining a model independent partial wave component Abstract base class for defining a model independent partial wave component. This model uses splines to produce a partial wave from two sets of real numbers that represent the amplitude at a series of points in the phase space. These real numbers at each point can be floated in the fit. Classes inheriting from this define whether these real numbers are e.g. the magnitude and phase or e.g. the real and imaginary part of the amplitude. */ #ifndef LAU_ABSMODINDPARTWAVE #define LAU_ABSMODINDPARTWAVE +#include #include #include "LauComplex.hh" #include "LauAbsResonance.hh" #include "Lau1DCubicSpline.hh" class LauAbsModIndPartWave : public LauAbsResonance { public: //! Constructor /*! \param [in] resInfo the object containing information on the resonance name, mass, width, spin, charge, etc. \param [in] resPairAmpInt the number of the daughter not produced by the resonance \param [in] daughters the daughter particles */ LauAbsModIndPartWave(LauResonanceInfo* resInfo, Int_t resPairAmpInt, const LauDaughters* daughters); - //! Destructor - virtual ~LauAbsModIndPartWave(); - //! Initialise the model virtual void initialise(); //! Define the knot positions /*! If absent from the set provided, knots are added automatically at the upper and lower kinematic limits \param [in] masses the mass values at which the knots are placed */ void defineKnots(const std::set& masses); //! Return the number of knots that have been defined (including those at the upper and lower kinematic limits) /*! \return the number of knots */ UInt_t nKnots() const { return nKnots_; } //! Set the values of the two real parameters that define the amplitude at a given knot /*! \param [in] knot the knot to be updated \param [in] ampVal1 the value of first real parameter representing the amplitude at the knot \param [in] ampVal2 the value of second real parameter representing the amplitude at the knot \param [in] fixAmpVal1 whether the first real parameter should be fixed \param [in] fixAmpVal2 whether the second real parameter should be fixed */ virtual void setKnotAmp(const UInt_t knot, const Double_t ampVal1, const Double_t ampVal2, const Bool_t fixAmpVal1, const Bool_t fixAmpVal2) = 0; //! Set whether the parameters should be floated only in the second-stage of a two-stage fit /*! By default, the parameters describing the amplitude at each knot will float from the outset of a fit. If, however, a good estimate of these is already known, it can be more efficient to initially fix them and then to float them only in a second stage (as is done for other resonance lineshape parameters). This function allows the toggling of this behaviour. \param secondStage whether the parameters should float only in the second stage */ void floatKnotsSecondStage(const Bool_t secondStage); //! Retrieve the value of the second stage flag Bool_t floatKnotsSecondStage() const {return secondStage_;}; //! Retrieve the resonance parameters, e.g. so that they can be loaded into a fit /*! \return floating parameters of the resonance */ virtual const std::vector& getFloatingParameters(); protected: //! Complex resonant amplitude /*! \param [in] mass appropriate invariant mass for the resonance \param [in] spinTerm spin term */ virtual LauComplex resAmp(Double_t mass, Double_t spinTerm); //! Evaluate the amplitude at the given point from the splines /*! \param [in] mass appropriate invariant mass for the resonance */ virtual void evaluateAmplitude(const Double_t mass) = 0; //! Method to check that the supplied knot positions are valid /*! \param [in] masses the mass values at which the knots are placed */ virtual std::set checkKnots(const std::set& masses); //! Method to create the parameter objects for the given knot /*! \param [in] iKnot the index of the knot */ virtual void createAmpParameters(const UInt_t iKnot) = 0; //! Method to set the type of interpolation used for the splines /*! \param [in] type1 the type of interpolation for the first spline \param [in] type2 the type of interpolation for the second spline */ - void setSplineType(Lau1DCubicSpline::LauSplineType type1, Lau1DCubicSpline::LauSplineType type2); + void setSplineType(Lau1DCubicSpline::SplineType type1, Lau1DCubicSpline::SplineType type2); //! Method to set the boundary conditions of the splines /*! \param [in] leftBound1 the type of boundary condition for the left edge of the first spline \param [in] rightBound1 the type of boundary condition for the right edge of the first spline \param [in] leftBound2 the type of boundary condition for the left edge of the second spline \param [in] rightBound2 the type of boundary condition for the right edge of the second spline \param [in] leftGrad1 the gradient at the left edge of the first spline if clamped \param [in] rightGrad1 the gradient at the right edge of the first spline if clamped \param [in] leftGrad2 the gradient at the left edge of the second spline if clamped \param [in] rightGrad2 the gradient at the right edge of the second spline if clamped */ - void setSplineBoundaryConditions(Lau1DCubicSpline::LauSplineBoundaryType leftBound1, - Lau1DCubicSpline::LauSplineBoundaryType rightBound1, - Lau1DCubicSpline::LauSplineBoundaryType leftBound2, - Lau1DCubicSpline::LauSplineBoundaryType rightBound2, + void setSplineBoundaryConditions(Lau1DCubicSpline::BoundaryType leftBound1, + Lau1DCubicSpline::BoundaryType rightBound1, + Lau1DCubicSpline::BoundaryType leftBound2, + Lau1DCubicSpline::BoundaryType rightBound2, Double_t leftGrad1 = 0.0, Double_t rightGrad1 = 0.0, Double_t leftGrad2 = 0.0, Double_t rightGrad2 = 0.0); //! Helper function to set the current amplitude value /*! \param [in] realPart the real part of the amplitude \param [in] imagPart the imaginary part of the amplitude */ void setAmp(const Double_t realPart, const Double_t imagPart) { amp_.setRealImagPart(realPart,imagPart); } //! Helper function to access the masses const std::vector& getMasses() {return masses_;} //! Helper function to access the 1st parameter set std::vector& getAmp1Vals() {return amp1Vals_;} //! Helper function to access the 2nd parameter set std::vector& getAmp2Vals() {return amp2Vals_;} //! Helper function to access the 1st parameter set std::vector& getAmp1Pars() {return amp1Pars_;} //! Helper function to access the 2nd parameter set std::vector& getAmp2Pars() {return amp2Pars_;} //! Helper function to access the 1st spline - const Lau1DCubicSpline* getSpline1() const {return spline1_;} + const Lau1DCubicSpline& getSpline1() const {return *spline1_;} //! Helper function to access the 1st spline - const Lau1DCubicSpline* getSpline2() const {return spline2_;} + const Lau1DCubicSpline& getSpline2() const {return *spline2_;} private: //! The number of knots - UInt_t nKnots_; + UInt_t nKnots_{0}; //! The masses at which knots are defined in the magnitude and phase splines std::vector masses_; //! The values of the first real parameter at each knot std::vector amp1Vals_; //! The values of the second real parameter at each knot std::vector amp2Vals_; //! The parameters for the first real value at the knots std::vector amp1Pars_; //! The parameters for the second real value at the knots std::vector amp2Pars_; //! The spline used to interpolate the values of the first real parameter - Lau1DCubicSpline* spline1_; + std::unique_ptr spline1_; //! The spline used to interpolate the values of the second real parameter - Lau1DCubicSpline* spline2_; + std::unique_ptr spline2_; //! The type of interpolation used for the first spline - Lau1DCubicSpline::LauSplineType type1_; + Lau1DCubicSpline::SplineType type1_{Lau1DCubicSpline::SplineType::StandardSpline}; //! The type of interpolation used for the second spline - Lau1DCubicSpline::LauSplineType type2_; + Lau1DCubicSpline::SplineType type2_{Lau1DCubicSpline::SplineType::StandardSpline}; //! The lower boundary condition type for the first spline - Lau1DCubicSpline::LauSplineBoundaryType leftBound1_; + Lau1DCubicSpline::BoundaryType leftBound1_{Lau1DCubicSpline::BoundaryType::NotAKnot}; //! The upper boundary condition type for the first spline - Lau1DCubicSpline::LauSplineBoundaryType rightBound1_; + Lau1DCubicSpline::BoundaryType rightBound1_{Lau1DCubicSpline::BoundaryType::NotAKnot}; //! The lower boundary condition type for the second spline - Lau1DCubicSpline::LauSplineBoundaryType leftBound2_; + Lau1DCubicSpline::BoundaryType leftBound2_{Lau1DCubicSpline::BoundaryType::NotAKnot}; //! The upper boundary condition type for the second spline - Lau1DCubicSpline::LauSplineBoundaryType rightBound2_; + Lau1DCubicSpline::BoundaryType rightBound2_{Lau1DCubicSpline::BoundaryType::NotAKnot}; //! The gradient at the left boundary for the first spline if clamped - Double_t leftGrad1_; + Double_t leftGrad1_{0.0}; //! The gradient at the right boundary for the first spline if clamped - Double_t rightGrad1_; + Double_t rightGrad1_{0.0}; //! The gradient at the left boundary for the second spline if clamped - Double_t leftGrad2_; + Double_t leftGrad2_{0.0}; //! The gradient at the right boundary for the second spline if clamped - Double_t rightGrad2_; + Double_t rightGrad2_{0.0}; //! Flag to determine if the parameters should be floated only in the second stage of the fit - Bool_t secondStage_; + Bool_t secondStage_{kFALSE}; //! The current value of the amplitude LauComplex amp_; ClassDef(LauAbsModIndPartWave,0) // model independent partial wave }; #endif diff --git a/inc/LauEFKLLMRes.hh b/inc/LauEFKLLMRes.hh index a05868b..888822c 100644 --- a/inc/LauEFKLLMRes.hh +++ b/inc/LauEFKLLMRes.hh @@ -1,147 +1,146 @@ /* Copyright 2015 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 LauEFKLLMRes.hh \brief File containing declaration of LauEFKLLMRes class. */ /*! \class LauEFKLLMRes \brief Class for defining the EFKLLM K-pi S-wave model Class for defining the EFKLLM form-factor model for the K-pi S-wave. The model consists of a tabulated form-factor, which is interpolated using cubic splines (one for magnitude values and one for phase values), multiplied by a mass-dependence (e.g. constant, 1/m^2, etc.). The massFactor resonance parameter is the power of the mass dependence - defaults to zero, i.e. constant. For more details see B. El-Bennich et al. Phys. Rev. D 79, 094005 (2009), arXiv:0902.3645 [hep-ph]. (The acronym EFKLLM is constructed from the surnames of the authors of the above paper.) */ #ifndef LAU_EFKLLM_RES #define LAU_EFKLLM_RES +#include + #include "TString.h" #include "LauComplex.hh" #include "LauAbsResonance.hh" class Lau1DCubicSpline; class LauEFKLLMRes : public LauAbsResonance { public: //! Constructor /*! \param [in] resInfo the object containing information on the resonance name, mass, width, spin, charge, etc. \param [in] resPairAmpInt the number of the daughter not produced by the resonance \param [in] daughters the daughter particles */ LauEFKLLMRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters); - //! Destructor - virtual ~LauEFKLLMRes(); - //! Read the form factor information from text file /*! Creates the splines from the tabulated form factor data. These are shared between all instances of this class. \param [in] inputFile the name of the file to be read */ static void setupFormFactor(const TString& inputFile); //! Initialise the model virtual void initialise(); //! Get the resonance model type /*! \return the resonance model type */ virtual LauAbsResonance::LauResonanceModel getResonanceModel() const {return LauAbsResonance::EFKLLM;} //! Set value of a resonance parameter /*! \param [in] name the name of the parameter to be changed \param [in] value the new parameter value */ virtual void setResonanceParameter(const TString& name, const Double_t value); //! Allow the various parameters to float in the fit /*! \param [in] name the name of the parameter to be floated */ virtual void floatResonanceParameter(const TString& name); //! Access the given resonance parameter /*! \param [in] name the name of the parameter \return the corresponding parameter */ virtual LauParameter* getResonanceParameter(const TString& name); //! Retrieve the resonance parameters, e.g. so that they can be loaded into a fit /*! \return floating parameters of the resonance */ virtual const std::vector& getFloatingParameters(); protected: //! Set the power of the mass dependence /*! \param [in] massFactor the new power of the mass dependence */ void setMassFactor(const Double_t massFactor); //! Get the power of the mass dependence /*! \return the power of the mass dependence */ - Double_t getMassFactor() const {return (massFactor_!=0) ? massFactor_->unblindValue() : 0.0;} + Double_t getMassFactor() const {return (massFactor_!=nullptr) ? massFactor_->unblindValue() : 0.0;} //! See if the mass factor parameter is fixed or floating /*! \return kTRUE if the mass factor parameter is fixed, kFALSE otherwise */ - Bool_t fixMassFactor() const {return (massFactor_!=0) ? massFactor_->fixed() : kTRUE;} + Bool_t fixMassFactor() const {return (massFactor_!=nullptr) ? massFactor_->fixed() : kTRUE;} //! Complex resonant amplitude /*! \param [in] mass appropriate invariant mass for the resonance \param [in] spinTerm spin term */ virtual LauComplex resAmp(Double_t mass, Double_t spinTerm); private: //! Spline describing the magnitude variation of the form-factor - static Lau1DCubicSpline* magSpline_; + static std::unique_ptr magSpline_; //! Spline describing the phase variation of the form-factor - static Lau1DCubicSpline* phaseSpline_; + static std::unique_ptr phaseSpline_; //! The power of the mass dependence - LauParameter* massFactor_; + LauParameter* massFactor_{nullptr}; ClassDef(LauEFKLLMRes,0) // EFKLLM resonance model }; #endif diff --git a/inc/LauModIndPartWaveMagPhase.hh b/inc/LauModIndPartWaveMagPhase.hh index 2149bd2..237b8b1 100644 --- a/inc/LauModIndPartWaveMagPhase.hh +++ b/inc/LauModIndPartWaveMagPhase.hh @@ -1,125 +1,122 @@ /* Copyright 2015 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 LauModIndPartWaveMagPhase.hh \brief File containing declaration of LauModIndPartWaveMagPhase class. */ /*! \class LauModIndPartWaveMagPhase \brief Class for defining a model independent partial wave component where the amplitudes are parameterised in terms of magnitude and phase This model uses splines to produce a partial wave from the magnitude and phase values of the amplitude at a series of points in the phase space. The magnitude and phase at each point can be floated in the fit. */ #ifndef LAU_MODINDPARTWAVE_MAGPHASE #define LAU_MODINDPARTWAVE_MAGPHASE #include "LauComplex.hh" #include "LauAbsModIndPartWave.hh" class TSpline3; class LauModIndPartWaveMagPhase : public LauAbsModIndPartWave { public: //! Constructor /*! \param [in] resInfo the object containing information on the resonance name, mass, width, spin, charge, etc. \param [in] resPairAmpInt the number of the daughter not produced by the resonance \param [in] daughters the daughter particles */ LauModIndPartWaveMagPhase(LauResonanceInfo* resInfo, Int_t resPairAmpInt, const LauDaughters* daughters); - //! Destructor - virtual ~LauModIndPartWaveMagPhase(); - //! Set the values of the two real parameters that define the amplitude at a given knot /*! \param [in] knot the knot to be updated \param [in] magVal the value of the magnitude at the knot \param [in] phaseVal the value of the phase at the knot \param [in] fixMagnitude whether the magnitude should be fixed \param [in] fixPhase whether the phase should be fixed */ virtual void setKnotAmp(const UInt_t knot, const Double_t magVal, const Double_t phaseVal, const Bool_t fixMagnitude, const Bool_t fixPhase); //! Method to set the type of interpolation used for the splines /*! \param [in] magType the type of interpolation for the magnitude spline \param [in] phaseType the type of interpolation for the phase spline */ - void setType(Lau1DCubicSpline::LauSplineType magType, Lau1DCubicSpline::LauSplineType phaseType) { + void setType(Lau1DCubicSpline::SplineType magType, Lau1DCubicSpline::SplineType phaseType) { this->setSplineType(magType,phaseType); } //! Method to set the boundary conditions of the splines /*! \param [in] magLeftBound the type of boundary condition for the left edge of the magnitude spline \param [in] magRightBound the type of boundary condition for the right edge of the magnitude spline \param [in] phaseLeftBound the type of boundary condition for the left edge of the phase spline \param [in] phaseRightBound the type of boundary condition for the right edge of the phase spline \param [in] magLeftGrad the gradient at the left edge of the magnitude spline if clamped \param [in] magRightGrad the gradient at the right edge of the magnitude spline if clamped \param [in] phaseLeftGrad the gradient at the left edge of the phase spline if clamped \param [in] phaseRightGrad the gradient at the right edge of the phase spline if clamped */ - void setBoundaryConditions(Lau1DCubicSpline::LauSplineBoundaryType magLeftBound, - Lau1DCubicSpline::LauSplineBoundaryType magRightBound, - Lau1DCubicSpline::LauSplineBoundaryType phaseLeftBound, - Lau1DCubicSpline::LauSplineBoundaryType phaseRightBound, + void setBoundaryConditions(Lau1DCubicSpline::BoundaryType magLeftBound, + Lau1DCubicSpline::BoundaryType magRightBound, + Lau1DCubicSpline::BoundaryType phaseLeftBound, + Lau1DCubicSpline::BoundaryType phaseRightBound, Double_t magLeftGrad = 0.0, Double_t magRightGrad = 0.0, Double_t phaseLeftGrad = 0.0, Double_t phaseRightGrad = 0.0) { this->setSplineBoundaryConditions(magLeftBound,magRightBound,phaseLeftBound,phaseRightBound,magLeftGrad,magRightGrad,phaseLeftGrad,phaseRightGrad); } //! Get the resonance model type /*! \return the resonance model type */ virtual LauAbsResonance::LauResonanceModel getResonanceModel() const {return LauAbsResonance::MIPW_MagPhase;} protected: //! Evaluate the amplitude at the given point from the splines /*! \param [in] mass appropriate invariant mass for the resonance */ virtual void evaluateAmplitude(const Double_t mass); //! Method to create the parameter objects for the given knot /*! \param [in] iKnot the index of the knot */ virtual void createAmpParameters(const UInt_t iKnot); private: ClassDef(LauModIndPartWaveMagPhase,0) // model independent partial wave }; #endif diff --git a/inc/LauModIndPartWaveRealImag.hh b/inc/LauModIndPartWaveRealImag.hh index e8da06d..75e0814 100644 --- a/inc/LauModIndPartWaveRealImag.hh +++ b/inc/LauModIndPartWaveRealImag.hh @@ -1,125 +1,122 @@ /* Copyright 2015 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 LauModIndPartWaveRealImag.hh \brief File containing declaration of LauModIndPartWaveRealImag class. */ /*! \class LauModIndPartWaveRealImag \brief Class for defining a model independent partial wave component where the amplitudes are parameterised in terms of real and imaginary parts This model uses splines to produce a partial wave from the values of the real and imaginary parts of the amplitude at a series of points in the phase space. The values at each point can be floated in the fit. */ #ifndef LAU_MODINDPARTWAVE_REALIMAG #define LAU_MODINDPARTWAVE_REALIMAG #include "LauComplex.hh" #include "LauAbsModIndPartWave.hh" class TSpline3; class LauModIndPartWaveRealImag : public LauAbsModIndPartWave { public: //! Constructor /*! \param [in] resInfo the object containing information on the resonance name, mass, width, spin, charge, etc. \param [in] resPairAmpInt the number of the daughter not produced by the resonance \param [in] daughters the daughter particles */ LauModIndPartWaveRealImag(LauResonanceInfo* resInfo, Int_t resPairAmpInt, const LauDaughters* daughters); - //! Destructor - virtual ~LauModIndPartWaveRealImag(); - //! Set the values of the two real parameters that define the amplitude at a given knot /*! \param [in] knot the knot to be updated \param [in] realVal the value of the real part at the knot \param [in] imagVal the value of the imaginary part at the knot \param [in] fixRealPart whether the real part should be fixed \param [in] fixImagPart whether the imaginary part should be fixed */ virtual void setKnotAmp(const UInt_t knot, const Double_t realVal, const Double_t imagVal, const Bool_t fixRealPart, const Bool_t fixImagPart); //! Method to set the type of interpolation used for the splines /*! \param [in] realType the type of interpolation for the real part spline \param [in] imagType the type of interpolation for the imaginary part spline */ - void setType(Lau1DCubicSpline::LauSplineType realType, Lau1DCubicSpline::LauSplineType imagType) { + void setType(Lau1DCubicSpline::SplineType realType, Lau1DCubicSpline::SplineType imagType) { this->setSplineType(realType,imagType); } //! Method to set the boundary conditions of the splines /*! \param [in] realLeftBound the type of boundary condition for the left edge of the real part spline \param [in] realRightBound the type of boundary condition for the right edge of the real part spline \param [in] imagLeftBound the type of boundary condition for the left edge of the imaginary part spline \param [in] imagRightBound the type of boundary condition for the right edge of the imaginary part spline \param [in] realLeftGrad the gradient at the left edge of the real part spline if clamped \param [in] realRightGrad the gradient at the right edge of the real part spline if clamped \param [in] imagLeftGrad the gradient at the left edge of the imaginary part spline if clamped \param [in] imagRightGrad the gradient at the right edge of the imaginary part spline if clamped */ - void setBoundaryConditions(Lau1DCubicSpline::LauSplineBoundaryType realLeftBound, - Lau1DCubicSpline::LauSplineBoundaryType realRightBound, - Lau1DCubicSpline::LauSplineBoundaryType imagLeftBound, - Lau1DCubicSpline::LauSplineBoundaryType imagRightBound, + void setBoundaryConditions(Lau1DCubicSpline::BoundaryType realLeftBound, + Lau1DCubicSpline::BoundaryType realRightBound, + Lau1DCubicSpline::BoundaryType imagLeftBound, + Lau1DCubicSpline::BoundaryType imagRightBound, Double_t realLeftGrad = 0.0, Double_t realRightGrad = 0.0, Double_t imagLeftGrad = 0.0, Double_t imagRightGrad = 0.0) { this->setSplineBoundaryConditions(realLeftBound,realRightBound,imagLeftBound,imagRightBound,realLeftGrad,realRightGrad,imagLeftGrad,imagRightGrad); } //! Get the resonance model type /*! \return the resonance model type */ virtual LauAbsResonance::LauResonanceModel getResonanceModel() const {return LauAbsResonance::MIPW_RealImag;} protected: //! Evaluate the amplitude at the given point from the splines /*! \param [in] mass appropriate invariant mass for the resonance */ virtual void evaluateAmplitude(const Double_t mass); //! Method to create the parameter objects for the given knot /*! \param [in] iKnot the index of the knot */ virtual void createAmpParameters(const UInt_t iKnot); private: ClassDef(LauModIndPartWaveRealImag,0) // model independent partial wave }; #endif diff --git a/inc/LauSplineDecayTimeEfficiency.hh b/inc/LauSplineDecayTimeEfficiency.hh index a0dbe57..f528566 100644 --- a/inc/LauSplineDecayTimeEfficiency.hh +++ b/inc/LauSplineDecayTimeEfficiency.hh @@ -1,424 +1,424 @@ /* 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 "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; //! 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 PosDef!" << std::endl; + 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(); } #endif diff --git a/src/Lau1DCubicSpline.cc b/src/Lau1DCubicSpline.cc index 61bcf4a..aecb5db 100644 --- a/src/Lau1DCubicSpline.cc +++ b/src/Lau1DCubicSpline.cc @@ -1,567 +1,587 @@ /* Copyright 2015 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 Lau1DCubicSpline.cc \brief File containing implementation of Lau1DCubicSpline class. */ #include #include #include #include #include #include "Lau1DCubicSpline.hh" #include "LauAbsRValue.hh" #include "LauParameter.hh" ClassImp(Lau1DCubicSpline) Lau1DCubicSpline::Lau1DCubicSpline(const std::vector& xs, const std::vector& ys, - LauSplineType type, - LauSplineBoundaryType leftBound, - LauSplineBoundaryType rightBound, - Double_t dydx0, Double_t dydxn) : - nKnots_(xs.size()), - x_(xs), - y_(ys), - type_(type), - leftBound_(leftBound), - rightBound_(rightBound), - dydx0_(dydx0), - dydxn_(dydxn) + const SplineType type, + const BoundaryType leftBound, + const BoundaryType rightBound, + const Double_t dydx0, const Double_t dydxn) : + nKnots_{xs.size()}, + x_{xs}, + y_{ys}, + type_{type}, + leftBound_{leftBound}, + rightBound_{rightBound}, + dydx0_{dydx0}, + dydxn_{dydxn} { init(); } -Double_t Lau1DCubicSpline::evaluate(Double_t x) const +Double_t Lau1DCubicSpline::evaluate(const Double_t x) const { // do not attempt to extrapolate the spline if( xx_[nKnots_-1] ) { std::cerr << "WARNING in Lau1DCubicSpline::evaluate : function is only defined between " << x_[0] << " and " << x_[nKnots_-1] << std::endl; std::cerr << " value at " << x << " returned as 0" << std::endl; return 0.; } // first determine which 'cell' of the spline x is in // cell i runs from knot i to knot i+1 - Int_t cell(0); + std::size_t cell{0}; while( x > x_[cell+1] ) { ++cell; } // obtain x- and y-values of the neighbouring knots - Double_t xLow = x_[cell]; - Double_t xHigh = x_[cell+1]; - Double_t yLow = y_[cell]; - Double_t yHigh = y_[cell+1]; + const Double_t xLow { x_[cell] }; + const Double_t xHigh { x_[cell+1] }; + const Double_t yLow { y_[cell] }; + const Double_t yHigh { y_[cell+1] }; - if(type_ == Lau1DCubicSpline::LinearInterpolation) { + if(type_ == SplineType::LinearInterpolation) { return yHigh*(x-xLow)/(xHigh-xLow) + yLow*(xHigh-x)/(xHigh-xLow); } // obtain t, the normalised x-coordinate within the cell, // and the coefficients a and b, which are defined in cell i as: // // a_i = k_i *(x_i+1 - x_i) - (y_i+1 - y_i), // b_i = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i) // // where k_i is (by construction) the first derivative at knot i - Double_t t = (x - xLow) / (xHigh - xLow); - Double_t a = dydx_[cell] * (xHigh - xLow) - (yHigh - yLow); - Double_t b = -1.*dydx_[cell+1] * (xHigh - xLow) + (yHigh - yLow); + const Double_t t { (x - xLow) / (xHigh - xLow) }; + const Double_t a { dydx_[cell] * (xHigh - xLow) - (yHigh - yLow) }; + const Double_t b { -1.*dydx_[cell+1] * (xHigh - xLow) + (yHigh - yLow) }; - Double_t retVal = (1 - t) * yLow + t * yHigh + t * (1 - t) * ( a * (1 - t) + b * t ); + const Double_t retVal { (1 - t) * yLow + t * yHigh + t * (1 - t) * ( a * (1 - t) + b * t ) }; return retVal; } void Lau1DCubicSpline::updateYValues(const std::vector& ys) { y_ = ys; this->calcDerivatives(); } void Lau1DCubicSpline::updateYValues(const std::vector& ys) { - for (std::size_t i=0; iunblindValue(); } this->calcDerivatives(); } void Lau1DCubicSpline::updateYValues(const std::vector& ys) { - for (std::size_t i=0; iunblindValue(); } this->calcDerivatives(); } -void Lau1DCubicSpline::updateType(LauSplineType type) +void Lau1DCubicSpline::updateType(const SplineType type) { if(type_ != type) { type_ = type; this->calcDerivatives(); } } -void Lau1DCubicSpline::updateBoundaryConditions(LauSplineBoundaryType leftBound, - LauSplineBoundaryType rightBound, - Double_t dydx0, Double_t dydxn) +void Lau1DCubicSpline::updateBoundaryConditions(const BoundaryType leftBound, + const BoundaryType rightBound, + const Double_t dydx0, const Double_t dydxn) { bool updateDerivatives{false}; if(leftBound_ != leftBound || rightBound_ != rightBound) { leftBound_ = leftBound; rightBound_ = rightBound; updateDerivatives = true; } if(dydx0_ != dydx0) { dydx0_ = dydx0; - if(leftBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = true; + if(leftBound_ == BoundaryType::Clamped) updateDerivatives = true; } if(dydxn_ != dydxn) { dydxn_ = dydxn; - if(rightBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = true; + if(rightBound_ == BoundaryType::Clamped) updateDerivatives = true; } if(updateDerivatives) { this->calcDerivatives(); } } std::array Lau1DCubicSpline::getCoefficients(const std::size_t segIndex, const bool normalise) const { std::array result = {0.,0.,0.,0.}; if(segIndex >= nKnots_-1) { std::cerr << "ERROR in Lau1DCubicSpline::getCoefficients requested for too high a knot value" << std::endl; return result; } Double_t xL = x_[segIndex] , xH = x_[segIndex+1]; Double_t yL = y_[segIndex] , yH = y_[segIndex+1]; Double_t h = xH-xL; //This number comes up a lot switch(type_) { - case Lau1DCubicSpline::StandardSpline: - case Lau1DCubicSpline::AkimaSpline: + case SplineType::StandardSpline: + case SplineType::AkimaSpline: { Double_t kL = dydx_[segIndex], kH = dydx_[segIndex+1]; //a and b based on definitions from https://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline Double_t a = kL*h-(yH-yL); Double_t b =-kH*h+(yH-yL); Double_t denom = -h*h*h;//The terms have a common demoninator result[0] = -b*xL*xL*xH + a*xL*xH*xH + h*h*(xL*yH - xH*yL); result[1] = -a*xH*(2*xL+xH) + b*xL*(xL+2*xH) + h*h*(yL-yH); result[2] = -b*(2*xL+xH) + a*(xL+2*xH); result[3] = -a+b; for(auto& res : result){res /= denom;} break; } -/* case Lau1DCubicSpline::AkimaSpline: // Double check the Akima description of splines (in evaluate) right now they're the same except for the first derivatives +/* case SplineType::AkimaSpline: // Double check the Akima description of splines (in evaluate) right now they're the same except for the first derivatives { //using fomulae from https://asmquantmacro.com/2015/09/01/akima-spline-interpolation-in-excel/ std::function m = [&](Int_t j) //formula to get the straight line gradient { if(j < 0){return 2*m(j+1)-m(j+2);} if(j >= nKnots_){return 2*m(j-1)-m(j-2);} return (y_[j+1]-y_[j]) / (x_[j+1]-x_[j]); }; auto t = [&](Int_t j) { Double_t res = 0.; //originally res was called 't' but that produced a shadow warning Double_t denom = TMath::Abs(m(j+1)-m(j)) + TMath::Abs(m(j-1)-m(j-2)); if(denom == 0){res = (m(j)-m(j-1))/2;} //Special case for when denom = 0 else { res = TMath::Abs(m(j+1)-m(j))*m(j-1) + TMath::Abs(m(j-1)-m(j-2))*m(j); res /= denom; } return res; }; //These are the p's to get the spline in the form p_k(x-xL)^k Double_t pDenom = x_[segIndex+1]/x_[segIndex]; //a denominator used for p[2] and p[3] std::array p = {y_[segIndex],t(segIndex),0.,0.}; //we'll do the last 2 below p[2] = 3*m(segIndex)-2*t(segIndex)-t(segIndex+1); p[2]/= pDenom; p[3] = t(segIndex)+t(segIndex+1)-2*m(segIndex); p[3]/= pDenom*pDenom; //Now finally rearranging the p's into the desired results result[0] = p[0]-p[1]*xL+p[2]*xL*xL-p[3]*xL*xL*xL; result[1] = p[1]-2*p[2]*xL+3*p[3]*xL*xL; result[2] = p[2]-3*p[3]*xL; result[3] = p[3]; break; }*/ - case Lau1DCubicSpline::LinearInterpolation: + case SplineType::LinearInterpolation: { result[0] = xH*yL-xL*yH; result[1] = yH-yL; for(auto& res : result){res /= h;} break; } } - if(normalise) - { - Double_t integral = this->integral(); + if (normalise) { + const Double_t integral { this->integral() }; for(auto& res : result){res /= integral;} } return result; } Double_t Lau1DCubicSpline::integral() const { Double_t integral{0.0}; for(std::size_t iKnot{0}; iKnot < nKnots_ -1; ++iKnot) { const Double_t minAbs { x_[iKnot] }; const Double_t maxAbs { x_[iKnot+1] }; const std::array coeffs { this->getCoefficients(iKnot, false) }; auto integralFunc = [&coeffs](Double_t x){return coeffs[0]*x + coeffs[1]*x*x/2.0 + coeffs[2]*x*x*x/3.0 + coeffs[3]*x*x*x*x/4.0;}; integral += integralFunc(maxAbs); integral -= integralFunc(minAbs); } return integral; } TF1* Lau1DCubicSpline::makeTF1(const bool normalise) const { TString functionString{""}; //make a long piecewise construction of all the spline pieces for(std::size_t segIndex{0}; segIndex < nKnots_-1; ++segIndex) { const std::array coeffs { this->getCoefficients(segIndex,normalise) }; functionString += Form("(x>%f && x<= %f)*",x_[segIndex],x_[segIndex+1]);//get the bounds of this piece functionString += Form("(%f + %f*x + %f*x^2 + %f*x^3)",coeffs[0],coeffs[1],coeffs[2],coeffs[3]); if(segIndex < nKnots_ -2){functionString += " + \n";}//add to all lines except the last } TF1* func = new TF1("SplineFunction", functionString, x_.front(), x_.back()); return func; } -TF1* Lau1DCubicSpline::normaliseToTH1(TH1* hist, LauOutputLevel printLevel) const +TF1* Lau1DCubicSpline::normaliseToTH1(TH1* hist, const LauOutputLevel printLevel) const { //first define the fit function auto fitf = [this](Double_t* x, Double_t* par) {//there is only 1 x (the abscissa) and 1 par (a scaling of the entire thing) return this->evaluate( x[0] ) * par[0]; }; //Make the function TF1* FittedFunc = new TF1("FittedFunc",fitf,x_.front(),x_.back(),1); //Set the parameter name FittedFunc -> SetParNames("Constant"); //Set the parameter limits and default value FittedFunc -> SetParLimits(0,0.,10.); FittedFunc -> SetParameter(0,1.); // Options to: // - N = do not store the graphics function // - 0 = do not plot the fit result TString fitOptions {"N 0"}; switch ( printLevel ) { case LauOutputLevel::Verbose : fitOptions += " V"; break; case LauOutputLevel::Quiet : fitOptions += " Q"; break; case LauOutputLevel::Standard : break; } hist->Fit( "FittedFunc", fitOptions ); return FittedFunc; } -TFitResultPtr Lau1DCubicSpline::fitToTH1(TH1* hist, LauOutputLevel printLevel, const bool doWL) +TFitResultPtr Lau1DCubicSpline::fitToTH1(TH1* hist, const LauOutputLevel printLevel, const bool doWL) { auto fitf = [this](Double_t* x, Double_t* par) { this -> updateYValues( std::vector(par, par + nKnots_) ); return this -> evaluate( x[0] ); }; TF1 FittedFunc("FittedFunc",fitf,x_.front(),x_.back(),nKnots_); const Double_t knotMax { hist->GetMaximum() * 1.5 }; for(std::size_t knot{0}; knot <= nKnots_ ; ++knot) { FittedFunc.SetParName(knot, Form("Knot%lu",knot)); FittedFunc.SetParLimits(knot, 0., knotMax); FittedFunc.SetParameter(knot, y_[knot]); } // Options to: // - N = do not store the graphics function // - 0 = do not plot the fit result // - S = save the result of the fit // - WL= do a weighted likelihood fit (not chi2) TString fitOptions {"N 0 S"}; if(doWL){fitOptions += " WL";} switch ( printLevel ) { case LauOutputLevel::Verbose : fitOptions += " V"; break; case LauOutputLevel::Quiet : fitOptions += " Q"; break; case LauOutputLevel::Standard : break; } return hist->Fit( "FittedFunc", fitOptions ); } void Lau1DCubicSpline::init() { if( y_.size() != x_.size()) { std::cerr << "ERROR in Lau1DCubicSpline::init : The number of y-values given does not match the number of x-values" << std::endl; std::cerr << " Found " << y_.size() << ", expected " << x_.size() << std::endl; gSystem->Exit(EXIT_FAILURE); } if( nKnots_ < 3 ) { std::cerr << "ERROR in Lau1DCubicSpline::init : The number of knots is too small" << std::endl; std::cerr << " Found " << nKnots_ << ", expected at least 3 (to have at least 1 internal knot)" << std::endl; gSystem->Exit(EXIT_FAILURE); } dydx_.assign(nKnots_,0.0); a_.assign(nKnots_,0.0); b_.assign(nKnots_,0.0); c_.assign(nKnots_,0.0); d_.assign(nKnots_,0.0); this->calcDerivatives(); } void Lau1DCubicSpline::calcDerivatives() { switch ( type_ ) { - case Lau1DCubicSpline::StandardSpline : + case SplineType::StandardSpline : this->calcDerivativesStandard(); break; - case Lau1DCubicSpline::AkimaSpline : + case SplineType::AkimaSpline : this->calcDerivativesAkima(); break; - case Lau1DCubicSpline::LinearInterpolation : + case SplineType::LinearInterpolation : //derivatives not needed for linear interpolation break; } } void Lau1DCubicSpline::calcDerivativesStandard() { // derivatives are determined such that the second derivative is continuous at internal knots // derivatives, k_i, are the solutions to a set of linear equations of the form: // a_i * k_i-1 + b_i * k+i + c_i * k_i+1 = d_i with a_0 = 0, c_n-1 = 0 // this is solved using the tridiagonal matrix algorithm as on en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm // first and last equations give boundary conditions // - for natural boundary, require f''(x) = 0 at end knot // - for 'not a knot' boundary, require f'''(x) continuous at second knot // - for clamped boundary, require predefined value of f'(x) at end knot // non-zero values of a_0 and c_n-1 would give cyclic boundary conditions a_[0] = 0.; c_[nKnots_-1] = 0.; // set left boundary condition - if(leftBound_ == Lau1DCubicSpline::Natural) { - b_[0] = 2./(x_[1]-x_[0]); - c_[0] = 1./(x_[1]-x_[0]); - d_[0] = 3.*(y_[1]-y_[0])/((x_[1]-x_[0])*(x_[1]-x_[0])); - - } else if(leftBound_ == Lau1DCubicSpline::NotAKnot) { - // define the width, h, and the 'slope', delta, of the first cell - Double_t h1(x_[1]-x_[0]), h2(x_[2]-x_[0]); - Double_t delta1((y_[1]-y_[0])/h1), delta2((y_[2]-y_[1])/h2); - - // these coefficients can be determined by requiring f'''_0(x_1) = f'''_1(x_1) - // the requirement f''_0(x_1) = f''_1(x_1) has been used to remove the dependence on k_2 - b_[0] = h2; - c_[0] = h1+h2; - d_[0] = delta1*(2.*h2*h2 + 3.*h1*h2)/(h1+h2) + delta2*5.*h1*h1/(h1+h2); - - } else { //Clamped - b_[0] = 1.; - c_[0] = 0.; - d_[0] = dydx0_; + switch ( leftBound_ ) { + case BoundaryType::Natural : + { + b_[0] = 2./(x_[1]-x_[0]); + c_[0] = 1./(x_[1]-x_[0]); + d_[0] = 3.*(y_[1]-y_[0])/((x_[1]-x_[0])*(x_[1]-x_[0])); + break; + } + + case BoundaryType::NotAKnot : + { + // define the width, h, and the 'slope', delta, of the first cell + const Double_t h1{x_[1]-x_[0]}; + const Double_t h2{x_[2]-x_[0]}; + const Double_t delta1{(y_[1]-y_[0])/h1}; + const Double_t delta2{(y_[2]-y_[1])/h2}; + // these coefficients can be determined by requiring f'''_0(x_1) = f'''_1(x_1) + // the requirement f''_0(x_1) = f''_1(x_1) has been used to remove the dependence on k_2 + b_[0] = h2; + c_[0] = h1+h2; + d_[0] = delta1*(2.*h2*h2 + 3.*h1*h2)/(h1+h2) + delta2*5.*h1*h1/(h1+h2); + break; + } + + case BoundaryType::Clamped : + { + b_[0] = 1.; + c_[0] = 0.; + d_[0] = dydx0_; + break; + } } // set right boundary condition - if(rightBound_ == Lau1DCubicSpline::Natural) { - a_[nKnots_-1] = 1./(x_[nKnots_-1]-x_[nKnots_-2]); - b_[nKnots_-1] = 2./(x_[nKnots_-1]-x_[nKnots_-2]); - d_[nKnots_-1] = 3.*(y_[nKnots_-1]-y_[nKnots_-2])/((x_[nKnots_-1]-x_[nKnots_-2])*(x_[nKnots_-1]-x_[nKnots_-2])); - - } else if(rightBound_ == Lau1DCubicSpline::NotAKnot) { - // define the width, h, and the 'slope', delta, of the last cell - Double_t hnm1(x_[nKnots_-1]-x_[nKnots_-2]), hnm2(x_[nKnots_-2]-x_[nKnots_-3]); - Double_t deltanm1((y_[nKnots_-1]-y_[nKnots_-2])/hnm1), deltanm2((y_[nKnots_-2]-y_[nKnots_-3])/hnm2); - - // these coefficients can be determined by requiring f'''_n-3(x_n-2) = f'''_n-2(x_n-2) - // the requirement f''_n-3(x_n-2) = f''_n-2(x_n-2) has been used to remove - // the dependence on k_n-3 - a_[nKnots_-1] = hnm2 + hnm1; - b_[nKnots_-1] = hnm1; - d_[nKnots_-1] = deltanm2*hnm1*hnm1/(hnm2+hnm1) + deltanm1*(2.*hnm2*hnm2 + 3.*hnm2*hnm1)/(hnm2+hnm1); - - } else { //Clamped - a_[nKnots_-1] = 0.; - b_[nKnots_-1] = 1.; - d_[nKnots_-1] = dydxn_; + switch ( rightBound_ ) { + case BoundaryType::Natural : + { + a_[nKnots_-1] = 1./(x_[nKnots_-1]-x_[nKnots_-2]); + b_[nKnots_-1] = 2./(x_[nKnots_-1]-x_[nKnots_-2]); + d_[nKnots_-1] = 3.*(y_[nKnots_-1]-y_[nKnots_-2])/((x_[nKnots_-1]-x_[nKnots_-2])*(x_[nKnots_-1]-x_[nKnots_-2])); + break; + } + case BoundaryType::NotAKnot : + { + // define the width, h, and the 'slope', delta, of the last cell + const Double_t hnm1{x_[nKnots_-1]-x_[nKnots_-2]}; + const Double_t hnm2(x_[nKnots_-2]-x_[nKnots_-3]); + const Double_t deltanm1{(y_[nKnots_-1]-y_[nKnots_-2])/hnm1}; + const Double_t deltanm2{(y_[nKnots_-2]-y_[nKnots_-3])/hnm2}; + // these coefficients can be determined by requiring f'''_n-3(x_n-2) = f'''_n-2(x_n-2) + // the requirement f''_n-3(x_n-2) = f''_n-2(x_n-2) has been used to remove + // the dependence on k_n-3 + a_[nKnots_-1] = hnm2 + hnm1; + b_[nKnots_-1] = hnm1; + d_[nKnots_-1] = deltanm2*hnm1*hnm1/(hnm2+hnm1) + deltanm1*(2.*hnm2*hnm2 + 3.*hnm2*hnm1)/(hnm2+hnm1); + break; + } + + case BoundaryType::Clamped : + { + a_[nKnots_-1] = 0.; + b_[nKnots_-1] = 1.; + d_[nKnots_-1] = dydxn_; + break; + } } // the remaining equations ensure that f_i-1''(x_i) = f''_i(x_i) for all internal knots for(std::size_t i{1}; i(nKnots_-2)}; i>=0; --i) { dydx_[i] = d_[i] - c_[i]*dydx_[i+1]; } } void Lau1DCubicSpline::calcDerivativesAkima() { //derivatives are calculated according to the Akima method // J.ACM vol. 17 no. 4 pp 589-602 - Double_t am1(0.), an(0.), anp1(0.); + Double_t am1{0.0}, an{0.0}, anp1{0.0}; // a[i] is the slope of the segment from point i-1 to point i // // n.b. segment 0 is before point 0 and segment n is after point n-1 // internal segments are numbered 1 - n-1 for(std::size_t i{1}; i #include #include "LauConstants.hh" #include "LauKinematics.hh" #include "LauAbsModIndPartWave.hh" #include "LauResonanceInfo.hh" ClassImp(LauAbsModIndPartWave) LauAbsModIndPartWave::LauAbsModIndPartWave(LauResonanceInfo* resInfo, Int_t resPairAmpInt, const LauDaughters* daughters) : - LauAbsResonance(resInfo, resPairAmpInt, daughters), - nKnots_(0), - spline1_(0), - spline2_(0), - type1_(Lau1DCubicSpline::StandardSpline), - type2_(Lau1DCubicSpline::StandardSpline), - leftBound1_(Lau1DCubicSpline::NotAKnot), - rightBound1_(Lau1DCubicSpline::NotAKnot), - leftBound2_(Lau1DCubicSpline::NotAKnot), - rightBound2_(Lau1DCubicSpline::NotAKnot), - leftGrad1_(0.), - rightGrad1_(0.), - leftGrad2_(0.), - rightGrad2_(0.), - secondStage_(kFALSE) + LauAbsResonance(resInfo, resPairAmpInt, daughters) { } -LauAbsModIndPartWave::~LauAbsModIndPartWave() -{ - delete spline1_; - delete spline2_; -} - void LauAbsModIndPartWave::floatKnotsSecondStage(const Bool_t secondStage) { secondStage_ = secondStage; // if the parameters have not yet been created we can now just return if ( amp1Pars_.size() != nKnots_ ) { return; } // otherwise we need to toggle their second-stage parameter for ( UInt_t i(0); i < nKnots_; ++i ) { amp1Pars_[i]->secondStage(secondStage_); amp2Pars_[i]->secondStage(secondStage_); } } std::set LauAbsModIndPartWave::checkKnots(const std::set& masses) { std::set knots = masses; const std::set::const_iterator first = knots.begin(); const std::set::const_reverse_iterator last = knots.rbegin(); const Double_t lower_limit = this->getMassDaug1() + this->getMassDaug2(); const Double_t upper_limit = this->getMassParent() - this->getMassBachelor(); // check whether we have been given knots at unphysical masses if ( *first < lower_limit ) { std::cerr << "WARNING in LauAbsModIndPartWave::checkKnots : Knot found at mass " << *first << " is below the lower kinematic limit." << std::endl; std::cerr << " : Lower kinematic limit is at mass " << lower_limit << std::endl; std::cerr << " : Aborting definition of knot positions." << std::endl; knots.clear(); return knots; } if ( *last > upper_limit ) { std::cerr << "WARNING in LauAbsModIndPartWave::checkKnots : Knot found at mass " << *last << " is above the upper kinematic limit." << std::endl; std::cerr << " : Upper kinematic limit is at mass " << upper_limit << std::endl; std::cerr << " : Aborting definition of knot positions." << std::endl; knots.clear(); return knots; } // check if we have knots at each extreme - if not, add them in if ( (*first) != lower_limit ) { knots.insert( lower_limit ); } if ( (*last) != upper_limit ) { knots.insert( upper_limit ); } return knots; } void LauAbsModIndPartWave::defineKnots(const std::set& masses) { if ( ! masses_.empty() ) { std::cerr << "WARNING in LauAbsModIndPartWave::defineKnots : Knot positions have already been defined, not making any changes." << std::endl; return; } const std::set knots = this->checkKnots( masses ); nKnots_ = knots.size(); if ( nKnots_ == 0 ) { return; } masses_.reserve(nKnots_); amp1Vals_.reserve(nKnots_); amp2Vals_.reserve(nKnots_); amp1Pars_.reserve(nKnots_); amp2Pars_.reserve(nKnots_); UInt_t counter(0); for ( std::set::const_iterator iter = knots.begin(); iter != knots.end(); ++iter ) { masses_.push_back( *iter ); amp1Vals_.push_back(1.0); amp2Vals_.push_back(1.0); this->createAmpParameters(counter); ++counter; } for ( std::vector::const_iterator iter = masses_.begin(); iter != masses_.end(); ++iter ) { std::cout << "INFO in LauAbsModIndPartWave::defineKnots : Knot added to resonance " << this->getResonanceName() << " at mass " << *iter << std::endl; } } void LauAbsModIndPartWave::initialise() { - if ( spline1_ != 0 ) { - delete spline1_; - spline1_ = 0; - } - if ( spline2_ != 0 ) { - delete spline2_; - spline2_ = 0; - } - for ( UInt_t i(0); i < nKnots_; ++i ) { amp1Vals_[i] = amp1Pars_[i]->unblindValue(); amp2Vals_[i] = amp2Pars_[i]->unblindValue(); } - spline1_ = new Lau1DCubicSpline(masses_, amp1Vals_, type1_, leftBound1_, rightBound1_, leftGrad1_, rightGrad1_); - spline2_ = new Lau1DCubicSpline(masses_, amp2Vals_, type2_, leftBound2_, rightBound2_, leftGrad2_, rightGrad2_); + spline1_ = std::make_unique(masses_, amp1Vals_, type1_, leftBound1_, rightBound1_, leftGrad1_, rightGrad1_); + spline2_ = std::make_unique(masses_, amp2Vals_, type2_, leftBound2_, rightBound2_, leftGrad2_, rightGrad2_); } LauComplex LauAbsModIndPartWave::resAmp(Double_t mass, Double_t spinTerm) { amp_.zero(); Bool_t paramChanged1(kFALSE), paramChanged2(kFALSE); for ( UInt_t i(0); i < nKnots_; ++i ) { if ( !amp1Pars_[i]->fixed() && amp1Pars_[i]->unblindValue() != amp1Vals_[i] ) { paramChanged1 = kTRUE; amp1Vals_[i] = amp1Pars_[i]->unblindValue(); } if ( !amp2Pars_[i]->fixed() && amp2Pars_[i]->unblindValue() != amp2Vals_[i] ) { paramChanged2 = kTRUE; amp2Vals_[i] = amp2Pars_[i]->unblindValue(); } } - if ( spline1_ == 0 || spline2_ == 0) { - std::cerr << "ERROR in LauAbsModIndPartWave::resAmp : One of the splines is null" << std::endl; + if ( spline1_ == nullptr || spline2_ == nullptr) { + std::cerr << "ERROR in LauAbsModIndPartWave::resAmp : One or both of the splines is null" << std::endl; return amp_; } if ( paramChanged1 ) { spline1_->updateYValues(amp1Vals_); } if ( paramChanged2 ) { spline2_->updateYValues(amp2Vals_); } this->evaluateAmplitude( mass ); amp_.rescale( spinTerm ); return amp_; } -void LauAbsModIndPartWave::setSplineType(Lau1DCubicSpline::LauSplineType type1, Lau1DCubicSpline::LauSplineType type2) +void LauAbsModIndPartWave::setSplineType(Lau1DCubicSpline::SplineType type1, Lau1DCubicSpline::SplineType type2) { type1_ = type1; type2_ = type2; } -void LauAbsModIndPartWave::setSplineBoundaryConditions(Lau1DCubicSpline::LauSplineBoundaryType leftBound1, - Lau1DCubicSpline::LauSplineBoundaryType rightBound1, - Lau1DCubicSpline::LauSplineBoundaryType leftBound2, - Lau1DCubicSpline::LauSplineBoundaryType rightBound2, +void LauAbsModIndPartWave::setSplineBoundaryConditions(Lau1DCubicSpline::BoundaryType leftBound1, + Lau1DCubicSpline::BoundaryType rightBound1, + Lau1DCubicSpline::BoundaryType leftBound2, + Lau1DCubicSpline::BoundaryType rightBound2, Double_t leftGrad1, - Double_t rightGrad1, + Double_t rightGrad1, Double_t leftGrad2, - Double_t rightGrad2) + Double_t rightGrad2) { leftBound1_ = leftBound1; rightBound1_ = rightBound1; leftBound2_ = leftBound2; rightBound2_ = rightBound2; leftGrad1_ = leftGrad1; rightGrad1_ = rightGrad1; leftGrad2_ = leftGrad2; rightGrad2_ = rightGrad2; } const std::vector& LauAbsModIndPartWave::getFloatingParameters() { this->clearFloatingParameters(); for ( UInt_t i(0); i < nKnots_; ++i ) { if ( !amp1Pars_[i]->fixed() ) { this->addFloatingParameter( amp1Pars_[i] ); } if ( !amp2Pars_[i]->fixed() ) { this->addFloatingParameter( amp2Pars_[i] ); } } return this->getParameters(); } diff --git a/src/LauEFKLLMRes.cc b/src/LauEFKLLMRes.cc index 1d9a7f8..8d8949e 100644 --- a/src/LauEFKLLMRes.cc +++ b/src/LauEFKLLMRes.cc @@ -1,175 +1,169 @@ /* Copyright 2015 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 LauEFKLLMRes.cc \brief File containing implementation of LauEFKLLMRes class. */ #include #include "Lau1DCubicSpline.hh" #include "LauKinematics.hh" #include "LauEFKLLMRes.hh" #include "LauResonanceInfo.hh" #include "LauTextFileParser.hh" ClassImp(LauEFKLLMRes); -Lau1DCubicSpline* LauEFKLLMRes::magSpline_ = 0; -Lau1DCubicSpline* LauEFKLLMRes::phaseSpline_ = 0; +std::unique_ptr LauEFKLLMRes::magSpline_; +std::unique_ptr LauEFKLLMRes::phaseSpline_; LauEFKLLMRes::LauEFKLLMRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) : - LauAbsResonance(resInfo, resPairAmpInt, daughters), - massFactor_(0) + LauAbsResonance(resInfo, resPairAmpInt, daughters) { - const Double_t massFactorVal = 0.; + const Double_t massFactorVal{0.0}; - const TString& parNameBase = this->getSanitisedName(); + const TString& parNameBase { this->getSanitisedName() }; - TString massFactorName(parNameBase); + TString massFactorName{parNameBase}; massFactorName += "_massFactor"; massFactor_ = resInfo->getExtraParameter( massFactorName ); - if ( massFactor_ == 0 ) { + if ( massFactor_ == nullptr ) { massFactor_ = new LauParameter( massFactorName, massFactorVal, -10.0, 10.0, kTRUE ); massFactor_->secondStage(kTRUE); resInfo->addExtraParameter( massFactor_ ); } } -LauEFKLLMRes::~LauEFKLLMRes() -{ -} - void LauEFKLLMRes::initialise() { } void LauEFKLLMRes::setResonanceParameter(const TString& name, const Double_t value) { if(name=="massFactor") { this->setMassFactor(value); std::cout << "INFO in LauEFKLLMRes::setResonanceParameter: Mass factor set to " << value << std::endl; } else { std::cerr << "WARNING in LauEFKLLMRes::setResonanceParameter: Parameter name not reconised." << std::endl; } } void LauEFKLLMRes::floatResonanceParameter(const TString& name) { if(name=="massFactor") { if ( massFactor_->fixed() ) { massFactor_->fixed( kFALSE ); this->addFloatingParameter( massFactor_ ); } else { std::cerr << "WARNING in LauEFKLLMRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl; } } else { std::cerr << "WARNING in LauEFKLLMRes::floatResonanceParameter: Parameter name not reconised." << std::endl; } } LauParameter* LauEFKLLMRes::getResonanceParameter(const TString& name) { if(name=="massFactor") { return massFactor_; } else { std::cerr << "WARNING in LauEFKLLMRes::getResonanceParameter: Parameter name not reconised." << std::endl; - return 0; + return nullptr; } } const std::vector& LauEFKLLMRes::getFloatingParameters() { this->clearFloatingParameters(); if ( ! this->fixMassFactor() ) { this->addFloatingParameter( massFactor_ ); } return this->getParameters(); } LauComplex LauEFKLLMRes::resAmp(Double_t mass, Double_t /*spinTerm*/) { - LauComplex amp(0.0, 0.0); + LauComplex amp{0.0, 0.0}; - if (magSpline_ == 0 || phaseSpline_ == 0) { - std::cerr << "ERROR in LauEFKLLMRes::resAmp : One of the splines is null." << std::endl; + if (magSpline_ == nullptr || phaseSpline_ == nullptr) { + std::cerr << "ERROR in LauEFKLLMRes::resAmp : One or both of the splines is null." << std::endl; return amp; } - const Double_t massSq = mass*mass; - const Double_t mag = magSpline_->evaluate(massSq); - const Double_t phase = TMath::DegToRad()*phaseSpline_->evaluate(massSq); - LauComplex ff(mag*TMath::Cos(phase), mag*TMath::Sin(phase)); + const Double_t massSq { mass * mass }; + const Double_t mag { magSpline_->evaluate(massSq) }; + const Double_t phase { TMath::DegToRad() * phaseSpline_->evaluate(massSq) }; + LauComplex ff{mag*TMath::Cos(phase), mag*TMath::Sin(phase)}; amp = ff.scale(TMath::Power(mass,this->getMassFactor())); return amp; } void LauEFKLLMRes::setupFormFactor(const TString& inputFile) { LauTextFileParser readFile(inputFile); readFile.processFile(); std::vector mSqVals; std::vector magVals; std::vector phaseVals; std::vector line; - line=readFile.getNextLine(); - while(!line.empty()) { - UInt_t length = line.size(); - if(length!=3) { + line = readFile.getNextLine(); + while ( ! line.empty() ) { + if ( line.size() != 3 ) { std::cerr << "ERROR in LauEFKLLMRes::setupFormFactor : Unexpected number of fields in text file, aborting reading of form-factor information." << std::endl; return; } - mSqVals.push_back( atof(line[0].c_str())); - magVals.push_back( atof(line[1].c_str())); - phaseVals.push_back(atof(line[2].c_str())); - line=readFile.getNextLine(); + + mSqVals.push_back( std::stod(line[0]) ); + magVals.push_back( std::stod(line[1]) ); + phaseVals.push_back( std::stod(line[2]) ); + + line = readFile.getNextLine(); } // Destroy any splines we already had defined but issue a warning just in case - if ( magSpline_ != 0 || phaseSpline_ != 0 ) { + if ( magSpline_ != nullptr || phaseSpline_ != nullptr ) { std::cerr << "WARNING in LauEFKLLMRes::setupFormFactor : Overwriting previous form-factor splines with newly read values." << std::endl; - delete magSpline_; - delete phaseSpline_; } - magSpline_ = new Lau1DCubicSpline(mSqVals, magVals, Lau1DCubicSpline::AkimaSpline); - phaseSpline_ = new Lau1DCubicSpline(mSqVals, phaseVals, Lau1DCubicSpline::AkimaSpline); + magSpline_ = std::make_unique(mSqVals, magVals, Lau1DCubicSpline::SplineType::AkimaSpline); + phaseSpline_ = std::make_unique(mSqVals, phaseVals, Lau1DCubicSpline::SplineType::AkimaSpline); } void LauEFKLLMRes::setMassFactor(const Double_t massFactor) { massFactor_->value( massFactor ); massFactor_->genValue( massFactor ); massFactor_->initValue( massFactor ); } diff --git a/src/LauModIndPartWaveMagPhase.cc b/src/LauModIndPartWaveMagPhase.cc index c9176eb..28959c9 100644 --- a/src/LauModIndPartWaveMagPhase.cc +++ b/src/LauModIndPartWaveMagPhase.cc @@ -1,141 +1,137 @@ /* Copyright 2015 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 LauModIndPartWaveMagPhase.cc \brief File containing implementation of LauModIndPartWaveMagPhase class. */ #include #include #include "Lau1DCubicSpline.hh" #include "LauConstants.hh" #include "LauKinematics.hh" #include "LauModIndPartWaveMagPhase.hh" #include "LauResonanceInfo.hh" ClassImp(LauModIndPartWaveMagPhase) LauModIndPartWaveMagPhase::LauModIndPartWaveMagPhase(LauResonanceInfo* resInfo, Int_t resPairAmpInt, const LauDaughters* daughters) : LauAbsModIndPartWave(resInfo, resPairAmpInt, daughters) { } -LauModIndPartWaveMagPhase::~LauModIndPartWaveMagPhase() -{ -} - void LauModIndPartWaveMagPhase::createAmpParameters(const UInt_t iKnot) { const TString& parNameBase = this->getSanitisedName(); const Bool_t secondStage = this->floatKnotsSecondStage(); std::vector& magnitudePars = this->getAmp1Pars(); std::vector& phasePars = this->getAmp2Pars(); TString magName(parNameBase); magName += "_A"; magName += iKnot; magnitudePars.push_back(this->getResInfo()->getExtraParameter( magName )); if( magnitudePars[iKnot] == 0) { magnitudePars[iKnot] = new LauParameter( magName, 1.0, 0.0, 10.0, kFALSE); magnitudePars[iKnot]->secondStage(secondStage); this->getResInfo()->addExtraParameter(magnitudePars[iKnot]); } TString phaseName(parNameBase); phaseName += "_d"; phaseName += iKnot; phasePars.push_back(this->getResInfo()->getExtraParameter( phaseName )); if( phasePars[iKnot] == 0) { phasePars[iKnot] = new LauParameter( phaseName, 1.0, -6.0*LauConstants::pi, 6.0*LauConstants::pi, kFALSE); phasePars[iKnot]->secondStage(secondStage); this->getResInfo()->addExtraParameter(phasePars[iKnot]); } } void LauModIndPartWaveMagPhase::setKnotAmp(const UInt_t knot, const Double_t magVal, const Double_t phaseVal, const Bool_t fixMagnitude, const Bool_t fixPhase) { const UInt_t nknots = this->nKnots(); if ( knot >= nknots ) { std::cerr << "WARNING in LauModIndPartWaveMagPhase::setKnotAmp : Index " << knot << " does not correspond to an existing knot in resonance " << this->getResonanceName() << std::endl; std::cerr << " : Index must be in range 0 to " << nknots-1 << std::endl; return; } const std::vector& masses = this->getMasses(); std::vector& magnitudes = this->getAmp1Vals(); std::vector& phases = this->getAmp2Vals(); std::vector& magnitudePars = this->getAmp1Pars(); std::vector& phasePars = this->getAmp2Pars(); magnitudes[knot] = magVal; magnitudePars[knot]->value(magVal); magnitudePars[knot]->genValue(magVal); magnitudePars[knot]->initValue(magVal); magnitudePars[knot]->fixed(fixMagnitude); phases[knot] = phaseVal; phasePars[knot]->value(phaseVal); phasePars[knot]->genValue(phaseVal); phasePars[knot]->initValue(phaseVal); phasePars[knot]->fixed(fixPhase); if ( knot == 0 ) { std::cout << "INFO in LauModIndPartWaveMagPhase::setKnotAmp : Knot updated in resonance " << this->getResonanceName() << " at lower kinematic limit (" << masses[knot] << ")" << std::endl; } else if ( knot == nknots-1 ) { std::cout << "INFO in LauModIndPartWaveMagPhase::setKnotAmp : Knot updated in resonance " << this->getResonanceName() << " at upper kinematic limit (" << masses[knot] << ")" << std::endl; } else { std::cout << "INFO in LauModIndPartWaveMagPhase::setKnotAmp : Knot updated in resonance " << this->getResonanceName() << " at mass " << masses[knot] << std::endl; } if ( fixMagnitude ) { std::cout << " : Magnitude fixed to " << magVal << std::endl; } else { std::cout << " : Magnitude set to " << magVal << std::endl; } if ( fixPhase ) { std::cout << " : Phase fixed to " << phaseVal << std::endl; } else { std::cout << " : Phase set to " << phaseVal << std::endl; } } void LauModIndPartWaveMagPhase::evaluateAmplitude(const Double_t mass) { - const Lau1DCubicSpline* splineMag = this->getSpline1(); - const Lau1DCubicSpline* splinePhase = this->getSpline2(); + const Lau1DCubicSpline& splineMag { this->getSpline1() }; + const Lau1DCubicSpline& splinePhase { this->getSpline2() }; - const Double_t mag = splineMag->evaluate(mass); - const Double_t phase = splinePhase->evaluate(mass); + const Double_t mag { splineMag.evaluate(mass) }; + const Double_t phase { splinePhase.evaluate(mass) }; this->setAmp(mag*TMath::Cos(phase), mag*TMath::Sin(phase)); } diff --git a/src/LauModIndPartWaveRealImag.cc b/src/LauModIndPartWaveRealImag.cc index 7e63ba0..7b923f7 100644 --- a/src/LauModIndPartWaveRealImag.cc +++ b/src/LauModIndPartWaveRealImag.cc @@ -1,141 +1,137 @@ /* Copyright 2015 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 LauModIndPartWaveRealImag.cc \brief File containing implementation of LauModIndPartWaveRealImag class. */ #include #include #include "Lau1DCubicSpline.hh" #include "LauConstants.hh" #include "LauKinematics.hh" #include "LauModIndPartWaveRealImag.hh" #include "LauResonanceInfo.hh" ClassImp(LauModIndPartWaveRealImag) LauModIndPartWaveRealImag::LauModIndPartWaveRealImag(LauResonanceInfo* resInfo, Int_t resPairAmpInt, const LauDaughters* daughters) : LauAbsModIndPartWave(resInfo, resPairAmpInt, daughters) { } -LauModIndPartWaveRealImag::~LauModIndPartWaveRealImag() -{ -} - void LauModIndPartWaveRealImag::createAmpParameters(const UInt_t iKnot) { const TString& parNameBase = this->getSanitisedName(); const Bool_t secondStage = this->floatKnotsSecondStage(); std::vector& rePars = this->getAmp1Pars(); std::vector& imPars = this->getAmp2Pars(); TString reName(parNameBase); reName += "_Re"; reName += iKnot; rePars.push_back(this->getResInfo()->getExtraParameter( reName )); if( rePars[iKnot] == 0) { rePars[iKnot] = new LauParameter( reName, 1.0, -10.0, 10.0, kFALSE); rePars[iKnot]->secondStage(secondStage); this->getResInfo()->addExtraParameter(rePars[iKnot]); } TString imName(parNameBase); imName += "_Im"; imName += iKnot; imPars.push_back(this->getResInfo()->getExtraParameter( imName )); if( imPars[iKnot] == 0) { imPars[iKnot] = new LauParameter( imName, 1.0, -10.0, 10.0, kFALSE); imPars[iKnot]->secondStage(secondStage); this->getResInfo()->addExtraParameter(imPars[iKnot]); } } void LauModIndPartWaveRealImag::setKnotAmp(const UInt_t knot, const Double_t realVal, const Double_t imagVal, const Bool_t fixRealPart, const Bool_t fixImagPart) { const UInt_t nknots = this->nKnots(); if ( knot >= nknots ) { std::cerr << "WARNING in LauModIndPartWaveRealImag::setKnotAmp : Index " << knot << " does not correspond to an existing knot in resonance " << this->getResonanceName() << std::endl; std::cerr << " : Index must be in range 0 to " << nknots-1 << std::endl; return; } const std::vector& masses = this->getMasses(); std::vector& realVals = this->getAmp1Vals(); std::vector& imagVals = this->getAmp2Vals(); std::vector& realPars = this->getAmp1Pars(); std::vector& imagPars = this->getAmp2Pars(); realVals[knot] = realVal; realPars[knot]->value(realVal); realPars[knot]->genValue(realVal); realPars[knot]->initValue(realVal); realPars[knot]->fixed(fixRealPart); imagVals[knot] = imagVal; imagPars[knot]->value(imagVal); imagPars[knot]->genValue(imagVal); imagPars[knot]->initValue(imagVal); imagPars[knot]->fixed(fixImagPart); if ( knot == 0 ) { std::cout << "INFO in LauModIndPartWaveRealImag::setKnotAmp : Knot updated in resonance " << this->getResonanceName() << " at lower kinematic limit (" << masses[knot] << ")" << std::endl; } else if ( knot == nknots-1 ) { std::cout << "INFO in LauModIndPartWaveRealImag::setKnotAmp : Knot updated in resonance " << this->getResonanceName() << " at upper kinematic limit (" << masses[knot] << ")" << std::endl; } else { std::cout << "INFO in LauModIndPartWaveRealImag::setKnotAmp : Knot updated in resonance " << this->getResonanceName() << " at mass " << masses[knot] << std::endl; } if ( fixRealPart ) { std::cout << " : Real part fixed to " << realVal << std::endl; } else { std::cout << " : Real part set to " << realVal << std::endl; } if ( fixImagPart ) { std::cout << " : Imaginary part fixed to " << imagVal << std::endl; } else { std::cout << " : Imaginary part set to " << imagVal << std::endl; } } void LauModIndPartWaveRealImag::evaluateAmplitude(const Double_t mass) { - const Lau1DCubicSpline* splineReal = this->getSpline1(); - const Lau1DCubicSpline* splineImag = this->getSpline2(); + const Lau1DCubicSpline& splineReal { this->getSpline1() }; + const Lau1DCubicSpline& splineImag { this->getSpline2() }; - const Double_t re = splineReal->evaluate(mass); - const Double_t im = splineImag->evaluate(mass); + const Double_t re { splineReal.evaluate(mass) }; + const Double_t im { splineImag.evaluate(mass) }; this->setAmp(re, im); }