diff --git a/src/LauDecayTimePdf.cc b/src/LauDecayTimePdf.cc index fa2f274..7581104 100644 --- a/src/LauDecayTimePdf.cc +++ b/src/LauDecayTimePdf.cc @@ -1,1677 +1,1563 @@ /* Copyright 2006 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 LauDecayTimePdf.cc \brief File containing implementation of LauDecayTimePdf class. */ #include #include using std::cout; using std::cerr; using std::endl; #include using std::complex; #include "TMath.h" #include "TRandom.h" #include "TSystem.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, TimeMeasurementMethod method) : 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), scaleMeans_(scale), scaleWidths_(scale), 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), pdfTerm_(0.0), effiTerm_(0.0), state_(Good), errHist_(0), pdfHist_(0), effiFun_(0), effiPars_(0), AProd_(0), ARaw_(0), flavTag_(0), asymTerm_(0), omegaTerm_(0), omegaTermh_(0) { this->initialise(); // Calculate the integrals of the decay time independent of the t this->storeIntegralTimeDep(minAbscissa_, maxAbscissa_); } 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, TimeMeasurementMethod method) : 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), scaleMeans_(scaleMeans), scaleWidths_(scaleWidths), 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), pdfTerm_(0.0), effiTerm_(0.0), state_(Good), errHist_(0), pdfHist_(0), effiFun_(0), effiPars_(0), AProd_(0), ARaw_(0), flavTag_(0), asymTerm_(0), omegaTerm_(0), omegaTermh_(0) { this->initialise(); // Calculate the integrals of the decay time independent of the t this->storeIntegralTimeDep(minAbscissa_, maxAbscissa_); } LauDecayTimePdf::~LauDecayTimePdf() { // Destructor delete errHist_; errHist_ = 0; delete pdfHist_; pdfHist_ = 0; delete effiFun_; effiFun_ = 0; } 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()) { cerr<<"ERROR in LauDecayTimePdf::initialise : scale vector size not the same as nGauss."<Exit(EXIT_FAILURE); } if (type_ == Hist){ if (this->nParameters() != 0){ 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)) { 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)) { cerr<<"ERROR in LauDecayTimePdf::initialise : Exp 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)) { 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)) { cerr<<"ERROR in LauDecayTimePdf::initialise : ExpHypTrig 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)) { cerr<<"ERROR in LauDecayTimePdf::initialise : DeltaExp type PDF requires:"<Exit(EXIT_FAILURE); } } else if (type_ == SimFitNormBd) { tau_ = this->findParameter("tau"); deltaM_ = this->findParameter("deltaM"); foundParams &= (tau_ != 0); foundParams &= (deltaM_ != 0); if ((this->nParameters() != (3*nGauss_-1+2)) || (!foundParams)) { cerr<<"ERROR in LauDecayTimePdf::initialise : SimFitNormBd type PDF requires:"<Exit(EXIT_FAILURE); } } else if (type_ == SimFitSigBd) { tau_ = this->findParameter("tau"); deltaM_ = this->findParameter("deltaM"); foundParams &= (tau_ != 0); foundParams &= (deltaM_ != 0); if ((this->nParameters() != (3*nGauss_-1+2)) || (!foundParams)) { cerr<<"ERROR in LauDecayTimePdf::initialise : SimFitSigBd type PDF requires:"<Exit(EXIT_FAILURE); } } else if (type_ == SimFitNormBs) { 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)) { cerr<<"ERROR in LauDecayTimePdf::initialise : SimFitNormBs type PDF requires:"<Exit(EXIT_FAILURE); } } else if (type_ == SimFitSigBs) { 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)) { cerr<<"ERROR in LauDecayTimePdf::initialise : SimFitSigBs type PDF requires:"<Exit(EXIT_FAILURE); } } } // Cache the normalisation factor //this->calcNorm(); } void LauDecayTimePdf::cacheInfo(const LauFitDataTree& inputData) { Bool_t hasBranch = inputData.haveBranch(this->varName()); if (!hasBranch) { cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Input data does not contain variable \""<varName()<<"\"."<varErrName()); if (!hasBranch) { 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); normTermsCosh_.clear(); normTermsCosh_.reserve(nEvents); normTermsSinh_.clear(); normTermsSinh_.reserve(nEvents); effiTerms_.clear(); effiTerms_.reserve(nEvents); for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) { const LauFitData& dataValues = inputData.getData(iEvt); LauFitData::const_iterator iter = dataValues.find(this->varName()); const Double_t abscissa = iter->second; if (abscissa > this->maxAbscissa() || abscissa < this->minAbscissa()) { cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Given value of the decay time: "<minAbscissa()<<","<maxAbscissa()<<"]."<Exit(EXIT_FAILURE); } abscissas_.push_back( abscissa ); iter = dataValues.find(this->varErrName()); Double_t abscissaErr = iter->second; if (abscissaErr > this->maxAbscissaError() || abscissaErr < this->minAbscissaError()) { cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Given value of the decay-time error: "<minAbscissaError()<<","<maxAbscissaError()<<"]."<Exit(EXIT_FAILURE); } abscissaErrors_.push_back(abscissaErr); this->calcLikelihoodInfo(abscissa, abscissaErr, iEvt); expTerms_.push_back(expTerm_); cosTerms_.push_back(cosTerm_); sinTerms_.push_back(sinTerm_); coshTerms_.push_back(coshTerm_); sinhTerms_.push_back(sinhTerm_); normTermsExp_.push_back(normTermExp_); normTermsCosh_.push_back(normTermCosh_); normTermsSinh_.push_back(normTermSinh_); if(effiFun_) effiTerms_.push_back(effiFun_->evaluate(abscissa)); else effiTerms_.push_back(1.0); } } } void LauDecayTimePdf::calcLikelihoodInfo(UInt_t iEvt) { 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]; normTermCosh_ = normTermsCosh_[iEvt]; normTermSinh_ = normTermsSinh_[iEvt]; } if ( errHist_ ) { errHist_->calcLikelihoodInfo(iEvt); errTerm_ = errHist_->getLikelihood(); } else { errTerm_ = 1.0; } const Double_t abscissa = abscissas_[iEvt]; //Parameters will change in some cases update things! if (type_ == SimFitNormBd || type_ == SimFitSigBd || type_ == SimFitNormBs || type_ == SimFitSigBs){ const Double_t abscissaErr = abscissaErrors_[iEvt]; this->calcLikelihoodInfo(abscissa,abscissaErr,iEvt); } if ( effiFun_ ) { this->updateEffiSpline(effiPars_); effiTerm_ = effiFun_->evaluate(abscissa); if (effiTerm_>1.0){effiTerm_=1.0;} if (effiTerm_<0.0){effiTerm_=0.0;} } else { effiTerm_ = 1.0; } // TODO need a check in here that none of the floating parameter values have changed // If they have, then we need to recalculate all or some of the terms /* if ( parsChanged ) { const Double_t abscissa = abscissas_[iEvt][0]; const Double_t abscissaErr = abscissaErrors_[iEvt]; this->calcLikelihoodInfo(abscissa, abscissaErr); } */ } void LauDecayTimePdf::calcLikelihoodInfo(Double_t abscissa, UInt_t iEvt) { // Check whether any of the gaussians should be scaled - if any of them should we need the per-event error 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) { cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfo : Per-event error on Delta t not provided, cannot calculate anything."<calcLikelihoodInfo(abscissa, 0.0,iEvt); } } void LauDecayTimePdf::calcLikelihoodInfo(Double_t abscissa, Double_t abscissaErr, UInt_t iEvt) { if (abscissa > this->maxAbscissa() || abscissa < this->minAbscissa()) { cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfo : Given value of the decay time: "<minAbscissa()<<","<maxAbscissa()<<"]."<Exit(EXIT_FAILURE); } if (abscissaErr > this->maxAbscissaError() || abscissaErr < this->minAbscissaError()) { cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfo : Given value of Delta t error: "<minAbscissaError()<<","<maxAbscissaError()<<"]."<Exit(EXIT_FAILURE); } // Initialise the various terms to zero if (type_ == Hist){ if ( pdfHist_ ) { pdfHist_->calcLikelihoodInfo(abscissa); pdfTerm_ = pdfHist_->getLikelihood(); } else { pdfTerm_ = 1.0; } }else{ // Reset the state to Good this->state(Good); // If we're not using the resolution function calculate the simple terms and return if (!this->doSmearing()) { this->calcNonSmearedTerms(abscissa, iEvt); return; } //TODO how much to be added below for SimFitNormBd/SimFitNormBs/SimFitSigBd/SimFitSigBs // Get all the up to date parameter values std::vector frac(nGauss_); std::vector mean(nGauss_); std::vector sigma(nGauss_); Double_t tau(0.0); Double_t deltaM(0.0); Double_t fracPrompt(0.0); Double_t Delta_gamma(0.0); frac[0] = 1.0; for (UInt_t i(0); ivalue(); sigma[i] = sigma_[i]->value(); if (i != 0) { frac[i] = frac_[i-1]->value(); frac[0] -= frac[i]; } } if (type_ != Delta) { tau = tau_->value(); if (type_ == ExpTrig) { deltaM = deltaM_->value(); } if (type_ == DeltaExp) { fracPrompt = fracPrompt_->value(); } if (type_ == ExpHypTrig){ deltaM = deltaM_->value(); Delta_gamma = deltaGamma_->value(); } } // 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) { std::vector expTerms(nGauss_); std::vector cosTerms(nGauss_); std::vector sinTerms(nGauss_); std::vector coshTerms(nGauss_); std::vector sinhTerms(nGauss_); std::vector expTermsNorm(nGauss_); // TODO - TEL changed this name to make it compile - please check! std::vector SinhTermsNorm(nGauss_); // Calculate values of the PDF convoluated with each Gaussian for a given value of the abscsissa for (UInt_t i(0); icalcTrigExponent(deltaM, tau, x[i], sigma[i], exponentTermRe, exponentTermIm); // Elements related to the trigonometric function, i.e. convolution of Exp*Sin or Cos with Gauss Double_t sinTrigTermRe, sinTrigTermIm, cosTrigTermRe, cosTrigTermIm; this->calcTrigConv(deltaM, tau, x[i], sigma[i], sinTrigTermRe, sinTrigTermIm, kFALSE); this->calcTrigConv(deltaM, tau, x[i], sigma[i], cosTrigTermRe, cosTrigTermIm, kTRUE); // Combining elements of the full pdf LauComplex zExp(exponentTermRe, exponentTermIm); LauComplex zTrigSin(sinTrigTermRe, sinTrigTermIm); LauComplex zTrigCos(cosTrigTermRe, cosTrigTermIm); LauComplex sinConv = zExp * zTrigSin; LauComplex cosConv = zExp * zTrigCos; sinConv.scale(1.0/4.0); cosConv.scale(1.0/4.0); // Cosine*Exp and Sine*Exp terms cosTerms[i] = cosConv.re(); sinTerms[i] = sinConv.im(); // Normalisation Double_t umax = xMax - mean[i]; Double_t umin = xMin - mean[i]; expTermsNorm[i] = (1.0/2.0) * tau * (-1.0 + TMath::Erf(umax/(LauConstants::root2 * sigma[i])) + TMath::Erfc(umin/(LauConstants::root2 * sigma[i])) + TMath::Exp((pow(sigma[i], 2) - 2.0 * tau * (xMax + xMin - mean[i]))/(2.0 * pow(tau, 2))) * (TMath::Exp(xMax/tau) * TMath::Erfc((pow(sigma[i], 2) - xMin)/(LauConstants::root2 * tau))) + (TMath::Exp(xMin/tau) * TMath::Erfc((pow(sigma[i], 2) - xMax)/(LauConstants::root2 * tau)))); } else { } } // Typical case (2): B0s/B0sbar if (type_ == ExpHypTrig) { // LHCb convention if (method_ == DecayTime) { // Convolution of Exp*cosh (Exp*sinh) with a gaussian //Double_t OverallExpFactor = 0.25*TMath::Exp(-(x[i]-mean[i])*(x[i]-mean[i])/(2*sigma[i]*sigma[i])); //Double_t ExpFirstTerm = TMath::Exp((2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))*(2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))/(8*sigma[i]*sigma[i]*tau*tau)); //Double_t ExpSecondTerm = TMath::Exp((2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))*(2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))/(8*sigma[i]*sigma[i]*tau*tau)); //Double_t ErfFirstTerm = TMath::Erf((2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); //Double_t ErfSecondTerm = TMath::Erf((2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); //Double_t sinhConv = OverallExpFactor*(ExpFirstTerm*(1+ErfFirstTerm) + ExpSecondTerm*(-1+ErfSecondTerm)); //Double_t coshConv = OverallExpFactor*(ExpFirstTerm*(1+ErfFirstTerm) - ExpSecondTerm*(-1+ErfSecondTerm)); //cosTerms[i] = sinhConv; // sinTerms[i] = coshConv; //TODO: check this formula and try to simplify it! double OverallExpTerm_max = (1/(2*(-4 + Delta_gamma*Delta_gamma*tau*tau)))*tau*TMath::Exp(-0.5*Delta_gamma*(xMax + mean[i]) - xMax/tau); double ErfTerm_max = -2*Delta_gamma*tau*TMath::Exp(0.5*Delta_gamma*(xMax+mean[i])+xMax/tau)*TMath::Erf((xMax-mean[i])/(TMath::Sqrt(2)*sigma[i])); double ExpFirstTerm_max = TMath::Exp(xMax*Delta_gamma+(sigma[i]*sigma[i]*(-2 + Delta_gamma*tau)*(-2 + Delta_gamma*tau))/(8*tau*tau)); double ErfcFirstTerm_max = TMath::Erfc((2*(-xMax + mean[i])*tau + sigma[i]*sigma[i]*(2 - Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); double ExpSecondTerm_max = TMath::Exp(Delta_gamma*mean[i] + (sigma[i]*sigma[i]*(2 + Delta_gamma*tau)*(2 + Delta_gamma*tau))/(8*tau*tau)); double ErfcSecondTerm_max = TMath::Erfc((2*(-xMax + mean[i])*tau + sigma[i]*sigma[i]*(2 + Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); double MaxVal= OverallExpTerm_max*(ErfTerm_max + TMath::Exp(mean[i]/tau)*(ExpFirstTerm_max*(2+Delta_gamma*tau)* ErfcFirstTerm_max + ExpSecondTerm_max*(-2+Delta_gamma*tau)* ErfcSecondTerm_max)); double OverallExpTerm_min = (1/(2*(-4 + Delta_gamma*Delta_gamma*tau*tau)))*tau*TMath::Exp(-0.5*Delta_gamma*(xMin + mean[i]) - xMin/tau); double ErfTerm_min = -2*Delta_gamma*tau*TMath::Exp(0.5*Delta_gamma*(xMin+mean[i])+xMin/tau)*TMath::Erf((xMin-mean[i])/(TMath::Sqrt(2)*sigma[i])); double ExpFirstTerm_min = TMath::Exp(xMin*Delta_gamma+(sigma[i]*sigma[i]*(-2 + Delta_gamma*tau)*(-2 + Delta_gamma*tau))/(8*tau*tau)); double ErfcFirstTerm_min = TMath::Erfc((2*(-xMin + mean[i])*tau + sigma[i]*sigma[i]*(2 - Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); // TODO - TEL added this (currently identical to ExpSecondTerm_max) to get this to compile - please check!! double ExpSecondTerm_min = TMath::Exp(Delta_gamma*mean[i] + (sigma[i]*sigma[i]*(2 + Delta_gamma*tau)*(2 + Delta_gamma*tau))/(8*tau*tau)); double ErfcSecondTerm_min = TMath::Erfc((2*(-xMin + mean[i])*tau + sigma[i]*sigma[i]*(2 + Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); double minVal= OverallExpTerm_min*(ErfTerm_min + TMath::Exp(mean[i]/tau)*(ExpFirstTerm_min*(2+Delta_gamma*tau)* ErfcFirstTerm_min + ExpSecondTerm_min*(-2+Delta_gamma*tau)* ErfcSecondTerm_min)); SinhTermsNorm[i] = MaxVal - minVal; } else { } } } for (UInt_t i(0); icalcLikelihoodInfo(abscissaErr); errTerm_ = errHist_->getLikelihood(); } else { errTerm_ = 1.0; } if ( effiFun_ ) { this->updateEffiSpline(effiPars_); effiTerm_ = effiFun_->evaluate(abscissa); if (effiTerm_>1.0){effiTerm_=1.0;} if (effiTerm_<0.0){effiTerm_=0.0;} } else { effiTerm_ = 1.0; } } void LauDecayTimePdf::calcLikelihoodInfoGen(Double_t abscissa, Double_t abscissaErr, Int_t curEvtTagFlv, Double_t curEvtMistag, Int_t curEvtTagCat, Int_t curEvtTrueTagFlv) { if (abscissa > this->maxAbscissa() || abscissa < this->minAbscissa()) { cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfoGen : Given value of the decay time: "<minAbscissa()<<","<maxAbscissa()<<"]."<Exit(EXIT_FAILURE); } if (abscissaErr > this->maxAbscissaError() || abscissaErr < this->minAbscissaError()) { cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfoGen : Given value of Delta t error: "<minAbscissaError()<<","<maxAbscissaError()<<"]."<Exit(EXIT_FAILURE); } // Initialise the various terms to zero if (type_ == Hist){ if ( pdfHist_ ) { pdfHist_->calcLikelihoodInfo(abscissa); pdfTerm_ = pdfHist_->getLikelihood(); } else { pdfTerm_ = 1.0; } }else{ // Reset the state to Good this->state(Good); // If we're not using the resolution function calculate the simple terms and return if (!this->doSmearing()) { this->calcNonSmearedTermsGen(abscissa, curEvtTagFlv, curEvtMistag, curEvtTagCat, curEvtTrueTagFlv); return; } //TODO how much to be added below for SimFitNormBd/SimFitNormBs/SimFitSigBd/SimFitSigBs // Get all the up to date parameter values std::vector frac(nGauss_); std::vector mean(nGauss_); std::vector sigma(nGauss_); Double_t tau(0.0); Double_t deltaM(0.0); Double_t fracPrompt(0.0); Double_t Delta_gamma(0.0); frac[0] = 1.0; for (UInt_t i(0); ivalue(); sigma[i] = sigma_[i]->value(); if (i != 0) { frac[i] = frac_[i-1]->value(); frac[0] -= frac[i]; } } if (type_ != Delta) { tau = tau_->value(); if (type_ == ExpTrig) { deltaM = deltaM_->value(); } if (type_ == DeltaExp) { fracPrompt = fracPrompt_->value(); } if (type_ == ExpHypTrig){ deltaM = deltaM_->value(); Delta_gamma = deltaGamma_->value(); } } // 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) { std::vector expTerms(nGauss_); std::vector cosTerms(nGauss_); std::vector sinTerms(nGauss_); std::vector coshTerms(nGauss_); std::vector sinhTerms(nGauss_); std::vector expTermsNorm(nGauss_); // TODO - TEL changed this name to make it compile - please check! std::vector SinhTermsNorm(nGauss_); // Calculate values of the PDF convoluated with each Gaussian for a given value of the abscsissa for (UInt_t i(0); icalcTrigExponent(deltaM, tau, x[i], sigma[i], exponentTermRe, exponentTermIm); // Elements related to the trigonometric function, i.e. convolution of Exp*Sin or Cos with Gauss Double_t sinTrigTermRe, sinTrigTermIm, cosTrigTermRe, cosTrigTermIm; this->calcTrigConv(deltaM, tau, x[i], sigma[i], sinTrigTermRe, sinTrigTermIm, kFALSE); this->calcTrigConv(deltaM, tau, x[i], sigma[i], cosTrigTermRe, cosTrigTermIm, kTRUE); // Combining elements of the full pdf LauComplex zExp(exponentTermRe, exponentTermIm); LauComplex zTrigSin(sinTrigTermRe, sinTrigTermIm); LauComplex zTrigCos(cosTrigTermRe, cosTrigTermIm); LauComplex sinConv = zExp * zTrigSin; LauComplex cosConv = zExp * zTrigCos; sinConv.scale(1.0/4.0); cosConv.scale(1.0/4.0); // Cosine*Exp and Sine*Exp terms cosTerms[i] = cosConv.re(); sinTerms[i] = sinConv.im(); // Normalisation Double_t umax = xMax - mean[i]; Double_t umin = xMin - mean[i]; expTermsNorm[i] = (1.0/2.0) * tau * (-1.0 + TMath::Erf(umax/(LauConstants::root2 * sigma[i])) + TMath::Erfc(umin/(LauConstants::root2 * sigma[i])) + TMath::Exp((pow(sigma[i], 2) - 2.0 * tau * (xMax + xMin - mean[i]))/(2.0 * pow(tau, 2))) * (TMath::Exp(xMax/tau) * TMath::Erfc((pow(sigma[i], 2) - xMin)/(LauConstants::root2 * tau))) + (TMath::Exp(xMin/tau) * TMath::Erfc((pow(sigma[i], 2) - xMax)/(LauConstants::root2 * tau)))); } else { } } // Typical case (2): B0s/B0sbar if (type_ == ExpHypTrig) { // LHCb convention if (method_ == DecayTime) { // Convolution of Exp*cosh (Exp*sinh) with a gaussian //Double_t OverallExpFactor = 0.25*TMath::Exp(-(x[i]-mean[i])*(x[i]-mean[i])/(2*sigma[i]*sigma[i])); //Double_t ExpFirstTerm = TMath::Exp((2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))*(2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))/(8*sigma[i]*sigma[i]*tau*tau)); //Double_t ExpSecondTerm = TMath::Exp((2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))*(2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))/(8*sigma[i]*sigma[i]*tau*tau)); //Double_t ErfFirstTerm = TMath::Erf((2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); //Double_t ErfSecondTerm = TMath::Erf((2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); //Double_t sinhConv = OverallExpFactor*(ExpFirstTerm*(1+ErfFirstTerm) + ExpSecondTerm*(-1+ErfSecondTerm)); //Double_t coshConv = OverallExpFactor*(ExpFirstTerm*(1+ErfFirstTerm) - ExpSecondTerm*(-1+ErfSecondTerm)); //cosTerms[i] = sinhConv; // sinTerms[i] = coshConv; //TODO: check this formula and try to simplify it! double OverallExpTerm_max = (1/(2*(-4 + Delta_gamma*Delta_gamma*tau*tau)))*tau*TMath::Exp(-0.5*Delta_gamma*(xMax + mean[i]) - xMax/tau); double ErfTerm_max = -2*Delta_gamma*tau*TMath::Exp(0.5*Delta_gamma*(xMax+mean[i])+xMax/tau)*TMath::Erf((xMax-mean[i])/(TMath::Sqrt(2)*sigma[i])); double ExpFirstTerm_max = TMath::Exp(xMax*Delta_gamma+(sigma[i]*sigma[i]*(-2 + Delta_gamma*tau)*(-2 + Delta_gamma*tau))/(8*tau*tau)); double ErfcFirstTerm_max = TMath::Erfc((2*(-xMax + mean[i])*tau + sigma[i]*sigma[i]*(2 - Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); double ExpSecondTerm_max = TMath::Exp(Delta_gamma*mean[i] + (sigma[i]*sigma[i]*(2 + Delta_gamma*tau)*(2 + Delta_gamma*tau))/(8*tau*tau)); double ErfcSecondTerm_max = TMath::Erfc((2*(-xMax + mean[i])*tau + sigma[i]*sigma[i]*(2 + Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); double MaxVal= OverallExpTerm_max*(ErfTerm_max + TMath::Exp(mean[i]/tau)*(ExpFirstTerm_max*(2+Delta_gamma*tau)* ErfcFirstTerm_max + ExpSecondTerm_max*(-2+Delta_gamma*tau)* ErfcSecondTerm_max)); double OverallExpTerm_min = (1/(2*(-4 + Delta_gamma*Delta_gamma*tau*tau)))*tau*TMath::Exp(-0.5*Delta_gamma*(xMin + mean[i]) - xMin/tau); double ErfTerm_min = -2*Delta_gamma*tau*TMath::Exp(0.5*Delta_gamma*(xMin+mean[i])+xMin/tau)*TMath::Erf((xMin-mean[i])/(TMath::Sqrt(2)*sigma[i])); double ExpFirstTerm_min = TMath::Exp(xMin*Delta_gamma+(sigma[i]*sigma[i]*(-2 + Delta_gamma*tau)*(-2 + Delta_gamma*tau))/(8*tau*tau)); double ErfcFirstTerm_min = TMath::Erfc((2*(-xMin + mean[i])*tau + sigma[i]*sigma[i]*(2 - Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); // TODO - TEL added this (currently identical to ExpSecondTerm_max) to get this to compile - please check!! double ExpSecondTerm_min = TMath::Exp(Delta_gamma*mean[i] + (sigma[i]*sigma[i]*(2 + Delta_gamma*tau)*(2 + Delta_gamma*tau))/(8*tau*tau)); double ErfcSecondTerm_min = TMath::Erfc((2*(-xMin + mean[i])*tau + sigma[i]*sigma[i]*(2 + Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau)); double minVal= OverallExpTerm_min*(ErfTerm_min + TMath::Exp(mean[i]/tau)*(ExpFirstTerm_min*(2+Delta_gamma*tau)* ErfcFirstTerm_min + ExpSecondTerm_min*(-2+Delta_gamma*tau)* ErfcSecondTerm_min)); SinhTermsNorm[i] = MaxVal - minVal; } else { } } } for (UInt_t i(0); icalcLikelihoodInfo(abscissaErr); errTerm_ = errHist_->getLikelihood(); } else { errTerm_ = 1.0; } if ( effiFun_ ) { this->updateEffiSpline(effiPars_); effiTerm_ = effiFun_->evaluate(abscissa); if (effiTerm_>1.0){effiTerm_=1.0;} if (effiTerm_<0.0){effiTerm_=0.0;} } else { effiTerm_ = 1.0; } } void LauDecayTimePdf::calcTrigExponent(Double_t deltaM, Double_t tau, Double_t x, Double_t sigma, Double_t& reTerm, Double_t& imTerm) { Double_t exponentTerm = TMath::Exp(-(2.0 * tau * x + pow(sigma, 2) * (pow(deltaM, 2) * pow(tau, 2) - 1.0))/(2.0 * pow(tau,2))); reTerm = exponentTerm * TMath::Cos(deltaM * (x - pow(sigma,2)/tau)); imTerm = - exponentTerm * TMath::Sin(deltaM * (x - pow(sigma,2)/tau)); } void LauDecayTimePdf::calcTrigConv(Double_t deltaM, Double_t tau, Double_t x, Double_t sigma, Double_t& reOutTerm, Double_t& imOutTerm, Bool_t trig) { Double_t reExpTerm, imExpTerm; LauComplex zExp; LauComplex zTrig1; LauComplex zTrig2; // Calculation for the sine or cosine term if (!trig) { reExpTerm = TMath::Sin(2.0 * deltaM * (x + pow(sigma,2)/tau)); imExpTerm = 2.0 * TMath::Sin(pow(deltaM * (x + pow(sigma,2)/tau), 2)); } else { reExpTerm = TMath::Cos(2.0 * deltaM * (x + pow(sigma,2)/tau)); imExpTerm = TMath::Sin(2.0 * deltaM * (x + pow(sigma,2)/tau)); } // Exponential term in front of Erfc/Erfi terms zExp.setRealPart(reExpTerm); zExp.setImagPart(imExpTerm); // Nominal Erfc term (common to both sine and cosine expressions zTrig1.setRealPart(-(tau * x - pow(sigma,2))/(LauConstants::root2 * tau * sigma)); zTrig1.setImagPart(-(deltaM * sigma)/ LauConstants::root2); // Second term for sine (Erfi) or cosine (Erfc) - notice the re-im swap and sign change zTrig2.setRealPart(-zTrig1.im()); zTrig2.setImagPart(-zTrig1.re()); // Calculation of Erfc and Erfi (if necessary) LauComplex term1 = ComplexErfc(zTrig1.re(), zTrig1.im()); LauComplex term2; if (!trig) { term2 = Erfi(zTrig2.re(), zTrig2.im()); } else { term2 = ComplexErfc(zTrig2.re(), zTrig2.im()); } // Multiplying all elemnets of the convolution LauComplex output = zExp * term1 + term2; reOutTerm = output.re(); imOutTerm = output.im(); } LauComplex LauDecayTimePdf::ComplexErf(Double_t x, Double_t y) { // Evaluate Erf(x + iy) using an infinite series approximation // From Abramowitz & Stegun (http://people.math.sfu.ca/~cbm/aands/page_299.htm) if (x==0){ // cout << "WARNING: Set x value to 1e-100 to avoid division by 0." << endl; x = 1e-100; } int n = 20; // this cotrols the number of iterations of the sum LauComplex ErfTerm(TMath::Erf(x),0.); LauComplex CosSineTerm(1-cos(2*x*y), sin(2*x*y)); CosSineTerm.rescale(TMath::Exp(-x*x)/(2*TMath::Pi()*x)); LauComplex firstPart = ErfTerm + CosSineTerm; LauComplex SumTerm(0,0); for (int k = 1; k<=n; k++){ Double_t f_k = 2*x*(1 - cos(2*x*y)*cosh(k*y)) + k*sin(2*x*y)*sinh(k*y); Double_t g_k = 2*x*sin(2*x*y)*cosh(k*y) + k*cos(2*x*y)*sinh(k*y); LauComplex fgTerm(f_k, g_k); fgTerm.rescale(TMath::Exp(-0.25*k*k)/(k*k + 4*x*x)); SumTerm += fgTerm; } SumTerm.rescale((2/TMath::Pi())*TMath::Exp(-x*x)); LauComplex result = firstPart + SumTerm; return result; } LauComplex LauDecayTimePdf::Erfi(Double_t x, Double_t y) { // Erfi(z) = -I*Erf(I*z) where z = x + iy double x_prime = -y; double y_prime = x; LauComplex a = ComplexErf(x_prime, y_prime); LauComplex result(a.im(), -a.re()); return result; } LauComplex LauDecayTimePdf::ComplexErfc(Double_t x, Double_t y) { // Erfc(z) = 1 - Erf(z) (z = x + iy) LauComplex one(1., 0.); LauComplex result = one - ComplexErf(x,y); return result; } void LauDecayTimePdf::calcNonSmearedTerms(Double_t abscissa, UInt_t iEvt) { if (type_ == Hist ){ cerr << "It is a histogrammed PDF" << endl; return; } if (type_ == Delta) { return; } Double_t tau = tau_->value(); + Double_t deltaM = deltaM_->value(); // Calculate the terms related to cosine and sine not normalised if (type_ == ExpTrig) { - Double_t deltaM = deltaM_->value(); if (method_ == DecayTime) { expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); } if (method_ == DecayTimeDiff) { expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); } cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_; coshTerm_ = expTerm_; sinhTerm_ = 0.0; } // Calculate the terms related to cosine not normalised - if (type_ == SimFitNormBd) { + if (type_ == SimFitNormBd || type_ == SimFitNormBs) { - Double_t deltaM = deltaM_->value(); if (method_ == DecayTime) { expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); } if (method_ == DecayTimeDiff) { expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); } - //asymTerm_ = (1-(r*Acp))*(1-(r*Af)); asymTerm_ = ARaw_->unblindValue(); Int_t curEvtTagFlv = flavTag_->getEvtTagFlvVals(iEvt); Int_t curEvtTrueTagFlv = flavTag_->getEvtTrueTagVals(iEvt); Double_t omega(0); Double_t omegabar(0); if (curEvtTagFlv==1){ omegabar = flavTag_->getOmega(iEvt); } else if (curEvtTagFlv==-1){ omega = flavTag_->getOmega(iEvt); } Double_t Ap = AProd_->unblindValue(); omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega); omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega); cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*asymTerm_*omegaTerm_*curEvtTrueTagFlv; sinTerm_ = 0.0; coshTerm_ = expTerm_*asymTerm_*omegaTermh_; sinhTerm_ = 0.0; + if (type_ == SimFitNormBs){ + Double_t deltaGamma = deltaGamma_->value(); + coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0); + } + } // Calculate the terms related to cosine and sine not normalised - if (type_ == SimFitSigBd) { + if (type_ == SimFitSigBd || type_ == SimFitSigBs) { - Double_t deltaM = deltaM_->value(); if (method_ == DecayTime) { expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); } if (method_ == DecayTimeDiff) { expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); } Int_t curEvtTagFlv = flavTag_->getEvtTagFlvVals(iEvt); Double_t omega(0); Double_t omegabar(0); if (curEvtTagFlv==1){ omegabar = flavTag_->getOmega(iEvt); } else if (curEvtTagFlv==-1){ omega = flavTag_->getOmega(iEvt); } Double_t Ap = AProd_->unblindValue(); omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega); omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega); cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*omegaTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_*omegaTerm_; coshTerm_ = expTerm_*omegaTermh_; sinhTerm_ = 0.0; - - } - - // Calculate the terms related to cosine and cosh not normalised - if (type_ == SimFitNormBs) { - - Double_t deltaM = deltaM_->value(); - Double_t deltaGamma = deltaGamma_->value(); - if (method_ == DecayTime) { - expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); - } - if (method_ == DecayTimeDiff) { - expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); - } - //asymTerm_ = (1-(r*Acp))*(1-(r*Af)); - asymTerm_ = ARaw_->unblindValue(); - - Int_t curEvtTagFlv = flavTag_->getEvtTagFlvVals(iEvt); - Int_t curEvtTrueTagFlv = flavTag_->getEvtTrueTagVals(iEvt); - - Double_t omega(0); - Double_t omegabar(0); - if (curEvtTagFlv==1){ - omegabar = flavTag_->getOmega(iEvt); - } else if (curEvtTagFlv==-1){ - omega = flavTag_->getOmega(iEvt); - } - Double_t Ap = AProd_->unblindValue(); - omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega); - omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega); - - cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*asymTerm_*omegaTerm_*curEvtTrueTagFlv; - sinTerm_ = 0.0; - coshTerm_ = TMath::CosH(deltaGamma*abscissa/2.0)*expTerm_*asymTerm_*omegaTermh_; - sinhTerm_ = 0.0; - - } - - // Calculate the terms related to cosine and sine not normalised - if (type_ == SimFitSigBs) { - - Double_t deltaM = deltaM_->value(); - Double_t deltaGamma = deltaGamma_->value(); - if (method_ == DecayTime) { - expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); - } - if (method_ == DecayTimeDiff) { - expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); - } - Int_t curEvtTagFlv = flavTag_->getEvtTagFlvVals(iEvt); - - Double_t omega(0); - Double_t omegabar(0); - if (curEvtTagFlv==1){ - omegabar = flavTag_->getOmega(iEvt); - } else if (curEvtTagFlv==-1){ - omega = flavTag_->getOmega(iEvt); + if (type_ == SimFitNormBs){ + Double_t deltaGamma = deltaGamma_->value(); + coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0); + sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_*omegaTermh_; } - - Double_t Ap = AProd_->unblindValue(); - omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega); - omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega); - - cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*omegaTerm_; - sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_*omegaTerm_; - coshTerm_ = TMath::CosH(deltaGamma*abscissa/2.0)*expTerm_*omegaTermh_; - sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_*omegaTermh_; } - + // Calculate the terms related to cosine, sine, cosh and sinh not normalised (no decayTimeDiff implemented) if (type_ == ExpHypTrig) { - Double_t deltaM = deltaM_->value(); Double_t deltaGamma = deltaGamma_->value(); expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_; coshTerm_ = TMath::CosH(deltaGamma*abscissa/2.0)*expTerm_; sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_; } if ( effiFun_ ) { this->updateEffiSpline(effiPars_); effiTerm_ = effiFun_->evaluate(abscissa); if (effiTerm_>1.0){effiTerm_=1.0;} if (effiTerm_<0.0){effiTerm_=0.0;} } else { effiTerm_ = 1.0; } } void LauDecayTimePdf::calcNonSmearedTermsGen(Double_t abscissa, Int_t curEvtTagFlv, Double_t curEvtMistag, Int_t curEvtTagCat, Int_t curEvtTrueTagFlv) { if (type_ == Hist ){ cerr << "It is a histogrammed PDF" << endl; return; } if (type_ == Delta) { return; } Double_t tau = tau_->value(); // Calculate the terms related to cosine and sine not normalised if (type_ == ExpTrig) { Double_t deltaM = deltaM_->value(); if (method_ == DecayTime) { expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); } if (method_ == DecayTimeDiff) { expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); } cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_; coshTerm_ = expTerm_; sinhTerm_ = 0.0; } // Calculate the terms related to cosine not normalised - if (type_ == SimFitNormBd) { + if (type_ == SimFitNormBd || type_ == SimFitNormBs) { Double_t deltaM = deltaM_->value(); if (method_ == DecayTime) { expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); } if (method_ == DecayTimeDiff) { expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); } - //asymTerm_ = (1-(r*Acp))*(1-(r*Af)); asymTerm_ = ARaw_->unblindValue(); Double_t omega(0); Double_t omegabar(0); if (curEvtTagFlv==1){ omegabar = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat); } else if (curEvtTagFlv==-1){ omega = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat); } Double_t Ap = AProd_->unblindValue(); omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega); omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega); cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*asymTerm_*omegaTerm_*curEvtTrueTagFlv; sinTerm_ = 0.0; coshTerm_ = expTerm_*asymTerm_*omegaTermh_; sinhTerm_ = 0.0; + if (type_ == SimFitNormBs){ + Double_t deltaGamma = deltaGamma_->value(); + coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0); + } + } // Calculate the terms related to cosine and sine not normalised - if (type_ == SimFitSigBd) { + if (type_ == SimFitSigBd || type_ == SimFitSigBs) { Double_t deltaM = deltaM_->value(); if (method_ == DecayTime) { expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); } if (method_ == DecayTimeDiff) { expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); } Double_t omega(0); Double_t omegabar(0); if (curEvtTagFlv==1){ omegabar = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat); } else if (curEvtTagFlv==-1){ omega = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat); } Double_t Ap = AProd_->unblindValue(); omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega); omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega); cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*omegaTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_*omegaTerm_; coshTerm_ = expTerm_*omegaTermh_; sinhTerm_ = 0.0; - } - - // Calculate the terms related to cosine and cosh not normalised - if (type_ == SimFitNormBs) { - - Double_t deltaM = deltaM_->value(); - Double_t deltaGamma = deltaGamma_->value(); - if (method_ == DecayTime) { - expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); - } - if (method_ == DecayTimeDiff) { - expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); - } - //asymTerm_ = (1-(r*Acp))*(1-(r*Af)); - asymTerm_ = ARaw_->unblindValue(); - - Double_t omega(0); - Double_t omegabar(0); - if (curEvtTagFlv==1){ - omegabar = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat); - } else if (curEvtTagFlv==-1){ - omega = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat); - } - Double_t Ap = AProd_->unblindValue(); - omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega); - omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega); - - cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*asymTerm_*omegaTerm_*curEvtTrueTagFlv; - sinTerm_ = 0.0; - coshTerm_ = TMath::CosH(deltaGamma*abscissa/2.0)*expTerm_*asymTerm_*omegaTermh_; - sinhTerm_ = 0.0; - - } - - // Calculate the terms related to cosine and sine not normalised - if (type_ == SimFitSigBs) { - - Double_t deltaM = deltaM_->value(); - Double_t deltaGamma = deltaGamma_->value(); - if (method_ == DecayTime) { - expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); - } - if (method_ == DecayTimeDiff) { - expTerm_ = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau); - } - Double_t omega(0); - Double_t omegabar(0); - if (curEvtTagFlv==1){ - omegabar = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat); - } else if (curEvtTagFlv==-1){ - omega = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat); + if (type_ == SimFitNormBs){ + Double_t deltaGamma = deltaGamma_->value(); + coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0); + sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_*omegaTermh_; } - - Double_t Ap = AProd_->unblindValue(); - omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega); - omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega); - - cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*omegaTerm_; - sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_*omegaTerm_; - coshTerm_ = TMath::CosH(deltaGamma*abscissa/2.0)*expTerm_*omegaTermh_; - sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_*omegaTermh_; - } - + // Calculate the terms related to cosine, sine, cosh and sinh not normalised (no decayTimeDiff implemented) if (type_ == ExpHypTrig) { Double_t deltaM = deltaM_->value(); Double_t deltaGamma = deltaGamma_->value(); expTerm_ = TMath::Exp(-abscissa/tau)/(2.0*tau); cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_; coshTerm_ = TMath::CosH(deltaGamma*abscissa/2.0)*expTerm_; sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_; } if ( effiFun_ ) { this->updateEffiSpline(effiPars_); effiTerm_ = effiFun_->evaluate(abscissa); if (effiTerm_>1.0){effiTerm_=1.0;} if (effiTerm_<0.0){effiTerm_=0.0;} } else { effiTerm_ = 1.0; } } Double_t LauDecayTimePdf::normExpHypTerm(Double_t Abs) { Double_t tau = tau_->value(); Double_t deltaGamma = deltaGamma_->value(); Double_t y = tau*deltaGamma/2; Double_t nonTrigTerm = -(TMath::Exp(-Abs/tau))/(1 - y*y); Double_t cosHTerm = TMath::CosH(deltaGamma*Abs/2); Double_t sinHTerm = TMath::SinH(deltaGamma*Abs/2); Double_t normTerm = nonTrigTerm*(cosHTerm + y*sinHTerm); return normTerm; } Double_t LauDecayTimePdf::normExpHypTermDep(Double_t Abs) { Double_t tau = tau_->value(); Double_t deltaGamma = deltaGamma_->value(); Double_t y = tau*deltaGamma/2; Double_t nonTrigTerm = -(TMath::Exp(-Abs/tau))/(1 - y*y); Double_t cosHTerm = TMath::CosH(deltaGamma*Abs/2); Double_t sinHTerm = TMath::SinH(deltaGamma*Abs/2); Double_t normTerm = nonTrigTerm*(sinHTerm + y*cosHTerm); return normTerm; } void LauDecayTimePdf::storeIntegralTimeDep(Double_t minAbs, Double_t maxAbs) { Double_t tau = tau_->value(); // Normalisation factor for B0 decays if (type_ == ExpTrig || type_ == SimFitNormBd || type_ == SimFitSigBd ) { if (method_ == DecayTime) { normTermExp_ = tau * (TMath::Exp(-minAbs/tau) - TMath::Exp(-maxAbs/tau)); } if (method_ == DecayTimeDiff) { normTermExp_ = tau * (2.0 - TMath::Exp(-maxAbs/tau) - TMath::Exp(-minAbs/tau)); } } // Normalisation factor for Bs decays if (type_ == ExpHypTrig) { normTermCosh_ = normExpHypTerm(maxAbs) - normExpHypTerm(minAbs); normTermSinh_ = normExpHypTermDep( maxAbs) - normExpHypTermDep( minAbs); } } Double_t LauDecayTimePdf::generateError(Bool_t forceNew) { if (errHist_ && (forceNew || !abscissaErrorGenerated_)) { LauFitData errData = errHist_->generate(0); abscissaError_ = errData.find(this->varErrName())->second; 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()) { 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_ != 0 ) { 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_ != 0 ) { cerr<<"WARNING in LauDecayTimePdf::setHistoPdf : PDF histogram already set, not doing it again."<varName(), hist, this->minAbscissa(), this->maxAbscissa()); } void LauDecayTimePdf::setEffiSpline(Lau1DCubicSpline* spline) { if ( effiFun_ != 0 ) { cerr<<"WARNING in LauDecayTimePdf::setEffiPdf : efficiency function already set, not doing it again."< effiPars) { if ( effiFun_ != 0 ) { cerr<<"WARNING in LauDecayTimePdf::setEffiPdf : efficiency function already set, not doing it again."<getnKnots()){ cerr<<"ERROR in LauDecayTimePdf::setEffiPdf : number of efficiency parameters is not equal to the number of spline knots."<Exit(EXIT_FAILURE); } effiPars_ = effiPars; } 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 params) { if (!params.empty()){ std::vector yvals; yvals.reserve(params.size()); for (ULong_t i(0); iunblindValue()); } effiFun_->updateYValues(yvals); } } void LauDecayTimePdf::setAsymmetries(LauParameter* AProd, LauParameter* ARaw) { AProd_ = AProd; ARaw_ = ARaw; } void LauDecayTimePdf::setFlavourTagging(LauFlavTag* flavTag) { flavTag_ = flavTag; }