diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc index 1673c41..fae3736 100644 --- a/examples/Test_Dpipi.cc +++ b/examples/Test_Dpipi.cc @@ -1,440 +1,447 @@ #include using std::cout; using std::cerr; using std::endl; #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 "LauDaughters.hh" #include "LauDecayTimePdf.hh" #include "LauEffModel.hh" #include "LauIsobarDynamics.hh" #include "LauMagPhaseCoeffSet.hh" #include "LauRandom.hh" #include "LauRealImagCoeffSet.hh" #include "LauTimeDepFitModel.hh" #include "LauVetoes.hh" #include "LauFlavTag.hh" #include "Lau1DHistPdf.hh" #include "Lau1DCubicSpline.hh" void usage(const TString& progName) { cerr<<"Usage:"< vector_argv{argv, argc+argv}; LauDecayTimePdf::EfficiencyMethod effiMethod = LauDecayTimePdf::EfficiencyMethod::Spline; auto histTest = std::find(std::begin(vector_argv), std::end(vector_argv), "--hist"); if(histTest != std::end(vector_argv)) { effiMethod = LauDecayTimePdf::EfficiencyMethod::Binned; vector_argv.erase(histTest); } histTest = std::find(std::begin(vector_argv), std::end(vector_argv), "--flat"); if(histTest != std::end(vector_argv)) { effiMethod = LauDecayTimePdf::EfficiencyMethod::Flat; vector_argv.erase(histTest); } const TString dtype = vector_argv[2]; Int_t port = 0; Int_t iFit = 0; Int_t firstExpt = 0; Int_t firstExptGen = 0; Int_t nExpt = 1; Int_t nExptGen = 1; LauTimeDepFitModel::CPEigenvalue eigenvalue = LauTimeDepFitModel::CPEven; if (dtype=="CPEven") { eigenvalue = LauTimeDepFitModel::CPEven; } else if (dtype=="CPOdd") { eigenvalue = LauTimeDepFitModel::CPOdd; } else if (dtype=="QFS") { eigenvalue = LauTimeDepFitModel::QFS; } else { usage(vector_argv[0]); return EXIT_FAILURE; } Bool_t fixPhiMix(kFALSE); Bool_t useSinCos(kTRUE); // check the command line arguments if (argc<1) { usage(vector_argv[0]); return EXIT_FAILURE; } TString command = vector_argv[1]; if (command != "gen" && command != "fit") { usage(vector_argv[0]); return EXIT_FAILURE; } if (command == "fit") { if (argc>3) { port = atoi(vector_argv[3]); if (argc>4) { iFit = atoi(vector_argv[4]); if (argc>5) { firstExpt = atoi(vector_argv[5]); if (argc>6) { nExpt = atoi(vector_argv[6]); if (argc>7) { nExptGen = atoi(vector_argv[7]); } } } } } for (firstExptGen = 0; firstExptGen<(firstExpt+nExpt); firstExptGen+=nExptGen) { } firstExptGen -= nExptGen; if ( (nExpt > nExptGen) || (nExptGen%nExpt != 0) ) { cerr<<"Error, nExpt must be a factor of nExptGen."<3) { firstExptGen = atoi(vector_argv[3]); if (argc>4) { nExptGen = atoi(vector_argv[4]); } } } cout<<"dtype "<addMassVeto( 2, 2.00776, 2.01276 ); daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0"); daughtersB0 = new LauDaughters("B0", "pi+", "pi-", "D0_bar"); // efficiency effModelB0bar = new LauEffModel(daughtersB0bar, vetoes); effModelB0 = new LauEffModel(daughtersB0, vetoes); LauIsobarDynamics* sigModelB0bar(0); LauIsobarDynamics* sigModelB0(0); LauTimeDepFitModel* fitModel(0); std::vector params; //Args for flavTag: useAveDelta - kFALSE and useEtaPrime - kFALSE LauFlavTag* flavTag = new LauFlavTag(kFALSE,kFALSE); flavTag->setTrueTagVarName("trueTag"); TFile* etafile(0); TH1* etahist(0); Lau1DHistPdf* etahistpdf(0); etafile = TFile::Open("histogram.root"); etahist = dynamic_cast(etafile->Get("htemp")); etahistpdf = new Lau1DHistPdf("eta",etahist,0.0,0.54,kTRUE,kFALSE); // signal dynamics 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); 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 fitModel = new LauTimeDepFitModel(sigModelB0bar,sigModelB0,flavTag); fitModel->setASqMaxValue(1.0); 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); } fitModel->setCPEigenvalue( eigenvalue ); fitModel->setPhiMix( 2.0*LauConstants::beta, fixPhiMix, useSinCos ); fitModel->setAsymmetries(0.0,kTRUE); std::pair tagEff (0.5, 0.5); std::pair calib0 (0.4, 0.4); std::pair calib1 (1.0, 1.0); flavTag->addTagger("OSTagger","tagVal_OS","mistagVal_OS", etahistpdf,tagEff,calib0,calib1); // Delta t PDFs const Double_t minDt(0.0); const Double_t maxDt(20.0); const Double_t minDtErr(0.0); const Double_t maxDtErr(2.5); const Int_t nGauss(3); std::vector scale(nGauss); scale[0] = kTRUE; scale[1] = kTRUE; scale[2] = kFALSE; std::vector dtPars(10); TString mean0Name("dt_"); mean0Name += "_mean_0"; TString mean1Name("dt_"); mean1Name += "_mean_1"; TString mean2Name("dt_"); mean2Name += "_mean_2"; TString sigma0Name("dt_"); sigma0Name += "_sigma_0"; TString sigma1Name("dt_"); sigma1Name += "_sigma_1"; TString sigma2Name("dt_"); sigma2Name += "_sigma_2"; TString frac1Name("dt_"); frac1Name += "_frac_1"; TString frac2Name("dt_"); frac2Name += "_frac_2"; TString tauName("dt_"); tauName += "_tau"; TString freqName("dt_"); freqName += "_deltaM"; LauParameter * mean1 = new LauParameter(mean1Name, -1.27); LauParameter * mean2 = new LauParameter(mean2Name, 0.0); LauParameter * sigma1 = new LauParameter(sigma1Name, 3.0); LauParameter * sigma2 = new LauParameter(sigma2Name, 8.0); LauParameter * frac1 = new LauParameter(frac1Name, 0.0930); LauParameter * frac2 = new LauParameter(frac2Name, 0.0036); LauParameter * tau = new LauParameter(tauName, 1.520); LauParameter * freq = new LauParameter(freqName, 0.5064); TString mean0tagcat63Name("dt_"); mean0tagcat63Name += 63; mean0tagcat63Name += "_mean_0"; TString sigma0tagcat63Name("dt_"); sigma0tagcat63Name += 63; sigma0tagcat63Name += "_sigma_0"; LauParameter * mean0tagcat63 = new LauParameter(mean0tagcat63Name, -0.031); LauParameter * sigma0tagcat63 = new LauParameter(sigma0tagcat63Name, 0.972); //Decay time acceptance spline - same for all tag cats (though doesn't have to be) std::vector dtvals; dtvals.push_back(0.0); dtvals.push_back(1.0); dtvals.push_back(2.0); dtvals.push_back(4.0); dtvals.push_back(6.0); dtvals.push_back(8.0); dtvals.push_back(11.0); dtvals.push_back(14.0); dtvals.push_back(17.0); dtvals.push_back(20.0); std::vector effvals; effvals.push_back(0.0); effvals.push_back(0.15); effvals.push_back(0.25); effvals.push_back(0.33); effvals.push_back(0.38); effvals.push_back(0.4); effvals.push_back(0.43); effvals.push_back(0.45); effvals.push_back(0.47); effvals.push_back(0.50); //Lau1DCubicSpline* dtEffSpline = new Lau1DCubicSpline(dtvals,effvals); Lau1DCubicSpline* dtEffSpline = nullptr; TH1D* dtEffHist = nullptr; switch(effiMethod) { case LauDecayTimePdf::EfficiencyMethod::Spline: { - dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::StandardSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); + dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); + //TCanvas canvas("canvas","canvas",600,400); + //TF1* function = dtEffSpline->makeTF1(); + //function -> Draw(); + //canvas.SaveAs("testTF1.pdf"); + //return 0; break; } case LauDecayTimePdf::EfficiencyMethod::Binned: { const UInt_t upscale = 5; //This value must be >= 1 this linearly interpolates a histogram between the points of the spline const Int_t nBins = upscale * (dtvals.size()-1) ; Double_t edges[nBins + 1]; for(Int_t i = 0; i < nBins + 1; ++i) { Double_t binWidth = dtvals[i/upscale + 1] - dtvals[i/upscale]; //This is the width of the bin if upscale is 1, it is the distance between the current 2 knots edges[i] = dtvals[i/upscale] + (i%upscale)*binWidth/upscale; } cout << "Bins: "; for (int i = 0; i < nBins + 1; ++i){cout << edges[i] << ", ";} cout << endl; Double_t binFilling[nBins]; dtEffHist = new TH1D("dtEffHist","Histogram of efficiencies", nBins, edges); for(Int_t i = 0; i < nBins; ++i) { Double_t lknot = effvals[i/upscale]; Double_t rknot = effvals[i/upscale + 1]; Double_t pos = 0.5*(edges[i]+edges[i+1]); //edges[i] + (edges[i+1]-edges[i]) * (1.0*(i%upscale)/upscale); Double_t weight = (pos - dtvals[i/upscale])/(dtvals[i/upscale + 1] - dtvals[i/upscale]); Double_t filling= lknot*(1-weight) + rknot*weight; cout << "lknot: " << lknot << "; rknot: " << rknot << "; pos: " << pos << "; filling: " << filling << endl; binFilling[i] = filling; dtEffHist->SetBinContent(i+1, binFilling[i]); } break; } case LauDecayTimePdf::EfficiencyMethod::Flat: {break;} } dtPars[0] = mean0tagcat63; dtPars[1] = mean1->createClone(); dtPars[2] = mean2->createClone(); dtPars[3] = sigma0tagcat63; dtPars[4] = sigma1->createClone(); dtPars[5] = sigma2->createClone(); dtPars[6] = frac1->createClone(); dtPars[7] = frac2->createClone(); dtPars[8] = tau->createClone(); dtPars[9] = freq->createClone(); if (dtype=="CPEven"){ LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::ExpTrig, nGauss, scale, LauDecayTimePdf::DecayTime, effiMethod ); dtPdf->doSmearing(kFALSE); //dtPdf->doSmearing(kTRUE); switch(effiMethod) { case LauDecayTimePdf::EfficiencyMethod::Spline: dtPdf->setEffiSpline(dtEffSpline); break; case LauDecayTimePdf::EfficiencyMethod::Binned: dtPdf->setEffiHist(dtEffHist); break; case LauDecayTimePdf::EfficiencyMethod::Flat: break; } fitModel->setSignalDtPdf( dtPdf ); } else { LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::ExpTrig, nGauss, scale, LauDecayTimePdf::DecayTime, effiMethod ); dtPdf->doSmearing(kFALSE); //dtPdf->doSmearing(kTRUE); switch(effiMethod) { case LauDecayTimePdf::EfficiencyMethod::Spline: dtPdf->setEffiSpline(dtEffSpline); break; case LauDecayTimePdf::EfficiencyMethod::Binned: dtPdf->setEffiHist(dtEffHist); break; case LauDecayTimePdf::EfficiencyMethod::Flat: break; } fitModel->setSignalDtPdf( dtPdf ); } // set the number of signal events cout<<"nSigEvents = "<setNSigEvents(nSigPar); // set the number of experiments if (command == "fit") { fitModel->setNExpts(nExpt, firstExpt); } else { fitModel->setNExpts(nExptGen, firstExptGen); } fitModel->useAsymmFitErrors(kFALSE); //fitModel->useRandomInitFitPars(kTRUE); fitModel->useRandomInitFitPars(kFALSE); fitModel->doPoissonSmearing(kFALSE); fitModel->doEMLFit(kFALSE); fitModel->writeLatexTable(kFALSE); TString dataFile(""); TString treeName("fitTree"); TString rootFileName(""); TString tableFileName(""); TString fitToyFileName(""); TString splotFileName(""); dataFile = "TEST-Dpipi"; switch(effiMethod) { case LauDecayTimePdf::EfficiencyMethod::Spline: dataFile += "_Spline"; break; case LauDecayTimePdf::EfficiencyMethod::Binned: dataFile += "_Hist"; break; case LauDecayTimePdf::EfficiencyMethod::Flat: dataFile += "_Flat"; break; } dataFile += "_"+dtype+"_expts"; dataFile += firstExptGen; dataFile += "-"; dataFile += firstExptGen+nExptGen-1; if ( dtype != "QFS" ) { dataFile += "_CP"; if ( eigenvalue == LauTimeDepFitModel::CPEven ) { dataFile += "even"; } else { dataFile += "odd"; } } dataFile += ".root"; if (command == "fit") { rootFileName = "fit"; rootFileName += iFit; rootFileName += "_Results_"; rootFileName += dtype; rootFileName += "_expts"; rootFileName += firstExpt; rootFileName += "-"; rootFileName += firstExpt+nExpt-1; rootFileName += ".root"; tableFileName = "fit"; tableFileName += iFit; tableFileName += "_Results_"; tableFileName += dtype; tableFileName += "_expts"; tableFileName += firstExpt; tableFileName += "-"; tableFileName += firstExpt+nExpt-1; fitToyFileName = "fit"; fitToyFileName += iFit; fitToyFileName += "_ToyMC_"; fitToyFileName += dtype; fitToyFileName += "_expts"; fitToyFileName += firstExpt; fitToyFileName += "-"; fitToyFileName += firstExpt+nExpt-1; fitToyFileName += ".root"; splotFileName = "fit"; splotFileName += iFit; splotFileName += "_sPlot_"; splotFileName += dtype; splotFileName += "_expts"; splotFileName += firstExpt; splotFileName += "-"; splotFileName += firstExpt+nExpt-1; splotFileName += ".root"; } else { rootFileName = "dummy.root"; tableFileName = "genResults"; } // 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 if ( command == "fit" ){ fitModel->runTask( dataFile, treeName, rootFileName, tableFileName, "localhost", port ); } else { fitModel->run( command, dataFile, treeName, rootFileName, tableFileName ); } return EXIT_SUCCESS; } diff --git a/inc/Lau1DCubicSpline.hh b/inc/Lau1DCubicSpline.hh index b1ee48a..ba922b9 100644 --- a/inc/Lau1DCubicSpline.hh +++ b/inc/Lau1DCubicSpline.hh @@ -1,213 +1,224 @@ /* 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 "LauParameter.hh" #include "Rtypes.h" +#include "TF1.h" class Lau1DCubicSpline { public: //! Define the allowed interpolation types enum LauSplineType { 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 { 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] 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); //! Destructor virtual ~Lau1DCubicSpline(); //! 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; //! Update the y-values of the knots /*! \param [in] ys the y-values of the knots */ void updateYValues(const std::vector& ys); void updateYValues(const std::vector& ys); //! Update the type of interpolation to perform /*! \param [in] type the type of interpolation */ void updateType(LauSplineType 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); //! Return the number of knots UInt_t getnKnots() const {return nKnots_;} //! Get y values std::vector getYValues(){return y_;} + //! Get the coefficients of spline section i in the form a + bx + cx^2 + dx^3 + /*! + /params [in] i refers to the left-hand index of the knot (i = 0 gets the params between x_0 and x_1) + */ + std::array getCoefficients(UInt_t i); + + //! Make a TF1 showing the spline with its current knot values. Mostly as a sanity check + TF1* makeTF1(); + private: //! Copy constructor - not implemented Lau1DCubicSpline( const Lau1DCubicSpline& rhs ); //! Copy assignment operator - not implemented Lau1DCubicSpline& operator=(const Lau1DCubicSpline& rhs); //! 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 UInt_t nKnots_; //! The x-value at each knot 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_; //! The left-hand boundary condition on the spline LauSplineBoundaryType leftBound_; //! The right-hand boundary condition on the spline LauSplineBoundaryType 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/src/Lau1DCubicSpline.cc b/src/Lau1DCubicSpline.cc index be9e5f0..e0c05a6 100644 --- a/src/Lau1DCubicSpline.cc +++ b/src/Lau1DCubicSpline.cc @@ -1,351 +1,455 @@ /* 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" 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) { init(); } Lau1DCubicSpline::~Lau1DCubicSpline() { } Double_t Lau1DCubicSpline::evaluate(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); 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]; if(type_ == Lau1DCubicSpline::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); 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 (UInt_t i=0; iunblindValue(); } this->calcDerivatives(); } void Lau1DCubicSpline::updateType(LauSplineType type) { if(type_ != type) { type_ = type; this->calcDerivatives(); } } void Lau1DCubicSpline::updateBoundaryConditions(LauSplineBoundaryType leftBound, LauSplineBoundaryType rightBound, Double_t dydx0, Double_t dydxn) { Bool_t updateDerivatives(kFALSE); if(leftBound_ != leftBound || rightBound_ != rightBound) { leftBound_ = leftBound; rightBound_ = rightBound; updateDerivatives = kTRUE; } if(dydx0_ != dydx0) { dydx0_ = dydx0; if(leftBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = kTRUE; } if(dydxn_ != dydxn) { dydxn_ = dydxn; if(rightBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = kTRUE; } if(updateDerivatives) { this->calcDerivatives(); } } +std::array Lau1DCubicSpline::getCoefficients(UInt_t i) +{ + std::array result = {0.,0.,0.,0.}; + if(i >= nKnots_) + { + std::cerr << "ERROR in Lau1DCubicSpline::getCoefficients requested for too high a knot value" << std::endl; + return result; + } + + Double_t xL = x_[i] , xH = x_[i+1]; + Double_t yL = y_[i] , yH = y_[i+1]; + Double_t h = xH-xL; //This number comes up a lot + + switch(type_) + { + case Lau1DCubicSpline::StandardSpline: + case Lau1DCubicSpline::AkimaSpline: + { + Double_t kL = dydx_[i], kH = dydx_[i+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 + { + //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_[i+1]/x_[i]; //a denominator used for p[2] and p[3] + std::array p = {y_[i],t(i),0.,0.}; //we'll do the last 2 below + + p[2] = 3*m(i)-2*t(i)-t(i+1); + p[2]/= pDenom; + + p[3] = t(i)+t(i+1)-2*m(i); + 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: + { + result[0] = xH*yL-xL*yH; + result[1] = yH-yL; + + for(auto& res : result){res /= h;} + break; + } + } + + return result; +} + +TF1* Lau1DCubicSpline::makeTF1() +{ + TString functionString = ""; //make a long piecewise construction of all the spline pieces + for(UInt_t i = 0; i < nKnots_-1; ++i) + { + functionString += Form("(x>%f && x<= %f)*",x_[i],x_[i+1]);//get the bounds of this piece + std::array coeffs = this->getCoefficients(i); + functionString += Form("(%f + %f*x + %f*x^2 + %f*x^3)",coeffs[0],coeffs[1],coeffs[2],coeffs[3]); + if(i < nKnots_ -2){functionString += " + \n";}//add to all lines except the last + } + + TF1* func = new TF1("SplineFunction", functionString, x_.front(), x_.back()); + + return func; +} + 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); } dydx_.insert(dydx_.begin(),nKnots_,0.); a_.insert(a_.begin(),nKnots_,0.); b_.insert(b_.begin(),nKnots_,0.); c_.insert(c_.begin(),nKnots_,0.); d_.insert(d_.begin(),nKnots_,0.); this->calcDerivatives(); } void Lau1DCubicSpline::calcDerivatives() { switch ( type_ ) { case Lau1DCubicSpline::StandardSpline : this->calcDerivativesStandard(); break; case Lau1DCubicSpline::AkimaSpline : this->calcDerivativesAkima(); break; case Lau1DCubicSpline::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_; } // 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_; } // the remaining equations ensure that f_i-1''(x_i) = f''_i(x_i) for all internal knots for(UInt_t i=1; 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.); // 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(UInt_t i=1; i