diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc index 2e9f967..e7ccc96 100644 --- a/examples/Test_Dpipi.cc +++ b/examples/Test_Dpipi.cc @@ -1,327 +1,332 @@ #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 "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" #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; } const Bool_t fixPhiMix{ settings.fixPhiMix || settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS }; const Bool_t useSinCos{kTRUE}; 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; } LauDaughters* daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0"); LauDaughters* daughtersB0 = new LauDaughters("B0", "pi+", "pi-", "D0_bar"); // efficiency LauVetoes* vetoes = new LauVetoes(); //vetoes->addMassVeto( 2, 2.00776, 2.01276 ); LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes); LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes); //Args for flavTag: useAveDelta - kFALSE and useEtaPrime - kFALSE LauFlavTag* flavTag = new LauFlavTag(kFALSE,kFALSE); flavTag->setTrueTagVarName("trueTag"); if (settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS) { flavTag->setDecayFlvVarName("decayFlv"); } TFile* etaFile = TFile::Open("ft-eta-hist.root"); TH1* etaHist = dynamic_cast(etaFile->Get("ft_eta_hist")); Lau1DHistPdf* etaHistPDF = new Lau1DHistPdf("eta",etaHist,0.0,0.5,kTRUE,kFALSE); const Double_t meanEta { etaHistPDF->getMean() }; // if the tagging is perfect then also make it perfectly efficient, otherwise 50% efficient const Double_t tagEffVal { (meanEta == 0.0) ? 1.0 : 0.5 }; std::pair tagEff {tagEffVal, tagEffVal}; // use a null calibration for the time being, so p0 = and p1 = 1 std::pair calib0 {meanEta, meanEta}; std::pair calib1 {1.0, 1.0}; flavTag->addTagger("OSTagger", "tagVal_OS", "mistagVal_OS", etaHistPDF, tagEff, calib0, calib1); // 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); } fitModel->setCPEigenvalue( settings.dType ); fitModel->setPhiMix( 2.0*LauConstants::beta, fixPhiMix, useSinCos ); // production asymmetry fitModel->setAsymmetries(0.0,kTRUE); // Delta t PDFs const Double_t minDt(0.0); const Double_t maxDt(15.0); const Double_t minDtErr(0.0); const Double_t maxDtErr(0.215); const std::vector scale { settings.perEventTimeErr && kTRUE, settings.perEventTimeErr && kTRUE, }; const UInt_t nGauss(scale.size()); LauParameter * mean0 = new LauParameter("dt_mean_0", scale[0] ? -1.63e-3 : -1.84e-03, -0.01, 0.01, kTRUE ); LauParameter * mean1 = new LauParameter("dt_mean_1", scale[1] ? -1.63e-3 : -3.62e-03, -0.01, 0.01, kTRUE ); LauParameter * sigma0 = new LauParameter("dt_sigma_0", scale[0] ? 0.991 : 3.05e-02, 0.0, 2.0, kTRUE ); LauParameter * sigma1 = new LauParameter("dt_sigma_1", scale[1] ? 1.80 : 6.22e-02, 0.0, 2.5, kTRUE ); LauParameter * frac1 = new LauParameter("dt_frac_1", scale[0] && scale[1] ? 0.065 : 0.761, 0.0, 1.0, kTRUE); LauParameter * tau = new LauParameter("dt_tau", 1.520, 0.5, 5.0, settings.fixLifetime); LauParameter * freq = new LauParameter("dt_deltaM", 0.5064, 0.0, 1.0, settings.fixDeltaM); std::vector dtPars { mean0, mean1, sigma0, sigma1, frac1, tau, freq }; // 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 y-values from the histogram contents std::vector dtvals; std::vector effvals; dtvals.push_back(minDt); effvals.push_back(0.0); for ( Int_t bin{0}; bin < dtaHist->GetNbinsX(); ++bin ) { dtvals.push_back( dtaHist->GetBinCenter(bin+1) ); effvals.push_back( dtaHist->GetBinContent(bin+1) ); } dtvals.push_back(maxDt); effvals.push_back(effvals.back()); + //DEBUG + effvals.clear(); + effvals.resize(dtvals.size()); + for(size_t i = 0; i < dtvals.size(); ++i){effvals[i] = 1.;} + // Decay time error histogram TFile* dteFile = TFile::Open("dte-hist.root"); TH1* dteHist = dynamic_cast(dteFile->Get("dte_hist")); LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "decayTime", "decayTimeErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::ExpTrig, nGauss, scale, LauDecayTimePdf::DecayTime, settings.timeEffModel ); dtPdf->doSmearing(settings.timeResolution); if ( settings.perEventTimeErr ) { dtPdf->setErrorHisto( dteHist ); } switch(settings.timeEffModel) { case LauDecayTimePdf::EfficiencyMethod::Spline: { fitModel->setASqMaxValue(0.06); - Lau1DCubicSpline* dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); + Lau1DCubicSpline* dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::StandardSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); dtPdf->setEffiSpline(dtEffSpline); break; } case LauDecayTimePdf::EfficiencyMethod::Binned: { fitModel->setASqMaxValue(0.06); dtPdf->setEffiHist(dtaHist); break; } case LauDecayTimePdf::EfficiencyMethod::Flat: { fitModel->setASqMaxValue(4.1); break; } } fitModel->setSignalDtPdf( dtPdf ); // set the number of signal events 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->doPoissonSmearing(kFALSE); fitModel->doEMLFit(kFALSE); fitModel->writeLatexTable(kFALSE); 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 LauDecayTimePdf::EfficiencyMethod::Spline: dataFile += "_Spline"; break; case LauDecayTimePdf::EfficiencyMethod::Binned: dataFile += "_Hist"; break; case LauDecayTimePdf::EfficiencyMethod::Flat: dataFile += "_Flat"; 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/src/Lau1DCubicSpline.cc b/src/Lau1DCubicSpline.cc index 3a1a4ac..845954a 100644 --- a/src/Lau1DCubicSpline.cc +++ b/src/Lau1DCubicSpline.cc @@ -1,485 +1,486 @@ /* 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(const UInt_t i, const bool normalise) const { std::array result = {0.,0.,0.,0.}; + if(i >= nKnots_-1) { 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; } } if(normalise) { Double_t integral = this->integral(); for(auto& res : result){res /= integral;} } return result; } Double_t Lau1DCubicSpline::integral() const { Double_t integral = 0.; for(UInt_t iKnot = 0; iKnot < nKnots_ -1; ++iKnot) { Double_t minAbs = x_[iKnot]; Double_t maxAbs = x_[iKnot+1]; std::array coeffs = this -> getCoefficients(iKnot, false); auto integralFunc = [&coeffs](Double_t x){return coeffs[0]*x + coeffs[1]*x*x/2 + coeffs[2]*x*x*x/3 + coeffs[3]*x*x*x*x/4;}; 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(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,normalise); 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); } 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 : 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 #include #include #include #include #include "TMath.h" #include "TRandom.h" #include "TSystem.h" #include "TH1.h" #include "RooMath.h" #include "Lau1DCubicSpline.hh" #include "Lau1DHistPdf.hh" #include "LauConstants.hh" #include "LauComplex.hh" #include "LauDecayTimePdf.hh" #include "LauFitDataTree.hh" #include "LauParameter.hh" #include "LauRandom.hh" ClassImp(LauDecayTimePdf) LauDecayTimePdf::LauDecayTimePdf(const TString& theVarName, const TString& theVarErrName, const std::vector& params, Double_t minAbscissaVal, Double_t maxAbscissaVal, Double_t minAbscissaErr, Double_t maxAbscissaErr, FuncType type, UInt_t nGauss, const std::vector& scale, const TimeMeasurementMethod method, const EfficiencyMethod effMethod) : varName_(theVarName), varErrName_(theVarErrName), param_(params), smear_(kTRUE), minAbscissa_(minAbscissaVal), maxAbscissa_(maxAbscissaVal), minAbscissaError_(minAbscissaErr), maxAbscissaError_(maxAbscissaErr), abscissaError_(0.0), abscissaErrorGenerated_(kFALSE), errorDistMPV_(0.230), // for signal 0.234, for qqbar 0.286 errorDistSigma_(0.075), // for signal 0.073, for qqbar 0.102 nGauss_(nGauss), mean_(nGauss_,0), sigma_(nGauss_,0), frac_(nGauss_-1,0), tau_(0), deltaM_(0), deltaGamma_(0), fracPrompt_(0), type_(type), method_(method), effMethod_(effMethod), scaleMeans_(scale), scaleWidths_(scale), scaleWithPerEventError_( std::accumulate( scale.begin(), scale.end(), kFALSE, std::logical_or() ) ), expTerm_(0.0), cosTerm_(0.0), sinTerm_(0.0), coshTerm_(0.0), sinhTerm_(0.0), normTermExp_(0.0), normTermCosh_(0.0), normTermSinh_(0.0), errTerm_(0.0), effiTerm_(0.0), pdfTerm_(0.0), errHist_(nullptr), pdfHist_(nullptr), effiFun_(nullptr), effiHist_(nullptr), effiPars_(0) { this->initialise(); } LauDecayTimePdf::LauDecayTimePdf(const TString& theVarName, const TString& theVarErrName, const std::vector& params, Double_t minAbscissaVal, Double_t maxAbscissaVal, Double_t minAbscissaErr, Double_t maxAbscissaErr, FuncType type, UInt_t nGauss, const std::vector& scaleMeans, const std::vector& scaleWidths, const TimeMeasurementMethod method, const EfficiencyMethod effMethod) : varName_(theVarName), varErrName_(theVarErrName), param_(params), smear_(kTRUE), minAbscissa_(minAbscissaVal), maxAbscissa_(maxAbscissaVal), minAbscissaError_(minAbscissaErr), maxAbscissaError_(maxAbscissaErr), abscissaError_(0.0), abscissaErrorGenerated_(kFALSE), errorDistMPV_(0.230), // for signal 0.234, for qqbar 0.286 errorDistSigma_(0.075), // for signal 0.073, for qqbar 0.102 nGauss_(nGauss), mean_(nGauss_,0), sigma_(nGauss_,0), frac_(nGauss_-1,0), tau_(0), deltaM_(0), deltaGamma_(0), fracPrompt_(0), type_(type), method_(method), effMethod_(effMethod), scaleMeans_(scaleMeans), scaleWidths_(scaleWidths), scaleWithPerEventError_( std::accumulate( scaleMeans.begin(), scaleMeans.end(), kFALSE, std::logical_or() ) || std::accumulate( scaleWidths.begin(), scaleWidths.end(), kFALSE, std::logical_or() ) ), expTerm_(0.0), cosTerm_(0.0), sinTerm_(0.0), coshTerm_(0.0), sinhTerm_(0.0), normTermExp_(0.0), normTermCosh_(0.0), normTermSinh_(0.0), errTerm_(0.0), effiTerm_(0.0), pdfTerm_(0.0), errHist_(nullptr), pdfHist_(nullptr), effiFun_(nullptr), effiHist_(nullptr), effiPars_(0) { this->initialise(); } LauDecayTimePdf::~LauDecayTimePdf() { // Destructor delete errHist_; errHist_ = nullptr; delete pdfHist_; pdfHist_ = nullptr; delete effiFun_; effiFun_ = nullptr; delete effiHist_; effiHist_ = nullptr; for( auto& par : effiPars_ ){ delete par; par = nullptr; } effiPars_.clear(); } void LauDecayTimePdf::initialise() { // The parameters are: // - the mean and the sigma (bias and spread in resolution) of the gaussian(s) // - the mean lifetime, denoted tau, of the exponential decay // - the frequency of oscillation, denoted Delta m, of the cosine and sine terms // - the decay width difference, denoted Delta Gamma, of the hyperbolic cosine and sine terms // // The next two arguments specify the range in which the PDF is defined, // and the PDF will be normalised w.r.t. these limits. // // The final three arguments define the type of Delta t PDF (Delta, Exp, ExpTrig or ExpHypTrig ), the number of gaussians // and whether or not the gaussian parameters should be scaled by the per-event errors on Delta t // First check whether the scale vector is nGauss in size if (nGauss_ != scaleMeans_.size() || nGauss_ != scaleWidths_.size()) { std::cerr<<"ERROR in LauDecayTimePdf::initialise : scale vector size not the same as nGauss."<Exit(EXIT_FAILURE); } if (type_ == Hist) { if (this->nParameters() != 0){ std::cerr<<"ERROR in LauDecayTimePdf::initialise : Hist PDF should have 0 parameters"<Exit(EXIT_FAILURE); } } else { TString meanName("mean_"); TString sigmaName("sigma_"); TString fracName("frac_"); Bool_t foundParams(kTRUE); for (UInt_t i(0); ifindParameter(tempName); foundParams &= (mean_[i] != 0); sigma_[i] = this->findParameter(tempName2); foundParams &= (sigma_[i] != 0); if (i!=0) { frac_[i-1] = this->findParameter(tempName3); foundParams &= (frac_[i-1] != 0); } } if (type_ == Delta) { if ((this->nParameters() != (3*nGauss_-1)) || (!foundParams)) { std::cerr<<"ERROR in LauDecayTimePdf::initialise : Delta type PDF requires:"<Exit(EXIT_FAILURE); } } else if (type_ == Exp) { tau_ = this->findParameter("tau"); foundParams &= (tau_ != 0); if ((this->nParameters() != (3*nGauss_-1+1)) || (!foundParams)) { std::cerr<<"ERROR in LauDecayTimePdf::initialise : Exp type PDF requires:"<Exit(EXIT_FAILURE); } } else if (type_ == DeltaExp) { tau_ = this->findParameter("tau"); fracPrompt_ = this->findParameter("frac_prompt"); foundParams &= (tau_ != 0); foundParams &= (fracPrompt_ != 0); if ((this->nParameters() != (3*nGauss_-1+2)) || (!foundParams)) { std::cerr<<"ERROR in LauDecayTimePdf::initialise : DeltaExp type PDF requires:"<Exit(EXIT_FAILURE); } } else if (type_ == ExpTrig) { tau_ = this->findParameter("tau"); deltaM_ = this->findParameter("deltaM"); foundParams &= (tau_ != 0); foundParams &= (deltaM_ != 0); if ((this->nParameters() != (3*nGauss_-1+2)) || (!foundParams)) { std::cerr<<"ERROR in LauDecayTimePdf::initialise : ExpTrig type PDF requires:"<Exit(EXIT_FAILURE); } } else if (type_ == ExpHypTrig) { tau_ = this->findParameter("tau"); deltaM_ = this->findParameter("deltaM"); deltaGamma_ = this->findParameter("deltaGamma"); foundParams &= (tau_ != 0); foundParams &= (deltaM_ != 0); foundParams &= (deltaGamma_ != 0); if ((this->nParameters() != (3*nGauss_-1+3)) || (!foundParams)) { std::cerr<<"ERROR in LauDecayTimePdf::initialise : ExpHypTrig type PDF requires:"<Exit(EXIT_FAILURE); } } } } Double_t LauDecayTimePdf::effectiveResolution() const { Double_t dilution = 0.; Double_t dMSq = deltaM_->unblindValue() * deltaM_->unblindValue(); // Might be cleaner to just append this to the vector in the init step, // the the consistency can also be checked Double_t fracSum = 0; for (auto f : frac_) fracSum += f->unblindValue(); Double_t lastFrac = 1. - fracSum; for (size_t i = 0; i < sigma_.size(); i++) { Double_t sigSq = sigma_[i]->unblindValue() * sigma_[i]->unblindValue(); Double_t thisFrac = lastFrac; if (i < sigma_.size() - 1) thisFrac = frac_[i]->unblindValue(); dilution += thisFrac * TMath::Exp(-dMSq * sigSq / 2.); } return TMath::Sqrt(-2. * TMath::Log(dilution)) / deltaM_->unblindValue(); } void LauDecayTimePdf::cacheInfo(const LauFitDataTree& inputData) { // Check that the input data contains the decay time variable Bool_t hasBranch = inputData.haveBranch(this->varName()); if (!hasBranch) { std::cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Input data does not contain variable \""<varName()<<"\"."<varErrName()); if (!hasBranch) { std::cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Input data does not contain variable \""<varErrName()<<"\"."<cacheInfo(inputData); } if (type_ == Hist) { // Pass the data to the decay-time PDF for caching if ( pdfHist_ ) { pdfHist_->cacheInfo(inputData); } } else { // determine whether we are caching our PDF value //TODO //Bool_t doCaching( this->nFixedParameters() == this->nParameters() ); //this->cachePDF( doCaching ); // clear the vectors and reserve enough space const UInt_t nEvents = inputData.nEvents(); abscissas_.clear(); abscissas_.reserve(nEvents); abscissaErrors_.clear(); abscissaErrors_.reserve(nEvents); expTerms_.clear(); expTerms_.reserve(nEvents); cosTerms_.clear(); cosTerms_.reserve(nEvents); sinTerms_.clear(); sinTerms_.reserve(nEvents); coshTerms_.clear(); coshTerms_.reserve(nEvents); sinhTerms_.clear(); sinhTerms_.reserve(nEvents); normTermsExp_.clear(); normTermsExp_.reserve(nEvents); normTermsCos_.clear(); normTermsCos_.reserve(nEvents); normTermsSin_.clear(); normTermsSin_.reserve(nEvents); normTermsCosh_.clear(); normTermsCosh_.reserve(nEvents); normTermsSinh_.clear(); normTermsSinh_.reserve(nEvents); effiTerms_.clear(); effiTerms_.reserve(nEvents); // If we're not using per-event information for the decay time // error, just calculate the normalisation terms once if ( ! scaleWithPerEventError_ ) { this->calcNorm(); } for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) { const LauFitData& dataValues = inputData.getData(iEvt); const Double_t abscissa { dataValues.at(this->varName()) }; if (abscissa > this->maxAbscissa() || abscissa < this->minAbscissa()) { std::cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Given value of the decay time: "<minAbscissa()<<","<maxAbscissa()<<"]."<Exit(EXIT_FAILURE); } abscissas_.push_back( abscissa ); const Double_t abscissaErr { scaleWithPerEventError_ ? dataValues.at(this->varErrName()) : 0.0 }; if ( scaleWithPerEventError_ && ( abscissaErr > this->maxAbscissaError() || abscissaErr < this->minAbscissaError() ) ) { std::cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Given value of the decay-time error: "<minAbscissaError()<<","<maxAbscissaError()<<"]."<Exit(EXIT_FAILURE); } abscissaErrors_.push_back(abscissaErr); // std::cout << "\033[1;31m IN CACHE INFO \033[0m" << std::endl; // std::cout << "\033[1;31m absErr: " << abscissaErr << "\033[0m" << std::endl; //DEBUG this->calcLikelihoodInfo(abscissa, abscissaErr); // If we are using per-event information for the decay // time error, need to calculate the normalisation // terms for every event if ( scaleWithPerEventError_ ) { this->calcNorm(abscissaErr); } expTerms_.push_back(expTerm_); cosTerms_.push_back(cosTerm_); sinTerms_.push_back(sinTerm_); coshTerms_.push_back(coshTerm_); sinhTerms_.push_back(sinhTerm_); normTermsExp_.push_back(normTermExp_); normTermsCos_.push_back(normTermCos_); normTermsSin_.push_back(normTermSin_); normTermsCosh_.push_back(normTermCosh_); normTermsSinh_.push_back(normTermSinh_); effiTerms_.push_back(effiTerm_); } } } void LauDecayTimePdf::calcLikelihoodInfo(const UInt_t iEvt) { // Extract all the terms and their normalisations if (type_ == Hist) { if ( pdfHist_ ) { pdfHist_->calcLikelihoodInfo(iEvt); pdfTerm_ = pdfHist_->getLikelihood(); } else { pdfTerm_ = 1.0; } } else { expTerm_ = expTerms_[iEvt]; cosTerm_ = cosTerms_[iEvt]; sinTerm_ = sinTerms_[iEvt]; coshTerm_ = coshTerms_[iEvt]; sinhTerm_ = sinhTerms_[iEvt]; normTermExp_ = normTermsExp_[iEvt]; normTermCos_ = normTermsCos_[iEvt]; normTermSin_ = normTermsSin_[iEvt]; normTermCosh_ = normTermsCosh_[iEvt]; normTermSinh_ = normTermsSinh_[iEvt]; } // Extract the decay time error PDF value if ( errHist_ ) { errHist_->calcLikelihoodInfo(iEvt); errTerm_ = errHist_->getLikelihood(); } else { errTerm_ = 1.0; } // Extract the decay time efficiency effiTerm_ = effiTerms_[iEvt]; // TODO - Parameters can change in some cases, so we'll need to update things! // - For the moment do the blunt force thing and recalculate everything for every event! // - Need to make this intelligent! const Double_t abscissa = abscissas_[iEvt]; const Double_t abscissaErr = abscissaErrors_[iEvt]; this->calcLikelihoodInfo(abscissa,abscissaErr); this->calcNorm(abscissaErr); } void LauDecayTimePdf::calcLikelihoodInfo(const Double_t abscissa) { // Check whether any of the gaussians should be scaled - if any of them should we need the per-event error if (scaleWithPerEventError_) { std::cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfo : Per-event error on decay time not provided, cannot calculate anything."<calcLikelihoodInfo(abscissa, 0.0); } void LauDecayTimePdf::calcLikelihoodInfo(const Double_t abscissa, const Double_t abscissaErr) { // Check that the decay time and the decay time error are in valid ranges if (abscissa > this->maxAbscissa() || abscissa < this->minAbscissa()) { std::cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfo : Given value of the decay time: "<minAbscissa()<<","<maxAbscissa()<<"]."<Exit(EXIT_FAILURE); } if ( scaleWithPerEventError_ && ( abscissaErr > this->maxAbscissaError() || abscissaErr < this->minAbscissaError() ) ) { std::cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfo : Given value of Delta t error: "<minAbscissaError()<<","<maxAbscissaError()<<"]."<Exit(EXIT_FAILURE); } // Determine the decay time efficiency switch( effMethod_ ) { case EfficiencyMethod::Spline : effiTerm_ = effiFun_ ? effiFun_ -> evaluate(abscissa) : 1.0 ; break; case EfficiencyMethod::Binned : effiTerm_ = effiHist_ ? effiHist_-> GetBinContent(effiHist_-> FindFixBin(abscissa)) : 1.0 ; break; case EfficiencyMethod::Flat : effiTerm_ = 1.0 ; break; } if ( effiTerm_ > 1.0 ) { effiTerm_ = 1.0; } else if ( effiTerm_ < 0.0 ) { effiTerm_ = 0.0; } // For the histogram PDF just calculate that term and return if (type_ == Hist){ if ( pdfHist_ ) { pdfHist_->calcLikelihoodInfo(abscissa); pdfTerm_ = pdfHist_->getLikelihood(); } else { pdfTerm_ = 1.0; } return; } // If we're not using the resolution function, calculate the simple terms and return if (!this->doSmearing()) { this->calcNonSmearedTerms(abscissa); return; } // Get all the up to date parameter values for the resolution function std::vector frac(nGauss_); std::vector mean(nGauss_); std::vector sigma(nGauss_); Double_t fracPrompt(0.0); // TODO - why do we do the fractions this way around? frac[0] = 1.0; for (UInt_t i(0); iunblindValue(); sigma[i] = sigma_[i]->unblindValue(); if (i != 0) { frac[i] = frac_[i-1]->unblindValue(); frac[0] -= frac[i]; } } if (type_ == DeltaExp) { fracPrompt = fracPrompt_->unblindValue(); } // Scale the gaussian parameters by the per-event error on Delta t (if appropriate) for (UInt_t i(0); i x(nGauss_); const Double_t xMax = this->maxAbscissa(); const Double_t xMin = this->minAbscissa(); for (UInt_t i(0); i 1e-10) { Double_t exponent(0.0); Double_t norm(0.0); Double_t scale = LauConstants::root2*sigma[i]; Double_t scale2 = LauConstants::rootPiBy2*sigma[i]; exponent = -0.5*x[i]*x[i]/(sigma[i]*sigma[i]); norm = scale2*(TMath::Erf((xMax - mean[i])/scale) - TMath::Erf((xMin - mean[i])/scale)); value += frac[i]*TMath::Exp(exponent)/norm; } } } if (type_ != Delta) { // Reset values of terms expTerm_ = 0.0; cosTerm_ = 0.0; sinTerm_ = 0.0; coshTerm_ = 0.0; sinhTerm_ = 0.0; // Calculate values of the PDF convoluted with each Gaussian for a given value of the abscsissa for (UInt_t i(0); i absErrVec = {abscissaErr}; //Otherwise seg fault errHist_->calcLikelihoodInfo(absErrVec); errTerm_ = errHist_->getLikelihood(); } else { errTerm_ = 1.0; } } void LauDecayTimePdf::calcNonSmearedTerms(Double_t abscissa) { // Reset values of terms errTerm_ = 1.0; expTerm_ = 0.0; cosTerm_ = 0.0; sinTerm_ = 0.0; coshTerm_ = 0.0; sinhTerm_ = 0.0; if ( type_ == Hist || type_ == Delta ){ return; } const Double_t tau { tau_->unblindValue() }; const Double_t gamma { 1.0 / tau }; if (method_ == DecayTime) { expTerm_ = TMath::Exp(-abscissa*gamma); } else if (method_ == DecayTimeDiff) { expTerm_ = TMath::Exp(-TMath::Abs(abscissa)*gamma); } // Calculate also the terms related to cosine and sine if (type_ == ExpTrig) { const Double_t deltaM = deltaM_->unblindValue(); coshTerm_ = expTerm_; sinhTerm_ = 0.0; cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_; } // Calculate also the terms related to cosh, sinh, cosine, and sine else if (type_ == ExpHypTrig) { const Double_t deltaM = deltaM_->unblindValue(); const Double_t deltaGamma = deltaGamma_->unblindValue(); coshTerm_ = TMath::CosH(0.5*deltaGamma*abscissa)*expTerm_; sinhTerm_ = TMath::SinH(0.5*deltaGamma*abscissa)*expTerm_; cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_; } } std::pair LauDecayTimePdf::smearedCosSinTerm(Double_t t, Double_t sigma, Double_t mu) { using namespace std::complex_literals; const Double_t gamma = 1. / this->tau_->unblindValue(); const Double_t x = (t - mu) / (LauConstants::root2 * sigma); const std::complex z = std::complex(gamma * sigma / LauConstants::root2, -this->deltaM_->unblindValue() * sigma / LauConstants::root2); const std::complex arg1 = std::complex(0., 1.) * (z - x); const std::complex arg2 { -(x*x) - (arg1 * arg1) }; const std::complex conv = arg1.imag() < -5.? 0.5 * std::exp(arg2) * RooMath::erfc( -1i * arg1 ) : 0.5 * TMath::Exp(-(x * x)) * RooMath::faddeeva(arg1) ; const Double_t cos_conv = conv.real(); const Double_t sin_conv = conv.imag(); return {cos_conv, sin_conv}; } std::pair LauDecayTimePdf::smearedCoshSinhTerm(Double_t t, Double_t sigma, Double_t mu) { using namespace std::complex_literals; Double_t gamma = 1. / this->tau_->unblindValue(); std::complex x((t - mu) / (LauConstants::root2 * sigma),0.); Double_t xRe = x.real(); Double_t z_H = ((gamma - deltaGamma_->unblindValue() / 2.) * sigma) / LauConstants::root2; Double_t z_L = ((gamma + deltaGamma_->unblindValue() / 2.) * sigma) / LauConstants::root2; //Doing H std::complex arg_H1(0., z_H - x.real()); std::complex arg_H2 = -(x*x) - (arg_H1 * arg_H1); std::complex conv_H = arg_H1.imag() < -5. ? (0.5 * std::exp(arg_H2)) * RooMath::erfc(-1i * arg_H1) : 0.5 * TMath::Exp(-( xRe * xRe )) * RooMath::faddeeva(arg_H1); //Doing L std::complex arg_L1(0., z_L - x.real()); std::complex arg_L2 = -(x*x) - (arg_L1 * arg_L1); std::complex conv_L = arg_L1.imag() < -5. ? (0.5 * std::exp(arg_L2)) * RooMath::erfc(-1i * arg_L1) : 0.5 * TMath::Exp(-( xRe * xRe )) * RooMath::faddeeva(arg_L1); std::complex cosh_conv = 0.5 * (conv_H + conv_L); std::complex sinh_conv = 0.5 * (conv_H - conv_L); return {cosh_conv.real(), sinh_conv.real()}; } Double_t LauDecayTimePdf::smearedExpTerm(Double_t t, Double_t sigma, Double_t mu) { using namespace std::complex_literals; const Double_t gamma = 1. / this->tau_->unblindValue(); const std::complex x((t - mu) / (LauConstants::root2 * sigma),0.); const Double_t xRe = x.real(); const Double_t z = (gamma * sigma) / LauConstants::root2; const std::complex arg1(0., z - x.real()); const std::complex arg2 = -(x * x) - (arg1 * arg1); const std::complex conv = arg1.imag() < -5. ? 0.5 * (std::exp(arg2)) * RooMath::erfc(-1i * arg1) : 0.5 * TMath::Exp(-(xRe * xRe)) * RooMath::faddeeva(arg1) ; return conv.real(); } std::pair LauDecayTimePdf::nonSmearedCosSinIntegral(Double_t minAbs, Double_t maxAbs) { // From 1407.0748, not clear whether complex is faster in this case Double_t gamma = 1. / this->tau_->unblindValue(); LauComplex denom = LauComplex(gamma, -this->deltaM_->unblindValue()); LauComplex exponent = LauComplex(-gamma, this->deltaM_->unblindValue()); LauComplex num0 = -exponent.scale(minAbs).exp(); LauComplex num1 = -exponent.scale(maxAbs).exp(); LauComplex integral = (num1 - num0) / denom; return {integral.re(), integral.im()}; } std::pair LauDecayTimePdf::smearedCosSinIntegral(Double_t minAbs, Double_t maxAbs, Double_t sigma, Double_t mu) { using namespace std::complex_literals; Double_t gamma = 1. / this->tau_->unblindValue(); Double_t x1 = (maxAbs - mu) / (LauConstants::root2 * sigma); Double_t x0 = (minAbs - mu) / (LauConstants::root2 * sigma); std::complex z = std::complex(gamma * sigma / LauConstants::root2, -this->deltaM_->unblindValue() * sigma / LauConstants::root2); std::complex arg1 = std::complex(0., 1.) * (z - x1); std::complex arg0 = std::complex(0., 1.) * (z - x0); std::complex integral = 0. + 0i; if(arg1.imag() < -5.) {integral = RooMath::erf(x1) - std::exp(-(x1 * x1) - (arg1 * arg1)) * RooMath::erfc(-1i * arg1);} else {integral = RooMath::erf(x1) - TMath::Exp(-(x1 * x1)) * RooMath::faddeeva(arg1);} if(arg0.imag() < -5.) {integral -= RooMath::erf(x0) - std::exp(-(x0 * x0) - (arg0 * arg0)) * RooMath::erfc(-1i * arg0);} else {integral -= RooMath::erf(x0) - TMath::Exp(-(x0 * x0)) * RooMath::faddeeva(arg0);} integral *= (sigma / (2. * LauConstants::root2 * z)); Double_t cos_integral = integral.real(); Double_t sin_integral = integral.imag(); return {cos_integral, sin_integral}; } Double_t LauDecayTimePdf::nonSmearedExpIntegral(Double_t minAbs, Double_t maxAbs) { const Double_t tau = tau_->unblindValue(); const Double_t Gamma = 1.0 / tau; return tau * ( TMath::Exp(-minAbs*Gamma) - TMath::Exp(-maxAbs*Gamma) ); } Double_t LauDecayTimePdf::smearedExpIntegral(Double_t minAbs, Double_t maxAbs, Double_t sigma, Double_t mu) { using namespace std::complex_literals; const Double_t gamma = 1. / this->tau_->unblindValue(); const Double_t x1 = (maxAbs - mu) / (LauConstants::root2 * sigma); const Double_t x0 = (minAbs - mu) / (LauConstants::root2 * sigma); const Double_t z = (gamma * sigma) / LauConstants::root2; std::complex arg1(0., z - x1); std::complex arg0(0., z - x0); std::complex integral = 0. + 0i; if(arg1.imag() < -5.) {integral = RooMath::erf(x1) - std::exp(-(x1 * x1) - (arg1 * arg1)) * RooMath::erfc(-1i * arg1);} else {integral = RooMath::erf(x1) - TMath::Exp(-(x1 * x1)) * RooMath::faddeeva(arg1);} if(arg0.imag() < -5.) {integral -= RooMath::erf(x0) - std::exp(-(x0 * x0) - (arg0 * arg0)) * RooMath::erfc(-1i * arg0);} else {integral -= RooMath::erf(x0) - TMath::Exp(-(x0 * x0)) * RooMath::faddeeva(arg0);} integral *= (sigma / (2. * LauConstants::root2 * z)); return integral.real(); } std::pair LauDecayTimePdf::nonSmearedCoshSinhIntegral(Double_t minAbs, Double_t maxAbs) { // Use exponential formualtion rather than cosh, sinh. // Fewer terms (reused for each), but not guaranteed to be faster. Double_t gamma = 1. / this->tau_->unblindValue(); Double_t gammaH = gamma - 0.5 * deltaGamma_->unblindValue(); Double_t gammaL = gamma - 0.5 * deltaGamma_->unblindValue(); Double_t nL1 = -TMath::Exp(-gammaL * maxAbs) / gammaL; Double_t nH1 = -TMath::Exp(-gammaH * maxAbs) / gammaH; Double_t nL0 = -TMath::Exp(-gammaL * minAbs) / gammaL; Double_t nH0 = -TMath::Exp(-gammaH * minAbs) / gammaH; Double_t cosh_integral = 0.5 * ( (nH1 + nL1) - (nH0 + nL0) ); Double_t sinh_integral = 0.5 * ( (nH1 - nL1) - (nH0 - nL0) ); return {cosh_integral, sinh_integral}; } std::pair LauDecayTimePdf::smearedCoshSinhIntegral(Double_t minAbs, Double_t maxAbs, Double_t sigma, Double_t mu) { using namespace std::complex_literals; Double_t gamma = 1. / this->tau_->unblindValue(); Double_t x1 = (maxAbs - mu) / (LauConstants::root2 * sigma); Double_t x0 = (minAbs - mu) / (LauConstants::root2 * sigma); Double_t z_H = ((gamma - deltaGamma_->unblindValue() / 2.) * sigma) / LauConstants::root2; std::complex arg1_H(0., z_H - x1); std::complex arg0_H(0., z_H - x0); std::complex integral_H = 0. + 0i; if(arg1_H.imag() < -5.) {integral_H = RooMath::erf(x1) - std::exp(-(x1 * x1) - (arg1_H * arg1_H)) * RooMath::erfc(-1i * arg1_H);} else {integral_H = RooMath::erf(x1) - TMath::Exp(-(x1 * x1)) * RooMath::faddeeva(arg1_H);} if(arg0_H.imag() < -5.) {integral_H -= RooMath::erf(x0) - std::exp(-(x0 * x0) - (arg0_H * arg0_H)) * RooMath::erfc(-1i * arg0_H);} else {integral_H -= RooMath::erf(x0) - TMath::Exp(-(x0 * x0)) * RooMath::faddeeva(arg0_H);} integral_H *= (sigma / (2. * LauConstants::root2 * z_H)); // Same for light (L) Double_t z_L = ((gamma + deltaGamma_->unblindValue() / 2.) * sigma) / LauConstants::root2; std::complex arg1_L(0., z_L - x1); std::complex arg0_L(0., z_L - x0); std::complex integral_L = 0. + 0i; if(arg1_L.imag() < -5.) {integral_L = RooMath::erf(x1) - std::exp(-(x1 * x1) - (arg1_L * arg1_L)) * RooMath::erfc(-1i * arg1_L);} else {integral_L = RooMath::erf(x1) - TMath::Exp(-(x1 * x1)) * RooMath::faddeeva(arg1_L);} if(arg0_L.imag() < -5.) {integral_L -= RooMath::erf(x0) - std::exp(-(x0 * x0) - (arg0_L * arg0_L)) * RooMath::erfc(-1i * arg0_L);} else {integral_L -= RooMath::erf(x0) - TMath::Exp(-(x0 * x0)) * RooMath::faddeeva(arg0_L);} integral_L *= (sigma / (2. * LauConstants::root2 * z_L)); std::complex cosh_integral = 0.5 * (integral_H + integral_L); std::complex sinh_integral = 0.5 * (integral_H - integral_L); return {cosh_integral.real(), sinh_integral.real()}; } void LauDecayTimePdf::calcNorm(const Double_t abscissaErr) { - if( abscissaErr <= 0. and scaleWithPerEventError_) + /*if( abscissaErr <= 0. and scaleWithPerEventError_) { std::cerr << "\033[1;31m IN CALCNORM: \33[0m" << std::endl; std::cerr << "\033[1;31m absErr: " << abscissaErr << "\033[0m" << std::endl; //DEBUG - } + }*/ // first reset integrals to zero normTermExp_ = 0.0; normTermCos_ = 0.0; normTermSin_ = 0.0; normTermCosh_ = 0.0; normTermSinh_ = 0.0; // Get all the up to date parameter values std::vector fracs(nGauss_); std::vector means(nGauss_); std::vector sigmas(nGauss_); // TODO - why do we do the fractions this way around? fracs[0] = 1.0; for (UInt_t i(0); iunblindValue(); sigmas[i] = sigma_[i]->unblindValue(); if (i != 0) { fracs[i] = frac_[i-1]->unblindValue(); fracs[0] -= fracs[i]; } } // Scale the gaussian parameters by the per-event error on decay time (if appropriate) for (UInt_t i(0); i doSmearing() ) {this->calcSmearedPartialIntegrals( minAbscissa_, maxAbscissa_ , 1.0, means, sigmas, fracs);} else {this->calcNonSmearedPartialIntegrals( minAbscissa_, maxAbscissa_, 1.0 );} break; case EfficiencyMethod::Binned : // Efficiency varies as piecewise constant // Total integral is sum of integrals in each bin, each weighted by efficiency in that bin for ( Int_t bin{1}; bin <= effiHist_->GetNbinsX(); ++bin ) { const Double_t loEdge {effiHist_->GetBinLowEdge(bin)}; const Double_t hiEdge {loEdge + effiHist_->GetBinWidth(bin)}; const Double_t effVal {effiHist_->GetBinContent(bin)}; if ( this -> doSmearing() ) {this->calcSmearedPartialIntegrals( loEdge, hiEdge, effVal, means, sigmas, fracs );} else {this->calcNonSmearedPartialIntegrals( loEdge, hiEdge, effVal );} } break; case EfficiencyMethod::Spline : // Efficiency varies as piecewise polynomial // TODO - to be worked out what to do here if(not effiFun_){std::cerr << "FATAL : no spline defined!"; gSystem->Exit(EXIT_FAILURE);} for(size_t i = 0; i < effiFun_ -> getnKnots()-1; ++i) { if( this -> doSmearing() ) {this -> calcSmearedSplinePartialIntegrals( i, means, sigmas, fracs );} else {this -> calcNonSmearedSplinePartialIntegrals( i );} } // std::cout << "\033[1;31m Normalisation values: \n" << normTermExp_ << "\n" << normTermCos_ << "\n" << normTermSin_ << "\n" << normTermCosh_ << "\n" << normTermSinh_ << "\n\n\033[0m" << std::endl; //DEBUG /* std::cerr << "WARNING in LauDecayTimePdf::calcNorm : normalisation integrals for spline acceptance not yet implemented - effect of acceptance will be neglected!" << std::endl; if ( this -> doSmearing() ) {this->calcSmearedPartialIntegrals( minAbscissa_, maxAbscissa_ , 1.0, mean, sigma, frac);} else {this->calcNonSmearedPartialIntegrals( minAbscissa_, maxAbscissa_, 1.0 );} */ break; } // TODO - should we check here that all terms we expect to use are now non-zero? +// std::cout << "\033[1;34m In calcNorm: \033[0m" << std::endl; //DEBUG +// std::cout << "\033[1;34m normTermExp : \033[0m" << normTermExp_ << std::endl; //DEBUG +// std::cout << "\033[1;34m normTerm[Cos,Sin] : \033[0m[" << normTermCos_ << ", " << normTermSin_ << "]" << std::endl; //DEBUG +// std::cout << "\033[1;34m normTerm[Cosh,Sinh]: \033[0m[" << normTermCosh_ << ", " << normTermSinh_ << "]" << std::endl; //DEBUG } // TODO - Mildly concerned this is void rather than returning the integrals // (but this would require refactoring for different return values). // As long as it doesn't get called outside of calcNorm() it should be fine - DPO void LauDecayTimePdf::calcNonSmearedPartialIntegrals(const Double_t minAbs, const Double_t maxAbs, const Double_t weight) { Double_t normTermExp {0.0}; if (method_ == DecayTime) { normTermExp = weight * this -> nonSmearedExpIntegral(minAbs, maxAbs); } else if (method_ == DecayTimeDiff) { const Double_t tau = tau_->unblindValue(); const Double_t Gamma = 1.0 / tau; // TODO - there should be some TMath::Abs here surely? normTermExp = weight * tau * (2.0 - TMath::Exp(-maxAbs*Gamma) - TMath::Exp(-minAbs*Gamma)); } normTermExp_ += normTermExp; // Normalisation factor for B0 decays if ( type_ == ExpTrig ) { normTermCosh_ += normTermExp; auto [cosIntegral, sinIntegral] = this->nonSmearedCosSinIntegral(minAbs, maxAbs); normTermCos_ += weight * cosIntegral; normTermSin_ += weight * sinIntegral; } // Normalisation factor for Bs decays else if ( type_ == ExpHypTrig ) { auto [coshIntegral, sinhIntegral] = this->nonSmearedCoshSinhIntegral(minAbs, maxAbs); normTermCosh_ += weight * coshIntegral; normTermSinh_ += weight * sinhIntegral; auto [cosIntegral, sinIntegral] = this->nonSmearedCosSinIntegral(minAbs, maxAbs); normTermCos_ += weight * cosIntegral; normTermSin_ += weight * sinIntegral; } - //std::cout << "\033[1;34m In calcNonSmearedPartialIntegrals: \033[0m" << std::endl; //DEBUG - //std::cout << "\033[1;34m normTermExp : \033[0m" << normTermExp_ << std::endl; //DEBUG - //std::cout << "\033[1;34m normTerm[Cos,Sin] : \033[0m[" << normTermCos_ << ", " << normTermSin_ << "]" << std::endl; //DEBUG - //std::cout << "\033[1;34m normTerm[Cosh,Sinh]: \033[0m[" << normTermCosh_ << ", " << normTermSinh_ << "]" << std::endl; //DEBUG } void LauDecayTimePdf::calcSmearedPartialIntegrals(const Double_t minAbs, const Double_t maxAbs, const Double_t weight, const std::vector& means, const std::vector& sigmas, const std::vector& fractions) { for (UInt_t i(0); i smearedExpIntegral(minAbs, maxAbs, sigmas[i], means[i]); } else if (method_ == DecayTimeDiff) { const Double_t tau = tau_->unblindValue(); const Double_t Gamma = 1.0 / tau; // TODO - this is neglecting resolution at the moment // TODO - there should be some TMath::Abs here surely? normTermExp = weight * tau * (2.0 - TMath::Exp(-maxAbs*Gamma) - TMath::Exp(-minAbs*Gamma)); } normTermExp_ += fractions[i] * normTermExp; // Normalisation factor for B0 decays if ( type_ == ExpTrig ) { normTermCosh_ += fractions[i] * normTermExp; auto [cosIntegral, sinIntegral] = this->smearedCosSinIntegral(minAbs, maxAbs, sigmas[i], means[i]); normTermCos_ += fractions[i] * weight * cosIntegral; normTermSin_ += fractions[i] * weight * sinIntegral; } // Normalisation factor for Bs decays else if ( type_ == ExpHypTrig ) { auto [coshIntegral, sinhIntegral] = this->smearedCoshSinhIntegral(minAbs, maxAbs, sigmas[i], means[i]); normTermCosh_ += fractions[i] * weight * coshIntegral; normTermSinh_ += fractions[i] * weight * sinhIntegral; auto [cosIntegral, sinIntegral] = this->smearedCosSinIntegral(minAbs, maxAbs, sigmas[i], means[i]); normTermCos_ += fractions[i] * weight * cosIntegral; normTermSin_ += fractions[i] * weight * sinIntegral; } } - //std::cout << "\033[1;34m In calcSmearedPartialIntegrals: \033[0m" << std::endl; //DEBUG - //std::cout << "\033[1;34m normTermExp : \033[0m" << normTermExp_ << std::endl; //DEBUG - //std::cout << "\033[1;34m normTerm[Cos,Sin] : \033[0m[" << normTermCos_ << ", " << normTermSin_ << "]" << std::endl; //DEBUG - //std::cout << "\033[1;34m normTerm[Cosh,Sinh]: \033[0m[" << normTermCosh_ << ", " << normTermSinh_ << "]" << std::endl; //DEBUG } void LauDecayTimePdf::calcSmearedSplinePartialIntegrals(const UInt_t splineIndex, const std::vector& means, const std::vector& sigmas, const std::vector& fractions) { using namespace std::complex_literals; const Double_t gamma = 1. / this->tau_->unblindValue(); for (UInt_t i(0); i z = (gamma * sigmas[i]) / LauConstants::root2; normTermExp = this -> smearedSplineNormalise(splineIndex, z, sigmas[i], means[i]).first; } else if (method_ == DecayTimeDiff) {//TODO this isn't implemented at all //const Double_t tau = tau_->unblindValue(); //const Double_t Gamma = 1.0 / tau; // TODO - this is neglecting resolution at the moment // TODO - there should be some TMath::Abs here surely? // normTermExp = tau * (2.0 - TMath::Exp(-maxAbs*Gamma) - TMath::Exp(-minAbs*Gamma)); // ; // nop so the compiler doesn't complain } normTermExp_ += fractions[i] * normTermExp; // Normalisation factor for B0 decays if ( type_ == ExpTrig ) { normTermCosh_ += fractions[i] * normTermExp; -// std::cout << "\033[1;31m exp: " << normTermExp << "\n\033[0m"; //DEBUG const std::complex z = std::complex(gamma * sigmas[i] / LauConstants::root2, -this->deltaM_->unblindValue() * sigmas[i] / LauConstants::root2); auto [cosIntegral, sinIntegral] = this -> smearedSplineNormalise(splineIndex, z, sigmas[i], means[i]); normTermCos_ += fractions[i] * cosIntegral; normTermSin_ += fractions[i] * sinIntegral; -// std::cout << "\033[1;31m cos: " << cosIntegral << "\n\033[0m"; //DEBUG -// std::cout << "\033[1;31m sin: " << sinIntegral << "\n\033[0m"; //DEBUG } // Normalisation factor for Bs decays else if ( type_ == ExpHypTrig ) { const std::complex z_H = ((gamma - deltaGamma_->unblindValue() / 2.) * sigmas[i]) / LauConstants::root2; const std::complex z_L = ((gamma + deltaGamma_->unblindValue() / 2.) * sigmas[i]) / LauConstants::root2; const Double_t N_H = this -> smearedSplineNormalise(splineIndex, z_H, sigmas[i], means[i]).first; const Double_t N_L = this -> smearedSplineNormalise(splineIndex, z_L, sigmas[i], means[i]).first; const Double_t coshIntegral = 0.5 * (N_H + N_L); const Double_t sinhIntegral = 0.5 * (N_H - N_L); -// std::cout << "\033[1;31m cosh: " << coshIntegral << "\n\033[0m"; //DEBUG -// std::cout << "\033[1;31m sinh: " << sinhIntegral << "\n\033[0m"; //DEBUG - normTermCosh_ += fractions[i] * coshIntegral; normTermSinh_ += fractions[i] * sinhIntegral; const std::complex z = std::complex(gamma * sigmas[i] / LauConstants::root2, -this->deltaM_->unblindValue() * sigmas[i] / LauConstants::root2); auto [cosIntegral, sinIntegral] = this -> smearedSplineNormalise(splineIndex, z, sigmas[i], means[i]); normTermCos_ += fractions[i] * cosIntegral; normTermSin_ += fractions[i] * sinIntegral; -// std::cout << "\033[1;31m cos: " << cosIntegral << "\n\033[0m"; //DEBUG -// std::cout << "\033[1;31m sin: " << sinIntegral << "\n\033[0m"; //DEBUG } } - //std::cout << "\033[1;34m In calcSmearedSplinePartialIntegrals: \033[0m" << std::endl; //DEBUG - //std::cout << "\033[1;34m normTermExp : \033[0m" << normTermExp_ << std::endl; //DEBUG - //std::cout << "\033[1;34m normTerm[Cos,Sin] : \033[0m[" << normTermCos_ << ", " << normTermSin_ << "]" << std::endl; //DEBUG - //std::cout << "\033[1;34m normTerm[Cosh,Sinh]: \033[0m[" << normTermCosh_ << ", " << normTermSinh_ << "]" << std::endl; //DEBUG } std::array,4> LauDecayTimePdf::generateKvector(const std::complex z) { std::array,4> K = {0.,0.,0.,0.}; const std::complex zr = 1./z; K[0] = 0.5*zr; K[1] = 0.5*zr*zr; K[2] = zr*(1.+zr*zr); K[3] = 3.*zr*zr*(1.+zr*zr); return K; } std::array,4> LauDecayTimePdf::generateMvector(const Double_t minAbs, const Double_t maxAbs, const std::complex z, const Double_t sigma, const Double_t mu) { using namespace std::complex_literals; std::array,4> M0 = {0.,0.,0.,0.}; std::array,4> M1 = {0.,0.,0.,0.}; std::array,4> M; const Double_t x1 = (maxAbs - mu) / (LauConstants::root2 * sigma); const Double_t x0 = (minAbs - mu) / (LauConstants::root2 * sigma); +// std::cout << "\033[1;31m x1 :" << x1 << ", x0:" << x0 <<" \33[0m\t"; //DEBUG //Values used a lot const Double_t ex2_1 = TMath::Exp(-(x1*x1)); const Double_t ex2_0 = TMath::Exp(-(x0*x0)); const Double_t sqrtPir = 1/LauConstants::rootPi; const std::complex arg1 = (0.+1.i) * (z-x1); const std::complex arg0 = (0.+1.i) * (z-x0); //fad = the faddeeva term times the ex2 value (done in different ways depending on the domain) std::complex fad1; std::complex fad0; if(arg1.imag() < -5.) {fad1 = std::exp(-(x1 * x1) - (arg1 * arg1)) * RooMath::erfc(-1i * arg1);} else {fad1 = ex2_1*RooMath::faddeeva(arg1);} if(arg0.imag() < -5.) {fad0 = std::exp(-(x0 * x0) - (arg0 * arg0)) * RooMath::erfc(-1i * arg0);} else {fad0 = ex2_0*RooMath::faddeeva(arg0);} //doing the actual functions for x1 M1[0] = RooMath::erf(x1) - fad1; M1[1] = 2. * (-sqrtPir*ex2_1 - x1*fad1); - M1[2] = 2. * (-2*x1*sqrtPir*ex2_1 - (2*x1 - 1)*fad1); + M1[2] = 2. * (-2*x1*sqrtPir*ex2_1 - (2*x1*x1 - 1)*fad1); M1[3] = 4. * (-(2*x1*x1 - 1)*sqrtPir*ex2_1 - x1*(2*x1*x1-3)*fad1); //doing them again for x0 M0[0] = RooMath::erf(x0) - fad0; M0[1] = 2. * (-sqrtPir*ex2_0 - x0*fad0); - M0[2] = 2. * (-2*x0*sqrtPir*ex2_0 - (2*x0 - 1)*fad0); + M0[2] = 2. * (-2*x0*sqrtPir*ex2_0 - (2*x0*x0 - 1)*fad0); M0[3] = 4. * (-(2*x0*x0 - 1)*sqrtPir*ex2_0 - x0*(2*x0*x0-3)*fad0); for(Int_t i = 0; i < 4; ++i){M[i] = M1[i] - M0[i];} return M; } std::pair LauDecayTimePdf::smearedSplineNormalise(const UInt_t splineIndex, std::complex z, const Double_t sigma, const Double_t mu) { using namespace std::complex_literals; const std::vector& xVals = effiFun_ -> getXValues(); const Double_t minAbs = xVals[splineIndex]; const Double_t maxAbs = xVals[splineIndex+1]; std::array coeffs = effiFun_ -> getCoefficients(splineIndex); std::array,4> K = this -> generateKvector(z); // std::cout << "\033[1;31m zVal: "<< z<<" Kvector: ["; for (std::complex k : K) std::cout << k << ", "; std::cout << "\b\b] \33[0m\t"; //DEBUG std::array,4> M = this -> generateMvector(minAbs, maxAbs, z, sigma, mu); // std::cout << "\033[1;31m Mvector: ["; for (std::complex m : M) std::cout << m << ", "; std::cout << "\b\b] \33[0m" << std::endl; //DEBUG +// gSystem -> Exit(EXIT_SUCCESS); //Double sum to get N (eqn 31 in https://arxiv.org/pdf/1407.0748.pdf) std::complex N = 0. + 0i; for(Int_t i = 0; i < 4; ++i) { for(Int_t j = 0; j <= i ; ++j) { Int_t ij = i+j; if(ij > 3){continue;} // sum component = 0 Double_t A = coeffs[ij] * TMath::Binomial(ij,j) / TMath::Power(2,ij); N += A * M[i] * K[j]; } } return std::make_pair( N.real(), N.imag() ); } void LauDecayTimePdf::calcNonSmearedSplinePartialIntegrals(const UInt_t splineIndex) { using namespace std::complex_literals; const Double_t gamma = 1. / this->tau_->unblindValue(); Double_t normTermExp {0.0}; if (method_ == DecayTime) { const std::complex u = gamma; normTermExp = this -> nonSmearedSplineNormalise(splineIndex, u).first; } else if (method_ == DecayTimeDiff) {//TODO this isn't implemented at all //const Double_t tau = tau_->unblindValue(); //const Double_t Gamma = 1.0 / tau; // TODO - this is neglecting resolution at the moment // TODO - there should be some TMath::Abs here surely? // normTermExp = tau * (2.0 - TMath::Exp(-maxAbs*Gamma) - TMath::Exp(-minAbs*Gamma)); // ; // nop so the compiler doesn't complain } normTermExp_ += normTermExp; // Normalisation factor for B0 decays if ( type_ == ExpTrig ) { normTermCosh_ += normTermExp; - // std::cout << "\033[1;31m exp: " << normTermExp << "\n\033[0m"; //DEBUG const std::complex u = std::complex(gamma, -this->deltaM_->unblindValue()); auto [cosIntegral, sinIntegral] = this -> nonSmearedSplineNormalise(splineIndex, u); normTermCos_ += cosIntegral; normTermSin_ += sinIntegral; - - // std::cout << "\033[1;31m cos: " << cosIntegral << "\n\033[0m"; //DEBUG - // std::cout << "\033[1;31m sin: " << sinIntegral << "\n\033[0m"; //DEBUG } // Normalisation factor for Bs decays - //TODO this uses u not z, which will be the same for both u_H and u_L? coshIntegral = cosIntegral and sinhIntegral = 0 ? else if ( type_ == ExpHypTrig ) { -/* - const std::complex z_H = ((gamma - deltaGamma_->unblindValue() / 2.) * sigmas[i]) / LauConstants::root2; - const std::complex z_L = ((gamma + deltaGamma_->unblindValue() / 2.) * sigmas[i]) / LauConstants::root2; - const Double_t N_H = this -> smearedSplineNormalise(splineIndex, z_H, sigmas[i], means[i]).first; - const Double_t N_L = this -> smearedSplineNormalise(splineIndex, z_L, sigmas[i], means[i]).first; + const std::complex u_H = (gamma - deltaGamma_->unblindValue() / 2.); + const std::complex u_L = (gamma + deltaGamma_->unblindValue() / 2.); + + const Double_t N_H = this -> nonSmearedSplineNormalise(splineIndex, u_H).first; + const Double_t N_L = this -> nonSmearedSplineNormalise(splineIndex, u_L).first; const Double_t coshIntegral = 0.5 * (N_H + N_L); const Double_t sinhIntegral = 0.5 * (N_H - N_L); - // std::cout << "\033[1;31m cosh: " << coshIntegral << "\n\033[0m"; //DEBUG - // std::cout << "\033[1;31m sinh: " << sinhIntegral << "\n\033[0m"; //DEBUG + normTermCosh_ += coshIntegral; + normTermSinh_ += sinhIntegral; - normTermCosh_ += fractions[i] * coshIntegral; - normTermSinh_ += fractions[i] * sinhIntegral; - - const std::complex z = std::complex(gamma * sigmas[i] / LauConstants::root2, -this->deltaM_->unblindValue() * sigmas[i] / LauConstants::root2); - auto [cosIntegral, sinIntegral] = this -> smearedSplineNormalise(splineIndex, z, sigmas[i], means[i]); - normTermCos_ += fractions[i] * cosIntegral; - normTermSin_ += fractions[i] * sinIntegral; - // std::cout << "\033[1;31m cos: " << cosIntegral << "\n\033[0m"; //DEBUG - // std::cout << "\033[1;31m sin: " << sinIntegral << "\n\033[0m"; //DEBUG - */ - std::cerr << "ERROR: cosh and sinh terms still being debugged!" << std::endl; + const std::complex u = std::complex(gamma, -this->deltaM_->unblindValue()); + auto [cosIntegral, sinIntegral] = this -> nonSmearedSplineNormalise(splineIndex, u); + + normTermCos_ += cosIntegral; + normTermSin_ += sinIntegral; } } std::complex LauDecayTimePdf::I_k(const Int_t k, const Double_t minAbs, const Double_t maxAbs, const std::complex u /*= Gamma - iDeltam*/) { //Taking mu = 0, this does not have to be the case in general auto G = [&u](const Int_t n){return -TMath::Factorial(n)/std::pow(u,n+1);};//power of n+1 used rather than n, this is due to maths error in the paper auto H = [&u](const Int_t n, const Double_t t){return std::pow(t,n)*std::exp(-u*t);}; std::complex ans = 0; for (Int_t j = 0; j <= k; ++j) {ans += TMath::Binomial(k,j)*G(j)*( H( k-j, maxAbs ) - H( k-j, minAbs ) );} return ans; } std::pair LauDecayTimePdf::nonSmearedSplineNormalise(const UInt_t splineIndex, std::complex u /*= Gamma - iDeltam*/) { using namespace std::complex_literals; const std::vector& xVals = effiFun_ -> getXValues(); const Double_t minAbs = xVals[splineIndex]; const Double_t maxAbs = xVals[splineIndex+1]; std::array coeffs = effiFun_ -> getCoefficients(splineIndex); //sum to get N (eqn 30 in https://arxiv.org/pdf/1407.0748.pdf, using I_k from Appendix B.1 with the corrected maths error) std::complex N = 0. + 0i; for(Int_t i = 0; i < 4; ++i){N += I_k(i, minAbs, maxAbs, u) * coeffs[i];} return std::make_pair( N.real(), N.imag() ); } Double_t LauDecayTimePdf::generateError(Bool_t forceNew) { if (errHist_ && (forceNew || !abscissaErrorGenerated_)) { LauFitData errData = errHist_->generate(nullptr); abscissaError_ = errData.at(this->varErrName()); abscissaErrorGenerated_ = kTRUE; } else { while (forceNew || !abscissaErrorGenerated_) { abscissaError_ = LauRandom::randomFun()->Landau(errorDistMPV_,errorDistSigma_); if (abscissaError_ < maxAbscissaError_ && abscissaError_ > minAbscissaError_) { abscissaErrorGenerated_ = kTRUE; forceNew = kFALSE; } } } return abscissaError_; } /* LauFitData LauDecayTimePdf::generate(const LauKinematics* kinematics) { // generateError SHOULD have been called before this // function but will call it here just to make sure // (has ns effect if has already been called) abscissaError_ = this->generateError(); // If the PDF is scaled by the per-event error then need to update the PDF height for each event Bool_t scale(kFALSE); for (std::vector::const_iterator iter = scaleMeans_.begin(); iter != scaleMeans_.end(); ++iter) { scale |= (*iter); } for (std::vector::const_iterator iter = scaleWidths_.begin(); iter != scaleWidths_.end(); ++iter) { scale |= (*iter); } if (scale || (!this->heightUpToDate() && !this->cachePDF())) { this->calcPDFHeight(kinematics); this->heightUpToDate(kTRUE); } // Generate the value of the abscissa. Bool_t gotAbscissa(kFALSE); Double_t genVal(0.0); Double_t genPDFVal(0.0); LauFitData genAbscissa; const Double_t xMin = this->minAbscissa(); const Double_t xMax = this->maxAbscissa(); const Double_t xRange = xMax - xMin; while (!gotAbscissa) { genVal = LauRandom::randomFun()->Rndm()*xRange + xMin; this->calcLikelihoodInfo(genVal, abscissaError_); genPDFVal = this->getUnNormLikelihood(); if (LauRandom::randomFun()->Rndm() <= genPDFVal/this->getMaxHeight()) {gotAbscissa = kTRUE;} if (genPDFVal > this->getMaxHeight()) { std::cerr<<"Warning in LauDecayTimePdf::generate()." <<" genPDFVal = "<getMaxHeight()<<" for the abscissa = "<varName()] = genVal; // mark that we need a new error to be generated next time abscissaErrorGenerated_ = kFALSE; return genAbscissa; } */ void LauDecayTimePdf::setErrorHisto(const TH1* hist) { if ( errHist_ != nullptr ) { std::cerr<<"WARNING in LauDecayTimePdf::setErrorHisto : Error histogram already set, not doing it again."<varErrName(), hist, this->minAbscissaError(), this->maxAbscissaError()); } void LauDecayTimePdf::setHistoPdf(const TH1* hist) { if ( pdfHist_ != nullptr ) { std::cerr<<"WARNING in LauDecayTimePdf::setHistoPdf : PDF histogram already set, not doing it again."<varName(), hist, this->minAbscissa(), this->maxAbscissa()); } void LauDecayTimePdf::setEffiHist(const TH1* hist) { if ( effiHist_ != nullptr ) { std::cerr << "WARNING in LauDecayTimePdf::setEffiHist : efficiency histogram already set, not doing it again." << std::endl; return; } if ( hist == nullptr ) { std::cerr << "WARNING in LauDecayTimePdf::setEffiHist : supplied efficiency histogram pointer is null." << std::endl; return; } // Check boundaries of histogram align with our abscissa's range const Double_t axisMin {hist->GetXaxis()->GetXmin()}; const Double_t axisMax {hist->GetXaxis()->GetXmax()}; if ( TMath::Abs(minAbscissa_ - axisMin)>1e-6 || TMath::Abs(maxAbscissa_ - axisMax)>1e-6 ) { std::cerr << "WARNING in LauDecayTimePdf::setEffiHist : mismatch in range between supplied histogram and abscissa\n" << " : histogram range: " << axisMin << " - " << axisMax << "\n" << " : abscissa range: " << minAbscissa_ << " - " << maxAbscissa_ << "\n" << " : Disregarding this histogram." << std::endl; return; } effiHist_ = dynamic_cast( hist->Clone() ); //Normalise the hist if the (relative) efficiencies have very large values if(effiHist_ -> GetMaximum() > 1.) { effiHist_ -> Scale( 1. / effiHist_->Integral() ); //Normalise std::cout << "INFO in LauDecayTimePdf::setEffiHist : Supplied histogram for Decay Time Acceptance has values too large: normalising..." << std::endl; } } void LauDecayTimePdf::setEffiSpline(Lau1DCubicSpline* spline) { if ( effiFun_ != 0 ) { std::cerr<<"WARNING in LauDecayTimePdf::setEffiPdf : efficiency function already set, not doing it again."< effis = effiFun_->getYValues(); effiPars_.resize( effis.size() ); size_t index = 0; for( Double_t& effi : effis ) { effiPars_[ index ] = new LauParameter( Form( "%s_Knot_%lu", varName_.Data() ,index ), effi, 0.0, 1.0, kTRUE ); ++index; } } LauAbsRValue* LauDecayTimePdf::findParameter(const TString& parName) { for ( std::vector::iterator iter = param_.begin(); iter != param_.end(); ++iter ) { if ((*iter)->name().Contains(parName)) { return (*iter); } } std::cerr << "ERROR in LauDecayTimePdf::findParameter : Parameter \"" << parName << "\" not found." << std::endl; return 0; } const LauAbsRValue* LauDecayTimePdf::findParameter(const TString& parName) const { for ( std::vector::const_iterator iter = param_.begin(); iter != param_.end(); ++iter ) { if ((*iter)->name().Contains(parName)) { return (*iter); } } std::cerr << "ERROR in LauDecayTimePdf::findParameter : Parameter \"" << parName << "\" not found." << std::endl; return 0; } void LauDecayTimePdf::updatePulls() { for ( std::vector::iterator iter = param_.begin(); iter != param_.end(); ++iter ) { std::vector params = (*iter)->getPars(); for (std::vector::iterator params_iter = params.begin(); params_iter != params.end(); ++params_iter ) { if (!(*iter)->fixed()) { (*params_iter)->updatePull(); } } } } void LauDecayTimePdf::updateEffiSpline(std::vector effiPars) { if (effiPars.size() != effiFun_->getnKnots()){ std::cerr<<"ERROR in LauDecayTimePdf::updateEffiSpline : number of efficiency parameters is not equal to the number of spline knots."<Exit(EXIT_FAILURE); } effiFun_->updateYValues(effiPars); }