diff --git a/inc/LauTimeDepFitModel.hh b/inc/LauTimeDepFitModel.hh index 826e3ba..9629293 100644 --- a/inc/LauTimeDepFitModel.hh +++ b/inc/LauTimeDepFitModel.hh @@ -1,673 +1,706 @@ /* 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 LauTimeDepFitModel.hh \brief File containing declaration of LauTimeDepFitModel class. */ /*! \class LauTimeDepFitModel \brief Class for defining a time-dependent fit model. LauTimeDepFitModel is a class that allows the user to define a three-body Dalitz plot according to the isobar model, i.e. defining a set of resonances that have complex amplitudes that can interfere with each other. It extends the LauSimpleFitModel and LauCPFitModel models in that it allows the fitting of time-dependent particle/antiparticle decays to flavour-conjugate Dalitz plots, including their interference through mixing. */ #ifndef LAU_TIMEDEP_FIT_MODEL #define LAU_TIMEDEP_FIT_MODEL #include #include #include #include #include "TString.h" #include "TStopwatch.h" #include "TSystem.h" #include "LauAbsFitModel.hh" #include "LauConstants.hh" #include "LauEmbeddedData.hh" #include "LauParameter.hh" #include "LauFlavTag.hh" class LauAbsBkgndDPModel; class LauAbsCoeffSet; class LauAbsPdf; class LauDecayTimePdf; class LauIsobarDynamics; class LauKinematics; class LauScfMap; class LauTimeDepFitModel : public LauAbsFitModel { public: //! Possible CP eigenvalues (the intrinsic CP of the final state particles) enum CPEigenvalue { CPOdd = -1, /*!< CP odd final state */ CPEven = 1 /*!< CP even final state */ }; //! Constructor /*! \param [in] modelB0bar DP model for the antiparticle \param [in] modelB0 DP model for the particle \param [in] flavTag flavour tagging information */ LauTimeDepFitModel(LauIsobarDynamics* modelB0bar, LauIsobarDynamics* modelB0, LauFlavTag* flavTag); //! Destructor virtual ~LauTimeDepFitModel(); //! Set the signal event yield /*! \param [in] nSigEvents contains the signal yield and option to fix it */ virtual void setNSigEvents(LauParameter* nSigEvents); //! Set the signal event yield and asymmetry /*! \param [in] nSigEvents contains the signal yield and option to fix it \param [in] sigAsym contains the signal asymmetry and option to fix it */ virtual void setNSigEvents(LauParameter* nSigEvents, LauParameter* sigAsym); //! Set the number of background events /*! The name of the parameter must be that of the corresponding background category (so that it can be correctly assigned) \param [in] nBkgndEvents contains the name, yield and option to fix the yield of the background */ virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents); //! Set the background event yield and asymmetry /*! \param [in] nBkgEvents contains the background yield and option to fix it \param [in] BkgAsym contains the background asymmetry and option to fix it */ virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym); + //! Set the background DP models + /*! + \param [in] bkgndClass the name of the background class + \param [in] model the DP model of the background + */ + void setBkgndDPModels(const TString& bkgndClass, LauAbsBkgndDPModel* model); + //! Switch on/off storage of amplitude info in generated ntuple void storeGenAmpInfo(Bool_t storeInfo) { storeGenAmpInfo_ = storeInfo; } //! Set CP eigenvalue /*! The CP eigenvalue can be supplied on an event-by-event basis, e.g. if the data contains daughters that are D0 mesons that can decay to either K+ K- (CP even) or KS pi0 (CP odd). This method allows you to set the default value that should be used if the data does not contain this information as well as the name of the variable in the data that will specify this information. If completely unspecified all events will be assumed to be CP even. \param defaultCPEV the default for the eigenvalue \param evVarName the variable name in the data tree that specifies the CP eigenvalue */ inline void setCPEigenvalue( const CPEigenvalue defaultCPEV, const TString& cpevVarName = "" ) { cpEigenValue_ = defaultCPEV; cpevVarName_ = cpevVarName; } //! Set the DP amplitude coefficients /*! \param [in] coeffSet the set of coefficients */ void setAmpCoeffSet(LauAbsCoeffSet* coeffSet); //! Set the decay time PDFs /*! \param [in] tagCat the tagging category for which the PDF should be used \param [in] pdf the signal decay time PDF */ void setSignalDtPdf(const Int_t tagCat, LauDecayTimePdf* pdf); + //! Set the decay time PDFs + /*! + \param [in] tagCat the tagging category for which the PDF should be used + \param [in] pdf the background decay time PDF + */ + void setBackgroundDtPdf(const Int_t tagCat, LauDecayTimePdf* pdf); + //! Set the signal PDF for a given variable /*! \param [in] tagCat the tagging category for which the PDF should be used \param [in] pdf the PDF to be added to the signal model */ void setSignalPdfs(const Int_t tagCat, LauAbsPdf* pdf); + //! Set the background PDF + /*! + \param [in] bkgndClass the name of the background class + \param [in] pdf the PDF to be added to the background model + */ + void setBkgndPdf(const TString& bkgndClass, LauAbsPdf* pdf); + void setSignalFlavTagPdfs( const Int_t tagCat, LauAbsPdf* pdf); void setBkgdFlavTagPdfs( const TString name, LauAbsPdf* pdf); //! Embed full simulation events for the signal, rather than generating toy from the PDFs /*! \param [in] tagCat the tagging category for which the file should be used \param [in] fileName the name of the file containing the events \param [in] treeName the name of the tree \param [in] reuseEventsWithinEnsemble \param [in] reuseEventsWithinExperiment \param [in] useReweighting */ void embedSignal(const Int_t tagCat, const TString& fileName, const TString& treeName, const Bool_t reuseEventsWithinEnsemble, const Bool_t reuseEventsWithinExperiment = kFALSE); + void embedBkgnd(Int_t tagCat, const TString& bkgndClass, const TString& fileName, const TString& treeName, + Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment = kFALSE); + //! Set the value of the mixing phase /*! \param [in] phiMix the value of the mixing phase \param [in] fixPhiMix whether the value should be fixed or floated \param [in] useSinCos whether to use the sine and cosine as separate parameters or to just use the mixing phase itself */ void setPhiMix(const Double_t phiMix, const Bool_t fixPhiMix, const Bool_t useSinCos = kFALSE); //! Initialise the fit virtual void initialise(); //! Initialise the signal DP model virtual void initialiseDPModels(); //! Recalculate Normalization the signal DP models virtual void recalculateNormalisation(); //! Update the coefficients virtual void updateCoeffs(); // Toy MC generation and fitting overloaded functions virtual Bool_t genExpt(); //! Set the maximum value of A squared to be used in the accept/reject /*! \param [in] value the new value */ inline void setASqMaxValue(const Double_t value) {aSqMaxSet_ = value;} //! Weight events based on the DP model /*! \param [in] dataFileName the name of the data file \param [in] dataTreeName the name of the data tree */ virtual void weightEvents( const TString& dataFileName, const TString& dataTreeName ); //! Calculate things that depend on the fit parameters after they have been updated by Minuit virtual void propagateParUpdates(); //! Read in the input fit data variables, e.g. m13Sq and m23Sq virtual void cacheInputFitVars(); //! Check the initial fit parameters virtual void checkInitFitParams(); //! Get the fit results and store them /*! \param [in] tablePrefixName prefix for the name of the output file */ virtual void finaliseFitResults(const TString& tablePrefixName); //! Save the pdf Plots for all the resonances of experiment number fitExp /*! TODO - not working in this model!! \param [in] label prefix for the file name to be saved */ virtual void savePDFPlots(const TString& label); //! Save the pdf Plots for the sum of ressonances correspondint to "sin" of experiment number fitExp /*! TODO - not working in this model!! \param [in] label prefix for the file name to be saved \param [in] spin spin of the wave to be saved */ virtual void savePDFPlotsWave(const TString& label, const Int_t& spin); //! Print the fit fractions, total DP rate and mean efficiency /*! \param [out] output the stream to which to print */ virtual void printFitFractions(std::ostream& output); //! Print the asymmetries /*! \param [out] output the stream to which to print */ virtual void printAsymmetries(std::ostream& output); //! Write the fit results in latex table format /*! \param [in] outputFile the name of the output file */ virtual void writeOutTable(const TString& outputFile); //! Store the per event likelihood values virtual void storePerEvtLlhds(); // Methods to do with calculating the likelihood functions // and manipulating the fitting parameters. //! Get the total likelihood for each event /*! \param [in] iEvt the event number */ virtual Double_t getTotEvtLikelihood(const UInt_t iEvt); //! Calculate the signal and background likelihoods for the DP for a given event /*! \param [in] iEvt the event number */ virtual void getEvtDPDtLikelihood(const UInt_t iEvt); //! Determine the signal and background likelihood for the extra variables for a given event /*! \param [in] iEvt the event number */ virtual void getEvtExtraLikelihoods(const UInt_t iEvt); virtual void getEvtFlavTagLikelihood(const UInt_t iEvt); //! Get the total number of events /*! \return the total number of events */ virtual Double_t getEventSum() const; //! Set the fit parameters for the DP model void setSignalDPParameters(); //! Set the fit parameters for the decay time PDFs void setDecayTimeParameters(); //! Set the fit parameters for the extra PDFs void setExtraPdfParameters(); //! Set the initial yields void setFitNEvents(); //! Set the calibration parameters void setCalibParams(); //! Set the efficiency parameters void setEffiParams(); //! Set the asymmetry parameters void setAsymParams(); //! Set the tagging asymmetry parameters void setFlavTagAsymParams(); //! Set-up other parameters that are derived from the fit results, e.g. fit fractions void setExtraNtupleVars(); //! Set production and detections asymmetries void setAsymmetries(const Double_t AProd, const Bool_t AProdFix, const Double_t ARaw, const Bool_t ARawFix); //! Randomise the initial fit parameters void randomiseInitFitPars(); //! Method to set up the storage for background-related quantities called by setBkgndClassNames virtual void setupBkgndVectors(); //! Calculate the CP asymmetries /*! \param [in] initValues is this before or after the fit */ void calcAsymmetries(const Bool_t initValues = kFALSE); //! Finalise value of mixing phase void checkMixingPhase(); //! Return the map of signal decay time PDFs typedef std::map< Int_t, LauDecayTimePdf*> LauTagCatDtPdfMap; LauTagCatDtPdfMap& getSignalDecayTimePdfs(){return signalDecayTimePdfs_;} + //! Return the map of background decay time PDFs + LauTagCatDtPdfMap& getBackgroundDecayTimePdfs(){return backgroundDecayTimePdfs_;} + protected: typedef std::map< Int_t, LauParameter> LauTagCatParamMap; typedef std::map< Int_t, LauPdfList > LauTagCatPdfListMap; typedef std::map< Int_t, LauAbsPdf* > LauTagCatPdfMap; typedef std::map< TString, LauAbsPdf* > LauBkgdPdfMap; typedef std::map< Int_t, Int_t > LauTagCatEventsMap; typedef std::map< Int_t, LauEmbeddedData* > LauTagCatEmbDataMap; typedef std::map< Int_t, std::pair > LauTagCatGenInfo; typedef std::map< std::pair, LauTagCatGenInfo > LauGenInfo; typedef std::vector LauTagCatEmbDataMapList; typedef std::vector LauBkgndDPModelList; typedef std::vector LauBkgndPdfsList; typedef std::vector LauBkgndYieldList; typedef std::vector LauBkgndReuseEventsList; //! Determine the number of events to generate for each hypothesis LauGenInfo eventsToGenerate(); //! Generate signal event Bool_t generateSignalEvent(); + //! Generate background event + /*! + \param [in] bgID ID number of the background class + */ + Bool_t generateBkgndEvent(UInt_t bgID); + //! Setup the required ntuple branches void setupGenNtupleBranches(); //! Store all of the DP and decay time information void setDPDtBranchValues(); //! Generate from the extra PDFs /*! \param [in] extraPdfs the list of extra PDFs \param [in] embeddedData the embedded data sample */ void generateExtraPdfValues(LauPdfList* extraPdfs, LauEmbeddedData* embeddedData); //! Add sPlot branches for the extra PDFs /*! \param [in] extraPdfs the list of extra PDFs \param [in] prefix the list of prefixes for the branch names */ void addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix); //! Set the branches for the sPlot ntuple with extra PDFs /*! \param [in] extraPdfs the list of extra PDFs \param [in] prefix the list of prefixes for the branch names \param [in] iEvt the event number */ Double_t setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt); //! Update the signal events after Minuit sets background parameters void updateSigEvents(); //! Add branches to store experiment number and the event number within the experiment virtual void setupSPlotNtupleBranches(); //! Returns the names of all variables in the fit virtual LauSPlot::NameSet variableNames() const; //! Returns the names and yields of species that are free in the fit virtual LauSPlot::NumbMap freeSpeciesNames() const; //! Returns the names and yields of species that are fixed in the fit virtual LauSPlot::NumbMap fixdSpeciesNames() const; //! Returns the species and variables for all 2D PDFs in the fit virtual LauSPlot::TwoDMap twodimPDFs() const; //! Check if the signal is split into well-reconstructed and mis-reconstructed types virtual Bool_t splitSignal() const { return kFALSE; } //! Check if the mis-reconstructed signal is to be smeared in the DP virtual Bool_t scfDPSmear() const { return kFALSE; } //! Add the parameters from each PDF into the fit /*! \param [in] theMap the container of PDFs */ UInt_t addParametersToFitList(LauTagCatPdfListMap& theMap); //! Add the parameters from each decay time PDF into the fit /*! \param [in] theMap the container of PDFs */ UInt_t addParametersToFitList(LauTagCatDtPdfMap& theMap); //! Calculate the component integrals of the interference term void calcInterferenceTermIntegrals(); //! Calculate the total integral of the interference term void calcInterTermNorm(); private: //! Dalitz plot PDF for the antiparticle LauIsobarDynamics* sigModelB0bar_; //! Dalitz plot PDF for the particle LauIsobarDynamics* sigModelB0_; //! Kinematics object for antiparticle LauKinematics* kinematicsB0bar_; //! Kinematics object for particle LauKinematics* kinematicsB0_; - //! The B0 background Dalitz plot models - LauBkgndDPModelList BkgndDPModelsB0_; + //! The background Dalitz plot models + LauBkgndDPModelList BkgndDPModels_; - //! The B0-bar background Dalitz plot models - LauBkgndDPModelList BkgndDPModelsB0bar_; + //! The background PDFs + LauBkgndPdfsList BkgndPdfs_; - //! The B0 background PDFs - LauBkgndPdfsList BkgndPdfsB0_; - - //! The B0bar background PDFs - LauBkgndPdfsList BkgndPdfsB0bar_; + //! Background boolean + Bool_t usingBkgnd_; //! LauFlavTag object for flavour tagging LauFlavTag* flavTag_; //! Flavour tag for current event Int_t curEvtTagFlv_; //! Tagging category for current event Int_t curEvtTagCat_; //! Per event mistag for current event Double_t curEvtMistag_; //! True tag for the current event Int_t curEvtTrueTagFlv_; //! Number of signal components UInt_t nSigComp_; //! Number of signal DP parameters UInt_t nSigDPPar_; //! Number of decay time PDF parameters UInt_t nDecayTimePar_; //! Number of extra PDF parameters UInt_t nExtraPdfPar_; //! Number of normalisation parameters (yields, asymmetries) UInt_t nNormPar_; //! Number of calibration parameters (p0, p1) UInt_t nCalibPar_; //! Number of efficiency parameters (p0, p1) UInt_t nEffiPar_; //! Number of asymmetry parameters UInt_t nAsymPar_; //! Number of tagging asymmetry parameters UInt_t nTagAsymPar_; //! The complex coefficients for antiparticle std::vector coeffsB0bar_; //! The complex coefficients for particle std::vector coeffsB0_; //! Magnitudes and Phases or Real and Imaginary Parts std::vector coeffPars_; //! The integrals of the efficiency corrected amplitude cross terms for each pair of amplitude components /*! Calculated as the sum of A* x Abar x efficiency */ std::vector< std::vector > fifjEffSum_; //! The normalisation for the term 2.0*Re(A*Abar*phiMix) Double_t interTermReNorm_; //! The normalisation for the term 2.0*Im(A*Abar*phiMix) Double_t interTermImNorm_; //! The antiparticle fit fractions LauParArray fitFracB0bar_; //! The particle fit fractions LauParArray fitFracB0_; //! The fit fraction asymmetries std::vector fitFracAsymm_; //! A_CP parameter std::vector acp_; //! The mean efficiency for the antiparticle LauParameter meanEffB0bar_; //! The mean efficiency for the particle LauParameter meanEffB0_; //! The average DP rate for the antiparticle LauParameter DPRateB0bar_; //! The average DP rate for the particle LauParameter DPRateB0_; //! Signal yields LauParameter* signalEvents_; //! Signal asymmetry LauParameter* signalAsym_; //! Background yield(s) LauBkgndYieldList bkgndEvents_; //! Background asymmetries(s) LauBkgndYieldList bkgndAsym_; //! CP eigenvalue variable name TString cpevVarName_; //! CP eigenvalue for current event CPEigenvalue cpEigenValue_; //! Vector to store event CP eigenvalues std::vector evtCPEigenVals_; //! The mass difference between the neutral mass eigenstates LauParameter deltaM_; //! The width difference between the neutral mass eigenstates LauParameter deltaGamma_; //! The lifetime LauParameter tau_; //! The mixing phase LauParameter phiMix_; //! The sine of the mixing phase LauParameter sinPhiMix_; //! The cosine of the mixing phase LauParameter cosPhiMix_; //! Flag whether to use the sine and cosine of the mixing phase or the phase itself as the fit parameters Bool_t useSinCos_; //! e^{-i*phiMix} LauComplex phiMixComplex_; - //! Decay time PDFs (one per tagging category) + //! Signal Decay time PDFs (one per tagging category) LauTagCatDtPdfMap signalDecayTimePdfs_; + //! Background Decay time PDFs (one per tagging category) + LauTagCatDtPdfMap backgroundDecayTimePdfs_; + //! Decay time for the current event Double_t curEvtDecayTime_; //! Decay time error for the current event Double_t curEvtDecayTimeErr_; //! PDFs for other variables LauTagCatPdfListMap sigExtraPdf_; //! eta PDFs for each TagCat LauTagCatPdfMap sigFlavTagPdf_; //! eta PDFs for each background LauBkgdPdfMap bkgdFlavTagPdf_; //! Production asymmetry between B0 and B0bar LauParameter AProd_; //! Detection asymmetry between B0 and B0bar LauParameter ARaw_; // Toy generation stuff //! The maximum allowed number of attempts when generating an event Int_t iterationsMax_; //! The number of unsucessful attempts to generate an event so far Int_t nGenLoop_; //! The value of A squared for the current event Double_t ASq_; //! The maximum value of A squared that has been seen so far while generating Double_t aSqMaxVar_; //! The maximum allowed value of A squared Double_t aSqMaxSet_; //! Flag for storage of amplitude info in generated ntuple Bool_t storeGenAmpInfo_; //! The signal event tree for embedding fully simulated events LauTagCatEmbDataMap signalTree_; //! The background event tree for embedding fully simulated events LauTagCatEmbDataMapList bkgndTree_; //! Boolean to control reuse of embedded signal events Bool_t reuseSignal_; //! Vector of booleans to reuse background events LauBkgndReuseEventsList reuseBkgnd_; // Likelihood values //! Signal DP likelihood value Double_t sigDPLike_; //! Signal likelihood from extra PDFs Double_t sigExtraLike_; Double_t sigFlavTagLike_; Double_t bkgdFlavTagLike_; //! Total signal likelihood Double_t sigTotalLike_; - //! Background DP likelihood value(s) + //! Background DP likelihood value(s) std::vector bkgndDPLike_; - //! Background likelihood value(s) from extra PDFs + //! Background likelihood value(s) from extra PDFs std::vector bkgndExtraLike_; //! Total background likelihood(s) std::vector bkgndTotalLike_; ClassDef(LauTimeDepFitModel,0) // Time-dependent neutral model }; #endif diff --git a/src/LauDecayTimePdf.cc b/src/LauDecayTimePdf.cc index 2ea493b..ce8a442 100644 --- a/src/LauDecayTimePdf.cc +++ b/src/LauDecayTimePdf.cc @@ -1,1594 +1,1595 @@ /* 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;} //effiTermN_ = effiFun_->getInt(); effiTermN_ = 1.0; } else { effiTerm_ = 1.0; effiTermN_ = 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) { 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 || type_ == SimFitNormBs) { 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_ = ARaw_->unblindValue(); Int_t curEvtTagFlv = flavTag_->getEvtTagFlvVals(iEvt); + // For the normalisation decays we know the true tag 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); } Int_t tagCat = flavTag_->getCurEvtTagCat(); Double_t At(0.0); if (tagCat != 0){ LauTagCatParamMap& tagAsyms = flavTag_->getTagAsym(); At = tagAsyms[tagCat].unblindValue(); } Double_t Ap = AProd_->unblindValue(); omegaTerm_ = ((1-Ap-At)*omegabar-(1+Ap+At)*omega); omegaTermh_ = ((1-Ap-At)*omegabar+(1+Ap+At)*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 || type_ == SimFitSigBs) { 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); } Int_t tagCat = flavTag_->getCurEvtTagCat(); Double_t At(0.0); if (tagCat != 0){ LauTagCatParamMap& tagAsyms = flavTag_->getTagAsym(); At = tagAsyms[tagCat].unblindValue(); } Double_t Ap = AProd_->unblindValue(); omegaTerm_ = ((1-Ap-At)*omegabar-(1+Ap+At)*omega); omegaTermh_ = ((1-Ap-At)*omegabar+(1+Ap+At)*omega); cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*omegaTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_*omegaTerm_; coshTerm_ = expTerm_*omegaTermh_; sinhTerm_ = 0.0; if (type_ == SimFitNormBs){ Double_t deltaGamma = deltaGamma_->value(); coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0); 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 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 || 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_ = 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); } Int_t tagCat = flavTag_->getCurEvtTagCat(); Double_t At(0.0); if (tagCat != 0){ LauTagCatParamMap& tagAsyms = flavTag_->getTagAsym(); At = tagAsyms[tagCat].unblindValue(); } Double_t Ap = AProd_->unblindValue(); omegaTerm_ = ((1-Ap-At)*omegabar-(1+Ap+At)*omega); omegaTermh_ = ((1-Ap-At)*omegabar+(1+Ap+At)*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 || 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); } Int_t tagCat = flavTag_->getCurEvtTagCat(); Double_t At(0.0); if (tagCat != 0){ LauTagCatParamMap& tagAsyms = flavTag_->getTagAsym(); At = tagAsyms[tagCat].unblindValue(); } Double_t Ap = AProd_->unblindValue(); omegaTerm_ = ((1-Ap-At)*omegabar-(1+Ap+At)*omega); omegaTermh_ = ((1-Ap-At)*omegabar+(1+Ap+At)*omega); cosTerm_ = TMath::Cos(deltaM*abscissa)*expTerm_*omegaTerm_; sinTerm_ = TMath::Sin(deltaM*abscissa)*expTerm_*omegaTerm_; coshTerm_ = expTerm_*omegaTermh_; sinhTerm_ = 0.0; if (type_ == SimFitNormBs){ Double_t deltaGamma = deltaGamma_->value(); coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0); 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;} //effiTermN_ = effiFun_->getInt(); effiTermN_ = 1.0; } else { effiTerm_ = 1.0; effiTermN_ = 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 effiPars) { if (effiPars.size() != effiFun_->getnKnots()){ cerr<<"ERROR in LauDecayTimePdf::updateEffiSpline : number of efficiency parameters is not equal to the number of spline knots."<Exit(EXIT_FAILURE); } effiFun_->updateYValues(effiPars); } void LauDecayTimePdf::setAsymmetries(LauParameter* AProd, LauParameter* ARaw) { AProd_ = AProd; ARaw_ = ARaw; } void LauDecayTimePdf::setFlavourTagging(LauFlavTag* flavTag) { flavTag_ = flavTag; } diff --git a/src/LauTimeDepFitModel.cc b/src/LauTimeDepFitModel.cc index c6f753c..c8c375d 100644 --- a/src/LauTimeDepFitModel.cc +++ b/src/LauTimeDepFitModel.cc @@ -1,2507 +1,2940 @@ /* 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 LauTimeDepFitModel.cc \brief File containing implementation of LauTimeDepFitModel class. */ #include #include #include #include #include #include "TFile.h" #include "TMinuit.h" #include "TRandom.h" #include "TSystem.h" #include "TVirtualFitter.h" #include "LauAbsBkgndDPModel.hh" #include "LauAbsCoeffSet.hh" #include "LauAbsPdf.hh" #include "LauAsymmCalc.hh" #include "LauComplex.hh" #include "LauConstants.hh" #include "LauDPPartialIntegralInfo.hh" #include "LauDaughters.hh" #include "LauDecayTimePdf.hh" #include "LauFitNtuple.hh" #include "LauGenNtuple.hh" #include "LauIsobarDynamics.hh" #include "LauKinematics.hh" #include "LauPrint.hh" #include "LauRandom.hh" #include "LauScfMap.hh" #include "LauTimeDepFitModel.hh" #include "LauFlavTag.hh" ClassImp(LauTimeDepFitModel) LauTimeDepFitModel::LauTimeDepFitModel(LauIsobarDynamics* modelB0bar, LauIsobarDynamics* modelB0, LauFlavTag* flavTag) : LauAbsFitModel(), sigModelB0bar_(modelB0bar), sigModelB0_(modelB0), kinematicsB0bar_(modelB0bar ? modelB0bar->getKinematics() : 0), kinematicsB0_(modelB0 ? modelB0->getKinematics() : 0), + usingBkgnd_(kFALSE), flavTag_(flavTag), nSigComp_(0), nSigDPPar_(0), nDecayTimePar_(0), nExtraPdfPar_(0), nNormPar_(0), nCalibPar_(0), nEffiPar_(0), nAsymPar_(0), nTagAsymPar_(0), coeffsB0bar_(0), coeffsB0_(0), coeffPars_(0), fitFracB0bar_(0), fitFracB0_(0), fitFracAsymm_(0), acp_(0), meanEffB0bar_("meanEffB0bar",0.0,0.0,1.0), meanEffB0_("meanEffB0",0.0,0.0,1.0), DPRateB0bar_("DPRateB0bar",0.0,0.0,100.0), DPRateB0_("DPRateB0",0.0,0.0,100.0), signalEvents_(0), signalAsym_(0), cpevVarName_(""), cpEigenValue_(CPEven), evtCPEigenVals_(0), deltaM_("deltaM",0.0), deltaGamma_("deltaGamma",0.0), tau_("tau",LauConstants::tauB0), phiMix_("phiMix", 2.0*LauConstants::beta, -LauConstants::threePi, LauConstants::threePi, kFALSE), sinPhiMix_("sinPhiMix", TMath::Sin(2.0*LauConstants::beta), -1.0, 1.0, kFALSE), cosPhiMix_("cosPhiMix", TMath::Cos(2.0*LauConstants::beta), -1.0, 1.0, kFALSE), useSinCos_(kFALSE), phiMixComplex_(TMath::Cos(-2.0*LauConstants::beta),TMath::Sin(-2.0*LauConstants::beta)), signalDecayTimePdfs_(), + backgroundDecayTimePdfs_(), curEvtDecayTime_(0.0), curEvtDecayTimeErr_(0.0), sigExtraPdf_(), sigFlavTagPdf_(), bkgdFlavTagPdf_(), AProd_("AProd",0.0,-1.0,1.0,kTRUE), ARaw_("ARaw",1.0,0.0,2.0,kTRUE), iterationsMax_(500000), nGenLoop_(0), ASq_(0.0), aSqMaxVar_(0.0), aSqMaxSet_(1.25), storeGenAmpInfo_(kFALSE), signalTree_(), reuseSignal_(kFALSE), sigDPLike_(0.0), sigExtraLike_(0.0), sigFlavTagLike_(0.0), bkgdFlavTagLike_(0.0), sigTotalLike_(0.0) { // Set up ftag here? // Make sure that the integration scheme will be symmetrised sigModelB0bar_->forceSymmetriseIntegration(kTRUE); sigModelB0_->forceSymmetriseIntegration(kTRUE); } LauTimeDepFitModel::~LauTimeDepFitModel() { for (LauTagCatEmbDataMap::iterator iter = signalTree_.begin(); iter != signalTree_.end(); ++iter){ delete iter->second; } for (LauTagCatEmbDataMapList::iterator iterlist = bkgndTree_.begin(); iterlist != bkgndTree_.end(); ++iterlist){ for (LauTagCatEmbDataMap::iterator iter = (*iterlist).begin(); iter != (*iterlist).end(); ++iter){ delete iter->second; } } } void LauTimeDepFitModel::setupBkgndVectors() { UInt_t nBkgnds = this->nBkgndClasses(); - BkgndDPModelsB0_.resize( nBkgnds ); - BkgndDPModelsB0bar_.resize( nBkgnds ); - BkgndPdfsB0_.resize( nBkgnds ); - BkgndPdfsB0bar_.resize( nBkgnds ); + BkgndDPModels_.resize( nBkgnds ); + BkgndPdfs_.resize( nBkgnds ); bkgndEvents_.resize( nBkgnds ); bkgndAsym_.resize( nBkgnds ); bkgndTree_.resize( nBkgnds ); reuseBkgnd_.resize( nBkgnds ); bkgndDPLike_.resize( nBkgnds ); bkgndExtraLike_.resize( nBkgnds ); bkgndTotalLike_.resize( nBkgnds ); } void LauTimeDepFitModel::setNSigEvents(LauParameter* nSigEvents) { if ( nSigEvents == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( signalEvents_ != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl; return; } if ( signalAsym_ != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl; return; } signalEvents_ = nSigEvents; signalEvents_->name("signalEvents"); Double_t value = nSigEvents->value(); signalEvents_->range(-2.0*(TMath::Abs(value)+1.0),2.0*(TMath::Abs(value)+1.0)); signalAsym_ = new LauParameter("signalAsym",0.0,-1.0,1.0,kTRUE); } void LauTimeDepFitModel::setNSigEvents(LauParameter* nSigEvents, LauParameter* sigAsym) { if ( nSigEvents == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The event LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( sigAsym == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The asym LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( signalEvents_ != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl; return; } if ( signalAsym_ != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl; return; } signalEvents_ = nSigEvents; signalEvents_->name("signalEvents"); Double_t value = nSigEvents->value(); signalEvents_->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0)); signalAsym_ = sigAsym; signalAsym_->name("signalAsym"); signalAsym_->range(-1.0,1.0); } void LauTimeDepFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents) { if ( nBkgndEvents == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBgkndEvents : The background yield LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; gSystem->Exit(EXIT_FAILURE); } UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() ); if ( bkgndEvents_[bkgndID] != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl; return; } if ( bkgndAsym_[bkgndID] != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl; return; } nBkgndEvents->name( nBkgndEvents->name()+"Events" ); if ( nBkgndEvents->isLValue() ) { Double_t value = nBkgndEvents->value(); LauParameter* yield = dynamic_cast( nBkgndEvents ); yield->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0)); } bkgndEvents_[bkgndID] = nBkgndEvents; bkgndAsym_[bkgndID] = new LauParameter(nBkgndEvents->name()+"Asym",0.0,-1.0,1.0,kTRUE); } void LauTimeDepFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym) { if ( nBkgndEvents == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : The background yield LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( bkgndAsym == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : The background asym LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; gSystem->Exit(EXIT_FAILURE); } UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() ); if ( bkgndEvents_[bkgndID] != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl; return; } if ( bkgndAsym_[bkgndID] != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl; return; } bkgndEvents_[bkgndID]->name( nBkgndEvents->name()+"Events" ); if ( nBkgndEvents->isLValue() ) { Double_t value = nBkgndEvents->value(); LauParameter* yield = dynamic_cast( nBkgndEvents ); yield->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0)); } bkgndEvents_[bkgndID] = nBkgndEvents; bkgndAsym_[bkgndID]->name( nBkgndEvents->name()+"Asym" ); if ( bkgndAsym->isLValue() ) { LauParameter* asym = dynamic_cast( bkgndAsym ); asym->range(-1.0, 1.0); } bkgndAsym_[bkgndID] = bkgndAsym; } void LauTimeDepFitModel::setSignalDtPdf(Int_t tagCat, LauDecayTimePdf* pdf) { if (!flavTag_->validTagCat(tagCat)) { std::cerr<<"ERROR in LauTimeDepFitModel::setSignalDtPdf : Tagging category \""<setAsymmetries(&AProd_, &ARaw_); } +void LauTimeDepFitModel::setBackgroundDtPdf(Int_t tagCat, LauDecayTimePdf* pdf) +{ + // TODO If these are all histograms shouldn't need to add much more code in other functions + if (!flavTag_->validTagCat(tagCat)) { + std::cerr<<"ERROR in LauTimeDepFitModel::setBackgroundDtPdf : Tagging category \""<setAsymmetries(&AProd_, &ARaw_); + usingBkgnd_ = kTRUE; +} + +void LauTimeDepFitModel::setBkgndDPModels(const TString& bkgndClass, LauAbsBkgndDPModel* model) +{ + if (model==0) { + std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDPModels : the model pointer is null." << std::endl; + return; + } + + // check that this background name is valid + if ( ! this->validBkgndClass( bkgndClass) ) { + std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDPModel : Invalid background class \"" << bkgndClass << "\"." << std::endl; + std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; + return; + } + + UInt_t bkgndID = this->bkgndClassID( bkgndClass ); + BkgndDPModels_[bkgndID] = model; + + usingBkgnd_ = kTRUE; +} + void LauTimeDepFitModel::setSignalPdfs(Int_t tagCat, LauAbsPdf* pdf) { // These "extra variables" are assumed to be purely kinematical, like mES and DeltaE //or making use of Rest of Event information, and therefore independent of whether //the parent is a B0 or a B0bar. If this assupmtion doesn't hold, do modify this part! if (!flavTag_->validTagCat(tagCat)) { std::cerr<<"ERROR in LauTimeDepFitModel::setSignalPdfs : Tagging category \""<validBkgndClass( bkgndClass ) ) { + std::cerr << "ERROR in LauTimeDepFitModel::setBkgndPdf : Invalid background class \"" << bkgndClass << "\"." << std::endl; + std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; + return; + } + + UInt_t bkgndID = this->bkgndClassID( bkgndClass ); + BkgndPdfs_[bkgndID].push_back(pdf); + + usingBkgnd_ = kTRUE; +} + void LauTimeDepFitModel::setPhiMix(const Double_t phiMix, const Bool_t fixPhiMix, const Bool_t useSinCos) { phiMix_.value(phiMix); phiMix_.initValue(phiMix); phiMix_.genValue(phiMix); phiMix_.fixed(fixPhiMix); const Double_t sinPhiMix = TMath::Sin(phiMix); sinPhiMix_.value(sinPhiMix); sinPhiMix_.initValue(sinPhiMix); sinPhiMix_.genValue(sinPhiMix); sinPhiMix_.fixed(fixPhiMix); const Double_t cosPhiMix = TMath::Cos(phiMix); cosPhiMix_.value(cosPhiMix); cosPhiMix_.initValue(cosPhiMix); cosPhiMix_.genValue(cosPhiMix); cosPhiMix_.fixed(fixPhiMix); useSinCos_ = useSinCos; phiMixComplex_.setRealPart(cosPhiMix); phiMixComplex_.setImagPart(-1.0*sinPhiMix); } void LauTimeDepFitModel::initialise() { // From the initial parameter values calculate the coefficients // so they can be passed to the signal model this->updateCoeffs(); // Initialisation if (this->useDP() == kTRUE) { this->initialiseDPModels(); } //Flavour tagging flavTag_->initialise(); if (!this->useDP() && sigExtraPdf_.empty()) { std::cerr<<"ERROR in LauTimeDepFitModel::initialise : Signal model doesn't exist for any variable."<Exit(EXIT_FAILURE); } if (this->useDP() == kTRUE) { // Check that we have all the Dalitz-plot models if ((sigModelB0bar_ == 0) || (sigModelB0_ == 0)) { std::cerr<<"ERROR in LauTimeDepFitModel::initialise : the pointer to one (particle or anti-particle) of the signal DP models is null."<Exit(EXIT_FAILURE); } } + // Next check that, if a given component is being used we've got the + // right number of PDFs for all the variables involved + // TODO - should probably check variable names and so on as well + + //UInt_t nsigpdfvars(0); + //for ( LauPdfList::const_iterator pdf_iter = sigExtraPdf_[tagCat]_.begin(); pdf_iter != sigExtraPdf_[tagCat]_.end(); ++pdf_iter ) { + // std::vector varNames = (*pdf_iter)->varNames(); + // for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { + // if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { + // ++nsigpdfvars; + // } + // } + //} + + //if (usingBkgnd_) { + // for (LauBkgndPdfsList::const_iterator bgclass_iter = BkgndPdfsB0_.begin(); bgclass_iter != BkgndPdfsB0_.end(); ++bgclass_iter) { + // UInt_t nbkgndpdfvars(0); + // const LauPdfList& pdfList = (*bgclass_iter); + // for ( LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter ) { + // std::vector varNames = (*pdf_iter)->varNames(); + // for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { + // if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { + // ++nbkgndpdfvars; + // } + // } + // } + // if (nbkgndpdfvars != nsigpdfvars) { + // std::cerr << "ERROR in LauTimeDepFitModel::initialise : There are " << nsigpdfvars << " signal PDF variables but " << nbkgndpdfvars << " bkgnd PDF variables." << std::endl; + // gSystem->Exit(EXIT_FAILURE); + // } + // } + //} + // Clear the vectors of parameter information so we can start from scratch this->clearFitParVectors(); // Set the fit parameters for signal and background models this->setSignalDPParameters(); // Set the fit parameters for the decay time models this->setDecayTimeParameters(); // Set the fit parameters for the extra PDFs this->setExtraPdfParameters(); // Set the initial bg and signal events this->setFitNEvents(); // Handle flavour-tagging calibration parameters this->setCalibParams(); // Add the efficiency parameters this->setEffiParams(); //Asymmetry terms AProd and ARaw added in setAsymmetries()? this->setAsymParams(); // Check that we have the expected number of fit variables const LauParameterPList& fitVars = this->fitPars(); if (fitVars.size() != (nSigDPPar_ + nDecayTimePar_ + nExtraPdfPar_ + nNormPar_ + nCalibPar_ + nEffiPar_ + nAsymPar_ + nTagAsymPar_)) { std::cerr<<"ERROR in LauTimeDepFitModel::initialise : Number of fit parameters not of expected size."<Exit(EXIT_FAILURE); } this->setExtraNtupleVars(); } void LauTimeDepFitModel::recalculateNormalisation() { sigModelB0bar_->recalculateNormalisation(); sigModelB0_->recalculateNormalisation(); sigModelB0bar_->modifyDataTree(); sigModelB0_->modifyDataTree(); this->calcInterferenceTermIntegrals(); } void LauTimeDepFitModel::initialiseDPModels() { if (sigModelB0bar_ == 0) { std::cerr<<"ERROR in LauTimeDepFitModel::initialiseDPModels : B0bar signal DP model doesn't exist"<Exit(EXIT_FAILURE); } if (sigModelB0_ == 0) { std::cerr<<"ERROR in LauTimeDepFitModel::initialiseDPModels : B0 signal DP model doesn't exist"<Exit(EXIT_FAILURE); } // Need to check that the number of components we have and that the dynamics has matches up const UInt_t nAmpB0bar = sigModelB0bar_->getnTotAmp(); const UInt_t nAmpB0 = sigModelB0_->getnTotAmp(); if ( nAmpB0bar != nAmpB0 ) { std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Unequal number of signal DP components in the particle and anti-particle models: " << nAmpB0bar << " != " << nAmpB0 << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( nAmpB0bar != nSigComp_ ) { std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Number of signal DP components in the model (" << nAmpB0bar << ") not equal to number of coefficients supplied (" << nSigComp_ << ")." << std::endl; gSystem->Exit(EXIT_FAILURE); } std::cout<<"INFO in LauTimeDepFitModel::initialiseDPModels : Initialising signal DP model"<initialise(coeffsB0bar_); sigModelB0_->initialise(coeffsB0_); fifjEffSum_.clear(); fifjEffSum_.resize(nSigComp_); for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) { fifjEffSum_[iAmp].resize(nSigComp_); } // calculate the integrals of the A*Abar terms this->calcInterferenceTermIntegrals(); this->calcInterTermNorm(); + // Add backgrounds + if (usingBkgnd_ == kTRUE) { + for (LauBkgndDPModelList::iterator iter = BkgndDPModels_.begin(); iter != BkgndDPModels_.end(); ++iter) { + (*iter)->initialise(); + } + } + } void LauTimeDepFitModel::calcInterferenceTermIntegrals() { const std::vector& integralInfoListB0bar = sigModelB0bar_->getIntegralInfos(); const std::vector& integralInfoListB0 = sigModelB0_->getIntegralInfos(); // TODO should check (first time) that they match in terms of number of entries in the vectors and that each entry has the same number of points, ranges, weights etc. LauComplex A, Abar, fifjEffSumTerm; for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) { for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) { fifjEffSum_[iAmp][jAmp].zero(); } } const UInt_t nIntegralRegions = integralInfoListB0bar.size(); for ( UInt_t iRegion(0); iRegion < nIntegralRegions; ++iRegion ) { const LauDPPartialIntegralInfo* integralInfoB0bar = integralInfoListB0bar[iRegion]; const LauDPPartialIntegralInfo* integralInfoB0 = integralInfoListB0[iRegion]; const UInt_t nm13Points = integralInfoB0bar->getnm13Points(); const UInt_t nm23Points = integralInfoB0bar->getnm23Points(); for (UInt_t m13 = 0; m13 < nm13Points; ++m13) { for (UInt_t m23 = 0; m23 < nm23Points; ++m23) { const Double_t weight = integralInfoB0bar->getWeight(m13,m23); const Double_t eff = integralInfoB0bar->getEfficiency(m13,m23); const Double_t effWeight = eff*weight; for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) { A = integralInfoB0->getAmplitude(m13, m23, iAmp); for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) { Abar = integralInfoB0bar->getAmplitude(m13, m23, jAmp); fifjEffSumTerm = Abar*A.conj(); fifjEffSumTerm.rescale(effWeight); fifjEffSum_[iAmp][jAmp] += fifjEffSumTerm; } } } } } } void LauTimeDepFitModel::calcInterTermNorm() { const std::vector& fNormB0bar = sigModelB0bar_->getFNorm(); const std::vector& fNormB0 = sigModelB0_->getFNorm(); LauComplex norm; for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) { for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) { LauComplex coeffTerm = coeffsB0bar_[jAmp]*coeffsB0_[iAmp].conj(); coeffTerm *= fifjEffSum_[iAmp][jAmp]; coeffTerm.rescale(fNormB0bar[jAmp] * fNormB0[iAmp]); norm += coeffTerm; } } norm *= phiMixComplex_; interTermReNorm_ = 2.0*norm.re(); interTermImNorm_ = 2.0*norm.im(); } void LauTimeDepFitModel::setAmpCoeffSet(LauAbsCoeffSet* coeffSet) { // Is there a component called compName in the signal models? TString compName = coeffSet->name(); TString conjName = sigModelB0bar_->getConjResName(compName); const LauDaughters* daughtersB0bar = sigModelB0bar_->getDaughters(); const LauDaughters* daughtersB0 = sigModelB0_->getDaughters(); const Bool_t conjugate = daughtersB0bar->isConjugate( daughtersB0 ); if ( ! sigModelB0bar_->hasResonance(compName) ) { if ( ! sigModelB0bar_->hasResonance(conjName) ) { std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0bar signal DP model doesn't contain component \""<name( compName ); } if ( conjugate ) { if ( ! sigModelB0_->hasResonance(conjName) ) { std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<hasResonance(compName) ) { std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<::const_iterator iter=coeffPars_.begin(); iter!=coeffPars_.end(); ++iter) { if ((*iter)->name() == compName) { std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : Have already set coefficients for \""<index(nSigComp_); coeffPars_.push_back(coeffSet); TString parName = coeffSet->baseName(); parName += "FitFracAsym"; fitFracAsymm_.push_back(LauParameter(parName, 0.0, -1.0, 1.0)); acp_.push_back(coeffSet->acp()); ++nSigComp_; std::cout<<"INFO in LauTimeDepFitModel::setAmpCoeffSet : Added coefficients for component \""<acp(); LauAsymmCalc asymmCalc(fitFracB0bar_[i][i].value(), fitFracB0_[i][i].value()); Double_t asym = asymmCalc.getAsymmetry(); fitFracAsymm_[i].value(asym); if (initValues) { fitFracAsymm_[i].genValue(asym); fitFracAsymm_[i].initValue(asym); } } } void LauTimeDepFitModel::setSignalDPParameters() { // Set the fit parameters for the signal model. nSigDPPar_ = 0; if ( ! this->useDP() ) { return; } std::cout << "INFO in LauTimeDepFitModel::setSignalDPParameters : Setting the initial fit parameters for the signal DP model." << std::endl; // Place isobar coefficient parameters in vector of fit variables LauParameterPList& fitVars = this->fitPars(); for (UInt_t i = 0; i < nSigComp_; ++i) { LauParameterPList pars = coeffPars_[i]->getParameters(); for (LauParameterPList::iterator iter = pars.begin(); iter != pars.end(); ++iter) { if ( !(*iter)->clone() ) { fitVars.push_back(*iter); ++nSigDPPar_; } } } // Obtain the resonance parameters and place them in the vector of fit variables and in a separate vector // Need to make sure that they are unique because some might appear in both DP models LauParameterPSet& resVars = this->resPars(); resVars.clear(); LauParameterPList& sigDPParsB0bar = sigModelB0bar_->getFloatingParameters(); LauParameterPList& sigDPParsB0 = sigModelB0_->getFloatingParameters(); for ( LauParameterPList::iterator iter = sigDPParsB0bar.begin(); iter != sigDPParsB0bar.end(); ++iter ) { if ( resVars.insert(*iter).second ) { fitVars.push_back(*iter); ++nSigDPPar_; } } for ( LauParameterPList::iterator iter = sigDPParsB0.begin(); iter != sigDPParsB0.end(); ++iter ) { if ( resVars.insert(*iter).second ) { fitVars.push_back(*iter); ++nSigDPPar_; } } } UInt_t LauTimeDepFitModel::addParametersToFitList(LauTagCatDtPdfMap& theMap) { UInt_t counter(0); LauParameterPList& fitVars = this->fitPars(); // loop through the map for (LauTagCatDtPdfMap::iterator iter = theMap.begin(); iter != theMap.end(); ++iter) { // grab the pdf and then its parameters LauDecayTimePdf* thePdf = (*iter).second; // The first one is the tagging category LauAbsRValuePList& rvalues = thePdf->getParameters(); // loop through the parameters for (LauAbsRValuePList::iterator pars_iter = rvalues.begin(); pars_iter != rvalues.end(); ++pars_iter) { LauParameterPList params = (*pars_iter)->getPars(); for (LauParameterPList::iterator params_iter = params.begin(); params_iter != params.end(); ++params_iter) { // for each "original" parameter add it to the list of fit parameters and increment the counter if ( !(*params_iter)->clone() && ( !(*params_iter)->fixed() || (this->twoStageFit() && (*params_iter)->secondStage()) ) ) { fitVars.push_back(*params_iter); ++counter; } } } } return counter; } UInt_t LauTimeDepFitModel::addParametersToFitList(LauTagCatPdfListMap& theMap) { UInt_t counter(0); // loop through the map for (LauTagCatPdfListMap::iterator iter = theMap.begin(); iter != theMap.end(); ++iter) { counter += this->addFitParameters(iter->second); // first is the tagging category } return counter; } void LauTimeDepFitModel::setDecayTimeParameters() { nDecayTimePar_ = 0; // Loop over the Dt PDFs nDecayTimePar_ += this->addParametersToFitList(signalDecayTimePdfs_); + if (usingBkgnd_){ + nDecayTimePar_ += this->addParametersToFitList(backgroundDecayTimePdfs_); + } + LauParameterPList& fitVars = this->fitPars(); if (useSinCos_) { fitVars.push_back(&sinPhiMix_); fitVars.push_back(&cosPhiMix_); nDecayTimePar_ += 2; } else { fitVars.push_back(&phiMix_); ++nDecayTimePar_; } } void LauTimeDepFitModel::setExtraPdfParameters() { // Include the parameters of the PDF for each tagging category in the fit // NB all of them are passed to the fit, even though some have been fixed through parameter.fixed(kTRUE) // With the new "cloned parameter" scheme only "original" parameters are passed to the fit. // Their clones are updated automatically when the originals are updated. nExtraPdfPar_ = 0; nExtraPdfPar_ += this->addParametersToFitList(sigExtraPdf_); + + if (usingBkgnd_ == kTRUE) { + for (LauBkgndPdfsList::iterator iter = BkgndPdfs_.begin(); iter != BkgndPdfs_.end(); ++iter) { + nExtraPdfPar_ += this->addFitParameters(*iter); + } + } } void LauTimeDepFitModel::setFitNEvents() { nNormPar_ = 0; // Initialise the total number of events to be the sum of all the hypotheses Double_t nTotEvts = signalEvents_->value(); this->eventsPerExpt(TMath::FloorNint(nTotEvts)); LauParameterPList& fitVars = this->fitPars(); // if doing an extended ML fit add the signal fraction into the fit parameters if (this->doEMLFit()) { std::cout<<"INFO in LauTimeDepFitModel::setFitNEvents : Initialising number of events for signal and background components..."<useDP() == kFALSE) { fitVars.push_back(signalAsym_); ++nNormPar_; } // TODO arguably should delegate this LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac(); // tagging-category fractions for signal events for (LauTagCatParamMap::iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) { if (iter == signalTagCatFrac.begin()) { continue; } LauParameter* par = &((*iter).second); fitVars.push_back(par); ++nNormPar_; } + + // Backgrounds + if (usingBkgnd_ == kTRUE) { + for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { + std::vector parameters = (*iter)->getPars(); + for ( LauParameter* parameter : parameters ) { + if(!parameter->clone()) { + fitVars.push_back(parameter); + ++nNormPar_; + } + } + } + for (LauBkgndYieldList::iterator iter = bkgndAsym_.begin(); iter != bkgndAsym_.end(); ++iter) { + std::vector parameters = (*iter)->getPars(); + for ( LauParameter* parameter : parameters ) { + if(!parameter->clone()) { + fitVars.push_back(parameter); + ++nNormPar_; + } + } + } + } } void LauTimeDepFitModel::setAsymParams() { LauParameterPList& fitVars = this->fitPars(); fitVars.push_back(&AProd_); fitVars.push_back(&ARaw_); nAsymPar_+=2; } void LauTimeDepFitModel::setCalibParams() { Bool_t useAltPars = flavTag_->getUseAveOmegaDeltaOmega(); if (useAltPars){ LauTagCatParamMap& p0pars_ave = flavTag_->getCalibP0Ave(); LauTagCatParamMap& p0pars_delta = flavTag_->getCalibP0Delta(); LauTagCatParamMap& p1pars_ave = flavTag_->getCalibP1Ave(); LauTagCatParamMap& p1pars_delta = flavTag_->getCalibP1Delta(); LauParameterPList& fitVars = this->fitPars(); for(LauTagCatParamMap::iterator iter = p0pars_ave.begin(); iter != p0pars_ave.end(); ++iter){ LauParameter* p0 = &((*iter).second); if (p0->fixed()){continue;} fitVars.push_back(p0); ++nCalibPar_; } for(LauTagCatParamMap::iterator iter = p0pars_delta.begin(); iter != p0pars_delta.end(); ++iter){ LauParameter* p0 = &((*iter).second); if (p0->fixed()){continue;} fitVars.push_back(p0); ++nCalibPar_; } for(LauTagCatParamMap::iterator iter = p1pars_ave.begin(); iter != p1pars_ave.end(); ++iter){ LauParameter* p1 = &((*iter).second); if (p1->fixed()){continue;} fitVars.push_back(p1); ++nCalibPar_; } for(LauTagCatParamMap::iterator iter = p1pars_delta.begin(); iter != p1pars_delta.end(); ++iter){ LauParameter* p1 = &((*iter).second); if (p1->fixed()){continue;} fitVars.push_back(p1); ++nCalibPar_; } } else { LauTagCatParamMap& p0pars_b0 = flavTag_->getCalibP0B0(); LauTagCatParamMap& p0pars_b0bar = flavTag_->getCalibP0B0bar(); LauTagCatParamMap& p1pars_b0 = flavTag_->getCalibP1B0(); LauTagCatParamMap& p1pars_b0bar = flavTag_->getCalibP1B0bar(); LauParameterPList& fitVars = this->fitPars(); for(LauTagCatParamMap::iterator iter = p0pars_b0.begin(); iter != p0pars_b0.end(); ++iter){ LauParameter* p0 = &((*iter).second); if (p0->fixed()){continue;} fitVars.push_back(p0); ++nCalibPar_; } for(LauTagCatParamMap::iterator iter = p0pars_b0bar.begin(); iter != p0pars_b0bar.end(); ++iter){ LauParameter* p0 = &((*iter).second); if (p0->fixed()){continue;} fitVars.push_back(p0); ++nCalibPar_; } for(LauTagCatParamMap::iterator iter = p1pars_b0.begin(); iter != p1pars_b0.end(); ++iter){ LauParameter* p1 = &((*iter).second); if (p1->fixed()){continue;} fitVars.push_back(p1); ++nCalibPar_; } for(LauTagCatParamMap::iterator iter = p1pars_b0bar.begin(); iter != p1pars_b0bar.end(); ++iter){ LauParameter* p1 = &((*iter).second); if (p1->fixed()){continue;} fitVars.push_back(p1); ++nCalibPar_; } } } void LauTimeDepFitModel::setFlavTagAsymParams(){ LauParameterPList& fitVars = this->fitPars(); LauTagCatParamMap& tagAsyms = flavTag_ -> getTagAsym(); for(LauTagCatParamMap::iterator iter = tagAsyms.begin(); iter != tagAsyms.end(); ++iter){ LauParameter* At = &((*iter).second); if (At->fixed()){continue;} fitVars.push_back(At); ++nTagAsymPar_; } } void LauTimeDepFitModel::setEffiParams() { LauParameterPList& fitVars = this->fitPars(); LauTagCatDtPdfMap& signalDecayTimePdfs = this->getSignalDecayTimePdfs(); // Loop over DecayTimePdfs for (LauTagCatDtPdfMap::iterator iter = signalDecayTimePdfs.begin(); iter != signalDecayTimePdfs.end(); ++iter){ LauDecayTimePdf* pdf = (*iter).second; std::vector& effiPars = pdf->getEffiPars(); for(std::vector::iterator iter2 = effiPars.begin(); iter2 != effiPars.end(); ++iter2){ LauParameter* par = *iter2; if (par->fixed()){continue;} fitVars.push_back(par); ++nEffiPar_; } } } void LauTimeDepFitModel::setExtraNtupleVars() { // Set-up other parameters derived from the fit results, e.g. fit fractions. if (this->useDP() != kTRUE) { return; } // First clear the vectors so we start from scratch this->clearExtraVarVectors(); LauParameterList& extraVars = this->extraPars(); // Add the B0 and B0bar fit fractions for each signal component fitFracB0bar_ = sigModelB0bar_->getFitFractions(); if (fitFracB0bar_.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); igetFitFractions(); if (fitFracB0_.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); icalcAsymmetries(kTRUE); // Add the Fit Fraction asymmetry for each signal component for (UInt_t i = 0; i < nSigComp_; i++) { extraVars.push_back(fitFracAsymm_[i]); } // Add the calculated CP asymmetry for each signal component for (UInt_t i = 0; i < nSigComp_; i++) { extraVars.push_back(acp_[i]); } // Now add in the DP efficiency values Double_t initMeanEffB0bar = sigModelB0bar_->getMeanEff().initValue(); meanEffB0bar_.value(initMeanEffB0bar); meanEffB0bar_.initValue(initMeanEffB0bar); meanEffB0bar_.genValue(initMeanEffB0bar); extraVars.push_back(meanEffB0bar_); Double_t initMeanEffB0 = sigModelB0_->getMeanEff().initValue(); meanEffB0_.value(initMeanEffB0); meanEffB0_.initValue(initMeanEffB0); meanEffB0_.genValue(initMeanEffB0); extraVars.push_back(meanEffB0_); // Also add in the DP rates Double_t initDPRateB0bar = sigModelB0bar_->getDPRate().initValue(); DPRateB0bar_.value(initDPRateB0bar); DPRateB0bar_.initValue(initDPRateB0bar); DPRateB0bar_.genValue(initDPRateB0bar); extraVars.push_back(DPRateB0bar_); Double_t initDPRateB0 = sigModelB0_->getDPRate().initValue(); DPRateB0_.value(initDPRateB0); DPRateB0_.initValue(initDPRateB0); DPRateB0_.genValue(initDPRateB0); extraVars.push_back(DPRateB0_); } void LauTimeDepFitModel::setAsymmetries(const Double_t AProd, const Bool_t AProdFix, const Double_t ARaw, const Bool_t ARawFix){ AProd_.value(AProd); AProd_.fixed(AProdFix); ARaw_.value(ARaw); ARaw_.fixed(ARawFix); //this->setAsymParams(); } void LauTimeDepFitModel::finaliseFitResults(const TString& tablePrefixName) { // Retrieve parameters from the fit results for calculations and toy generation // and eventually store these in output root ntuples/text files // Now take the fit parameters and update them as necessary // i.e. to make mag > 0.0, phase in the right range. // This function will also calculate any other values, such as the // fit fractions, using any errors provided by fitParErrors as appropriate. // Also obtain the pull values: (measured - generated)/(average error) if (this->useDP() == kTRUE) { for (UInt_t i = 0; i < nSigComp_; ++i) { // Check whether we have "a > 0.0", and phases in the right range coeffPars_[i]->finaliseValues(); } } // update the pulls on the event fractions and asymmetries if (this->doEMLFit()) { signalEvents_->updatePull(); } if (this->useDP() == kFALSE) { signalAsym_->updatePull(); } // Finalise the pulls on the decay time parameters for (LauTagCatDtPdfMap::iterator iter = signalDecayTimePdfs_.begin(); iter != signalDecayTimePdfs_.end(); ++iter) { LauDecayTimePdf* pdf = (*iter).second; pdf->updatePulls(); } + // and for backgrounds if required + if (usingBkgnd_){ + for (LauTagCatDtPdfMap::iterator iter = backgroundDecayTimePdfs_.begin(); iter != backgroundDecayTimePdfs_.end(); ++iter) { + LauDecayTimePdf* pdf = (*iter).second; + pdf->updatePulls(); + } + } + if (useSinCos_) { cosPhiMix_.updatePull(); sinPhiMix_.updatePull(); } else { this->checkMixingPhase(); } + if (usingBkgnd_ == kTRUE) { + for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { + std::vector parameters = (*iter)->getPars(); + for ( LauParameter* parameter : parameters ) { + parameter->updatePull(); + } + } + for (LauBkgndYieldList::iterator iter = bkgndAsym_.begin(); iter != bkgndAsym_.end(); ++iter) { + std::vector parameters = (*iter)->getPars(); + for ( LauParameter* parameter : parameters ) { + parameter->updatePull(); + } + } + } + // Update the pulls on all the extra PDFs' parameters for (LauTagCatPdfListMap::iterator iter = sigExtraPdf_.begin(); iter != sigExtraPdf_.end(); ++iter) { this->updateFitParameters(iter->second); } + if (usingBkgnd_ == kTRUE) { + for (LauBkgndPdfsList::iterator iter = BkgndPdfs_.begin(); iter != BkgndPdfs_.end(); ++iter) { + this->updateFitParameters(*iter); + } + } LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac(); // Tagging-category fractions for signal and background events Double_t firstCatFrac(1.0); Int_t firstCat(0); for (LauTagCatParamMap::iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) { if (iter == signalTagCatFrac.begin()) { firstCat = iter->first; continue; } LauParameter& par = (*iter).second; firstCatFrac -= par.value(); // update the parameter pull par.updatePull(); } signalTagCatFrac[firstCat].value(firstCatFrac); signalTagCatFrac[firstCat].updatePull(); // Fill the fit results to the ntuple // update the coefficients and then calculate the fit fractions and ACP's if (this->useDP() == kTRUE) { this->updateCoeffs(); sigModelB0bar_->updateCoeffs(coeffsB0bar_); sigModelB0bar_->calcExtraInfo(); sigModelB0_->updateCoeffs(coeffsB0_); sigModelB0_->calcExtraInfo(); LauParArray fitFracB0bar = sigModelB0bar_->getFitFractions(); if (fitFracB0bar.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } LauParArray fitFracB0 = sigModelB0_->getFitFractions(); if (fitFracB0.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); igetMeanEff().value()); meanEffB0_.value(sigModelB0_->getMeanEff().value()); DPRateB0bar_.value(sigModelB0bar_->getDPRate().value()); DPRateB0_.value(sigModelB0_->getDPRate().value()); this->calcAsymmetries(); // Then store the final fit parameters, and any extra parameters for // the signal model (e.g. fit fractions, FF asymmetries, ACPs, mean efficiency and DP rate) this->clearExtraVarVectors(); LauParameterList& extraVars = this->extraPars(); for (UInt_t i(0); iprintFitFractions(std::cout); this->printAsymmetries(std::cout); } const LauParameterPList& fitVars = this->fitPars(); const LauParameterList& extraVars = this->extraPars(); LauFitNtuple* ntuple = this->fitNtuple(); ntuple->storeParsAndErrors(fitVars, extraVars); // find out the correlation matrix for the parameters ntuple->storeCorrMatrix(this->iExpt(), this->fitStatus(), this->covarianceMatrix()); // Fill the data into ntuple ntuple->updateFitNtuple(); // Print out the partial fit fractions, phases and the // averaged efficiency, reweighted by the dynamics (and anything else) if (this->writeLatexTable()) { TString sigOutFileName(tablePrefixName); sigOutFileName += "_"; sigOutFileName += this->iExpt(); sigOutFileName += "Expt.tex"; this->writeOutTable(sigOutFileName); } } void LauTimeDepFitModel::printFitFractions(std::ostream& output) { // Print out Fit Fractions, total DP rate and mean efficiency // First for the B0bar events for (UInt_t i = 0; i < nSigComp_; i++) { const TString compName(coeffPars_[i]->name()); output<<"B0bar FitFraction for component "<useDP() == kTRUE) { // print the fit coefficients in one table coeffPars_.front()->printTableHeading(fout); for (UInt_t i = 0; i < nSigComp_; i++) { coeffPars_[i]->printTableRow(fout); } fout<<"\\hline"<name(); resName = resName.ReplaceAll("_", "\\_"); fout< =$ & $"; print.printFormat(fout, meanEffB0bar_.value()); fout << "$ & $"; print.printFormat(fout, meanEffB0_.value()); fout << "$ & & \\\\" << std::endl; if (useSinCos_) { fout << "$\\sinPhiMix =$ & $"; print.printFormat(fout, sinPhiMix_.value()); fout << " \\pm "; print.printFormat(fout, sinPhiMix_.error()); fout << "$ & & & & & & & \\\\" << std::endl; fout << "$\\cosPhiMix =$ & $"; print.printFormat(fout, cosPhiMix_.value()); fout << " \\pm "; print.printFormat(fout, cosPhiMix_.error()); fout << "$ & & & & & & & \\\\" << std::endl; } else { fout << "$\\phiMix =$ & $"; print.printFormat(fout, phiMix_.value()); fout << " \\pm "; print.printFormat(fout, phiMix_.error()); fout << "$ & & & & & & & \\\\" << std::endl; } fout << "\\hline \n\\end{tabular}" << std::endl; } if (!sigExtraPdf_.empty()) { fout<<"\\begin{tabular}{|l|c|}"<printFitParameters(iter->second, fout); } + if (usingBkgnd_ == kTRUE && !BkgndPdfs_.empty()) { + fout << "\\hline" << std::endl; + fout << "\\Extra Background PDFs' Parameters: & \\\\" << std::endl; + for (LauBkgndPdfsList::const_iterator iter = BkgndPdfs_.begin(); iter != BkgndPdfs_.end(); ++iter) { + this->printFitParameters(*iter, fout); + } + } fout<<"\\hline \n\\end{tabular}"<updateSigEvents(); // Check whether we want to have randomised initial fit parameters for the signal model if (this->useRandomInitFitPars() == kTRUE) { this->randomiseInitFitPars(); } } void LauTimeDepFitModel::randomiseInitFitPars() { // Only randomise those parameters that are not fixed! std::cout<<"INFO in LauTimeDepFitModel::randomiseInitFitPars : Randomising the initial values of the coefficients of the DP components (and phiMix)..."<randomiseInitValues(); } phiMix_.randomiseValue(-LauConstants::pi, LauConstants::pi); if (useSinCos_) { sinPhiMix_.initValue(TMath::Sin(phiMix_.initValue())); cosPhiMix_.initValue(TMath::Cos(phiMix_.initValue())); } } LauTimeDepFitModel::LauGenInfo LauTimeDepFitModel::eventsToGenerate() { // Determine the number of events to generate for each hypothesis // If we're smearing then smear each one individually // NB this individual smearing has to be done individually per tagging category as well LauGenInfo nEvtsGen; LauTagCatGenInfo eventsB0, eventsB0bar; // Signal // If we're including the DP and decay time we can't decide on the tag // yet, since it depends on the whole DP+dt PDF, however, if // we're not then we need to decide. Double_t evtWeight(1.0); Double_t nEvts = signalEvents_->genValue(); if ( nEvts < 0.0 ) { evtWeight = -1.0; nEvts = TMath::Abs( nEvts ); } LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac(); Double_t sigAsym(0.0); if (this->useDP() == kFALSE) { sigAsym = signalAsym_->genValue(); for (LauTagCatParamMap::const_iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) { const LauParameter& par = iter->second; Double_t eventsbyTagCat = par.value() * nEvts; Double_t eventsB0byTagCat = TMath::Nint(eventsbyTagCat/2.0 * (1.0 - sigAsym)); Double_t eventsB0barbyTagCat = TMath::Nint(eventsbyTagCat/2.0 * (1.0 + sigAsym)); if (this->doPoissonSmearing()) { eventsB0byTagCat = LauRandom::randomFun()->Poisson(eventsB0byTagCat); eventsB0barbyTagCat = LauRandom::randomFun()->Poisson(eventsB0barbyTagCat); } eventsB0[iter->first] = std::make_pair( TMath::Nint(eventsB0byTagCat), evtWeight ); eventsB0bar[iter->first] = std::make_pair( TMath::Nint(eventsB0barbyTagCat), evtWeight ); } // CONVENTION WARNING nEvtsGen[std::make_pair("signal",-1)] = eventsB0; nEvtsGen[std::make_pair("signal",+1)] = eventsB0bar; } else { Double_t rateB0bar = sigModelB0bar_->getDPRate().value(); Double_t rateB0 = sigModelB0_->getDPRate().value(); if ( rateB0bar+rateB0 > 1e-30) { sigAsym = (rateB0bar-rateB0)/(rateB0bar+rateB0); } for (LauTagCatParamMap::const_iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) { const LauParameter& par = iter->second; Double_t eventsbyTagCat = par.value() * nEvts; if (this->doPoissonSmearing()) { eventsbyTagCat = LauRandom::randomFun()->Poisson(eventsbyTagCat); } eventsB0[iter->first] = std::make_pair( TMath::Nint(eventsbyTagCat), evtWeight ); } nEvtsGen[std::make_pair("signal",0)] = eventsB0; // generate signal event, decide tag later. } std::cout<<"INFO in LauTimeDepFitModel::eventsToGenerate : Generating toy MC with:"<first.first); + + // Type + const TString& type(iter->first.first); // find the true flavour curEvtTagFlv_ = iter->first.second; // loop through each tagging category const LauTagCatGenInfo& eventsByTagCat = iter->second; for (LauTagCatGenInfo::const_iterator tagCatIter = eventsByTagCat.begin(); tagCatIter != eventsByTagCat.end(); ++tagCatIter) { // get the tagging category curEvtTagCat_ = tagCatIter->first; // get the event weight for this category const Double_t evtWeight( tagCatIter->second.second ); // get the number of events to generate in this configuration const Int_t nEvtsGen( tagCatIter->second.first ); for (Int_t iEvt(0); iEvtsetGenNtupleDoubleBranchValue( "evtWeight", evtWeight ); if (evtCategory == "signal") { this->setGenNtupleIntegerBranchValue("genSig",1); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], 0 ); + } // All the generate*Event() methods have to fill in curEvtDecayTime_ and curEvtDecayTimeErr_ // In addition, generateSignalEvent has to decide on the tag and fill in curEvtTagFlv_ genOK = this->generateSignalEvent(); } else { - genOK = kFALSE; + this->setGenNtupleIntegerBranchValue("genSig",0); + UInt_t bkgndID(0); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + Int_t gen(0); + if ( bkgndClassNames[iBkgnd] == type ) { + gen = 1; + bkgndID = iBkgnd; + } + this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], gen ); + } + genOK = this->generateBkgndEvent(bkgndID); } if (!genOK) { // If there was a problem with the generation then break out and return. // The problem model will have adjusted itself so that all should be OK next time. break; } if (this->useDP() == kTRUE) { this->setDPDtBranchValues(); // store DP, decay time and tagging variables in the ntuple } // If the tagging category is 0, we don't know what the flavour is... if (curEvtTagCat_==0){ curEvtTagFlv_=0; } // Store the event's tag and tagging category this->setGenNtupleIntegerBranchValue("cpEigenvalue", cpEigenValue_); this->setGenNtupleIntegerBranchValue("tagCat",curEvtTagCat_); this->setGenNtupleIntegerBranchValue("tagFlv",curEvtTagFlv_); this->setGenNtupleDoubleBranchValue(flavTag_->getMistagVarName(),curEvtMistag_); this->setGenNtupleDoubleBranchValue(flavTag_->getTrueTagVarName(),curEvtTrueTagFlv_); //Store calibration parameters and asymmetries during debugging at least //this->setGenNtupleDoubleBranchValue("AProd", AProd_->unblindValue()); //this->setGenNtupleDoubleBranchValue("ARaw", ARaw_->unblindValue()); // Store the event number (within this experiment) // and then increment it this->setGenNtupleIntegerBranchValue("iEvtWithinExpt",evtNum); ++evtNum; // Write the values into the tree this->fillGenNtupleBranches(); // Print an occasional progress message if (iEvt%1000 == 0) {std::cout<<"INFO in LauTimeDepFitModel::genExpt : Generated event number "<useDP() && genOK) { sigModelB0bar_->checkToyMC(kTRUE); sigModelB0_->checkToyMC(kTRUE); std::cout<<"aSqMaxSet = "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } LauParArray fitFracB0 = sigModelB0_->getFitFractions(); if (fitFracB0.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); igetMeanEff().value()); meanEffB0_.value(sigModelB0_->getMeanEff().value()); DPRateB0bar_.value(sigModelB0bar_->getDPRate().value()); DPRateB0_.value(sigModelB0_->getDPRate().value()); } } // If we're reusing embedded events or if the generation is being // reset then clear the lists of used events //if (!signalTree_.empty() && (reuseSignal_ || !genOK)) { if (reuseSignal_ || !genOK) { for(LauTagCatEmbDataMap::const_iterator iter = signalTree_.begin(); iter != signalTree_.end(); ++iter) { (iter->second)->clearUsedList(); } } + for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { + if (reuseBkgnd_[bkgndID] || !genOK) { + for (LauTagCatEmbDataMapList::iterator iterlist = bkgndTree_.begin(); iterlist != bkgndTree_.end(); ++iterlist){ + for (LauTagCatEmbDataMap::iterator iter = (*iterlist).begin(); iter != (*iterlist).end(); ++iter){ + (iter->second)->clearUsedList(); + } + } + } + } return genOK; } Bool_t LauTimeDepFitModel::generateSignalEvent() { // Generate signal event, including SCF if necessary. // DP:DecayTime generation follows. // If it's ok, we then generate mES, DeltaE, Fisher/NN... Bool_t genOK(kTRUE); Bool_t generatedEvent(kFALSE); Bool_t doSquareDP = kinematicsB0bar_->squareDP(); doSquareDP &= kinematicsB0_->squareDP(); LauKinematics* kinematics(kinematicsB0bar_); // find the right decay time PDF for the current tagging category LauTagCatDtPdfMap::const_iterator dt_iter = signalDecayTimePdfs_.find(curEvtTagCat_); LauDecayTimePdf* decayTimePdf = (dt_iter != signalDecayTimePdfs_.end()) ? dt_iter->second : 0; // add flavour tagging for time pdf decayTimePdf->setFlavourTagging(flavTag_); // find the right embedded data for the current tagging category LauTagCatEmbDataMap::const_iterator emb_iter = signalTree_.find(curEvtTagCat_); LauEmbeddedData* embeddedData = (emb_iter != signalTree_.end()) ? emb_iter->second : 0; // find the right extra PDFs for the current tagging category LauTagCatPdfListMap::iterator extra_iter = sigExtraPdf_.find(curEvtTagCat_); LauPdfList* extraPdfs = (extra_iter != sigExtraPdf_.end()) ? &(extra_iter->second) : 0; if (this->useDP()) { if (embeddedData) { embeddedData->getEmbeddedEvent(kinematics); curEvtTagFlv_ = TMath::Nint(embeddedData->getValue("tagFlv")); curEvtDecayTimeErr_ = embeddedData->getValue(decayTimePdf->varErrName()); curEvtDecayTime_ = embeddedData->getValue(decayTimePdf->varName()); if (embeddedData->haveBranch("mcMatch")) { Int_t match = TMath::Nint(embeddedData->getValue("mcMatch")); if (match) { this->setGenNtupleIntegerBranchValue("genTMSig",1); this->setGenNtupleIntegerBranchValue("genSCFSig",0); } else { this->setGenNtupleIntegerBranchValue("genTMSig",0); this->setGenNtupleIntegerBranchValue("genSCFSig",1); } } } else { nGenLoop_ = 0; // generate the decay time error (NB the kTRUE forces the generation of a new value) curEvtDecayTimeErr_ = decayTimePdf->generateError(kTRUE); while (generatedEvent == kFALSE && nGenLoop_ < iterationsMax_) { // Calculate the unnormalised truth-matched signal likelihood // First let define the tag flavour CONVENTION WARNING Double_t randNo = LauRandom::randomFun()->Rndm(); if (randNo < 0.5) { curEvtTagFlv_ = +1; // B0bar tag } else { curEvtTagFlv_ = -1; // B0 tag } curEvtTrueTagFlv_ = 0; if (decayTimePdf->getFuncType() == LauDecayTimePdf::SimFitNormBd || decayTimePdf->getFuncType() == LauDecayTimePdf::SimFitNormBs || curEvtTagCat_==0){ // First let define the tag flavour CONVENTION WARNING Double_t rand = LauRandom::randomFun()->Rndm(); if (rand < 0.5) { curEvtTrueTagFlv_ = +1; // B0bar tag } else { curEvtTrueTagFlv_ = -1; // B0 tag } } // Calculate event quantities that depend only on the tagCat and tagFlv Double_t qD(0.); Double_t qDDo2(0.); if(flavTag_->getUsePerEvtMistag() && curEvtTagCat_!=0 ){ curEvtMistag_ = flavTag_->getEtaGen(sigFlavTagPdf_[curEvtTagCat_]); Double_t omega = flavTag_->getOmegaGen(curEvtMistag_,curEvtTagFlv_,curEvtTagCat_); qD = curEvtTagFlv_*(1.0-2.0*omega); //qD = curEvtTagFlv_*(1-2*curEvtMistag_); }else{ curEvtMistag_ = 0.5; LauTagCatParamMap dilution_ = flavTag_->getDilution(); LauTagCatParamMap deltaDilution_ = flavTag_->getDeltaDilution(); qD = curEvtTagFlv_*dilution_[curEvtTagCat_].unblindValue(); qDDo2 = curEvtTagFlv_*0.5*deltaDilution_[curEvtTagCat_].unblindValue(); } // Generate the DP position Double_t m13Sq(0.0), m23Sq(0.0); kinematicsB0bar_->genFlatPhaseSpace(m13Sq, m23Sq); // Next, calculate the total A and Abar for the given DP position sigModelB0_->calcLikelihoodInfo(m13Sq, m23Sq); sigModelB0bar_->calcLikelihoodInfo(m13Sq, m23Sq); // Retrieve the amplitudes and efficiency from the dynamics const LauComplex& Abar = sigModelB0bar_->getEvtDPAmp(); const LauComplex& A = sigModelB0_->getEvtDPAmp(); Double_t eff = sigModelB0bar_->getEvtEff(); // Next calculate the DP terms Double_t aSqSum = A.abs2() + Abar.abs2(); Double_t aSqDif = A.abs2() - Abar.abs2(); LauComplex inter = Abar * A.conj() * phiMixComplex_; Double_t interTermIm = 2.0 * inter.im(); Double_t interTermRe = 2.0 * inter.re(); // Generate decay time const Double_t tMin = decayTimePdf->minAbscissa(); const Double_t tMax = decayTimePdf->maxAbscissa(); curEvtDecayTime_ = LauRandom::randomFun()->Rndm()*(tMax-tMin) + tMin; // Calculate all the decay time info //decayTimePdf->calcLikelihoodInfo(curEvtDecayTime_,curEvtDecayTimeErr_); decayTimePdf->calcLikelihoodInfoGen(curEvtDecayTime_,curEvtDecayTimeErr_,curEvtTagFlv_,curEvtMistag_,curEvtTagCat_,curEvtTrueTagFlv_); // ...and check that the calculation went ok, otherwise loop again if (decayTimePdf->state() != LauDecayTimePdf::Good) { std::cout<<"decayTimePdf state is bad"<getEffiTerm(); // First get all the decay time terms Double_t dtCos = decayTimePdf->getCosTerm(); Double_t dtSin = decayTimePdf->getSinTerm(); Double_t dtCosh = decayTimePdf->getCoshTerm(); Double_t dtSinh = decayTimePdf->getSinhTerm(); // Combine all terms Double_t cosTerm = dtCos * qD * aSqDif; Double_t sinTerm = dtSin * qD * interTermIm; Double_t coshTerm = dtCosh * (1.0 + qDDo2) * aSqSum; Double_t sinhTerm = dtSinh * (1.0 + qDDo2) * interTermRe; if ( cpEigenValue_ == CPOdd ) { sinTerm *= -1.0; sinhTerm *= -1.0; } // ... to get the total and multiply by the efficiency Double_t ASq = coshTerm + cosTerm - sinTerm + sinhTerm; //if (curEvtTrueTagFlv_ !=0 && curEvtTagCat_>0 && curEvtTagFlv_==-1){ // std::cout<<"##### Event details #####"<getNormTerm(); ASq *= eff; ASq *= dtEff; //std::cout << "Total Amplitude Eff: " << ASq << std::endl; //Finally we throw the dice to see whether this event should be generated //We make a distinction between the likelihood of TM and SCF to tag the SCF events as such randNo = LauRandom::randomFun()->Rndm(); if (randNo <= ASq/aSqMaxSet_ ) { generatedEvent = kTRUE; nGenLoop_ = 0; if (ASq > aSqMaxVar_) {aSqMaxVar_ = ASq;} } else { nGenLoop_++; } } // end of while !generatedEvent loop } // end of if (embeddedData) else control } else { if ( embeddedData ) { embeddedData->getEmbeddedEvent(0); curEvtTagFlv_ = TMath::Nint(embeddedData->getValue("tagFlv")); curEvtDecayTimeErr_ = embeddedData->getValue(decayTimePdf->varErrName()); curEvtDecayTime_ = embeddedData->getValue(decayTimePdf->varName()); } } // Check whether we have generated the toy MC OK. if (nGenLoop_ >= iterationsMax_) { aSqMaxSet_ = 1.01 * aSqMaxVar_; genOK = kFALSE; std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Hit max iterations: setting aSqMaxSet_ to "< aSqMaxSet_) { aSqMaxSet_ = 1.01 * aSqMaxVar_; genOK = kFALSE; std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Found a larger ASq value: setting aSqMaxSet_ to "<updateKinematics(kinematicsB0bar_->getm13Sq(), kinematicsB0bar_->getm23Sq() ); this->generateExtraPdfValues(extraPdfs, embeddedData); } // Check for problems with the embedding if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) { std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events."<clearUsedList(); } return genOK; } +Bool_t LauTimeDepFitModel::generateBkgndEvent(UInt_t bkgndID) +{ + // Generate Bkgnd event + Bool_t genOK(kTRUE); + + LauAbsBkgndDPModel* model(0); + LauEmbeddedData* embeddedData(0); + LauPdfList* extraPdfs(0); + LauKinematics* kinematics(0); + + model = BkgndDPModels_[bkgndID]; + if (this->enableEmbedding()) { + // find the right embedded data for the current tagging category + LauTagCatEmbDataMap::const_iterator emb_iter = bkgndTree_[bkgndID].find(curEvtTagCat_); + embeddedData = (emb_iter != bkgndTree_[bkgndID].end()) ? emb_iter->second : 0; + } + extraPdfs = &BkgndPdfs_[bkgndID]; + kinematics = kinematicsB0bar_; + + if (this->useDP()) { + if (embeddedData) { + embeddedData->getEmbeddedEvent(kinematics); + } else { + if (model == 0) { + const TString& bkgndClass = this->bkgndClassName(bkgndID); + std::cerr << "ERROR in LauCPFitModel::generateBkgndEvent : Can't find the DP model for background class \"" << bkgndClass << "\"." << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + genOK = model->generate(); + } + } else { + if (embeddedData) { + embeddedData->getEmbeddedEvent(0); + } + } + if (genOK) { + this->generateExtraPdfValues(extraPdfs, embeddedData); + } + + // Check for problems with the embedding + if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) { + const TString& bkgndClass = this->bkgndClassName(bkgndID); + std::cerr << "WARNING in LauCPFitModel::generateBkgndEvent : Source of embedded " << bkgndClass << " events used up, clearing the list of used events." << std::endl; + embeddedData->clearUsedList(); + } + + return genOK; +} + void LauTimeDepFitModel::setupGenNtupleBranches() { // Setup the required ntuple branches this->addGenNtupleDoubleBranch("evtWeight"); this->addGenNtupleIntegerBranch("genSig"); this->addGenNtupleIntegerBranch("cpEigenvalue"); this->addGenNtupleIntegerBranch("tagFlv"); this->addGenNtupleIntegerBranch("tagCat"); if (this->useDP() == kTRUE) { // Let's add the decay time variables. if (signalDecayTimePdfs_.begin() != signalDecayTimePdfs_.end()) { LauDecayTimePdf* pdf = signalDecayTimePdfs_.begin()->second; this->addGenNtupleDoubleBranch(pdf->varName()); this->addGenNtupleDoubleBranch(pdf->varErrName()); } this->addGenNtupleDoubleBranch("m12"); this->addGenNtupleDoubleBranch("m23"); this->addGenNtupleDoubleBranch("m13"); this->addGenNtupleDoubleBranch("m12Sq"); this->addGenNtupleDoubleBranch("m23Sq"); this->addGenNtupleDoubleBranch("m13Sq"); this->addGenNtupleDoubleBranch("cosHel12"); this->addGenNtupleDoubleBranch("cosHel23"); this->addGenNtupleDoubleBranch("cosHel13"); if (kinematicsB0bar_->squareDP() && kinematicsB0_->squareDP()) { this->addGenNtupleDoubleBranch("mPrime"); this->addGenNtupleDoubleBranch("thPrime"); } // Can add the real and imaginary parts of the B0 and B0bar total // amplitudes seen in the generation (restrict this with a flag // that defaults to false) if ( storeGenAmpInfo_ ) { this->addGenNtupleDoubleBranch("reB0Amp"); this->addGenNtupleDoubleBranch("imB0Amp"); this->addGenNtupleDoubleBranch("reB0barAmp"); this->addGenNtupleDoubleBranch("imB0barAmp"); } } // Let's look at the extra variables for signal in one of the tagging categories if ( ! sigExtraPdf_.empty() ) { LauPdfList oneTagCatPdfList = sigExtraPdf_.begin()->second; for (LauPdfList::const_iterator pdf_iter = oneTagCatPdfList.begin(); pdf_iter != oneTagCatPdfList.end(); ++pdf_iter) { for ( std::vector::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { this->addGenNtupleDoubleBranch( (*var_iter) ); } } } } } void LauTimeDepFitModel::setDPDtBranchValues() { // Store the decay time variables. if (signalDecayTimePdfs_.begin() != signalDecayTimePdfs_.end()) { LauDecayTimePdf* pdf = signalDecayTimePdfs_.begin()->second; this->setGenNtupleDoubleBranchValue(pdf->varName(),curEvtDecayTime_); this->setGenNtupleDoubleBranchValue(pdf->varErrName(),curEvtDecayTimeErr_); } // CONVENTION WARNING LauKinematics* kinematics(0); if (curEvtTagFlv_<0) { kinematics = kinematicsB0_; } else { kinematics = kinematicsB0bar_; } // Store all the DP information this->setGenNtupleDoubleBranchValue("m12", kinematics->getm12()); this->setGenNtupleDoubleBranchValue("m23", kinematics->getm23()); this->setGenNtupleDoubleBranchValue("m13", kinematics->getm13()); this->setGenNtupleDoubleBranchValue("m12Sq", kinematics->getm12Sq()); this->setGenNtupleDoubleBranchValue("m23Sq", kinematics->getm23Sq()); this->setGenNtupleDoubleBranchValue("m13Sq", kinematics->getm13Sq()); this->setGenNtupleDoubleBranchValue("cosHel12", kinematics->getc12()); this->setGenNtupleDoubleBranchValue("cosHel23", kinematics->getc23()); this->setGenNtupleDoubleBranchValue("cosHel13", kinematics->getc13()); if (kinematics->squareDP()) { this->setGenNtupleDoubleBranchValue("mPrime", kinematics->getmPrime()); this->setGenNtupleDoubleBranchValue("thPrime", kinematics->getThetaPrime()); } // Can add the real and imaginary parts of the B0 and B0bar total // amplitudes seen in the generation (restrict this with a flag // that defaults to false) if ( storeGenAmpInfo_ ) { if ( this->getGenNtupleIntegerBranchValue("genSig")==1 ) { LauComplex Abar = sigModelB0bar_->getEvtDPAmp(); LauComplex A = sigModelB0_->getEvtDPAmp(); this->setGenNtupleDoubleBranchValue("reB0Amp", A.re()); this->setGenNtupleDoubleBranchValue("imB0Amp", A.im()); this->setGenNtupleDoubleBranchValue("reB0barAmp", Abar.re()); this->setGenNtupleDoubleBranchValue("imB0barAmp", Abar.im()); } else { this->setGenNtupleDoubleBranchValue("reB0Amp", 0.0); this->setGenNtupleDoubleBranchValue("imB0Amp", 0.0); this->setGenNtupleDoubleBranchValue("reB0barAmp", 0.0); this->setGenNtupleDoubleBranchValue("imB0barAmp", 0.0); } } } void LauTimeDepFitModel::generateExtraPdfValues(LauPdfList* extraPdfs, LauEmbeddedData* embeddedData) { // CONVENTION WARNING LauKinematics* kinematics(0); if (curEvtTagFlv_<0) { kinematics = kinematicsB0_; } else { kinematics = kinematicsB0bar_; } // Generate from the extra PDFs if (extraPdfs) { for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) { LauFitData genValues; if (embeddedData) { genValues = embeddedData->getValues( (*pdf_iter)->varNames() ); } else { genValues = (*pdf_iter)->generate(kinematics); } for ( LauFitData::const_iterator var_iter = genValues.begin(); var_iter != genValues.end(); ++var_iter ) { TString varName = var_iter->first; if ( varName != "m13Sq" && varName != "m23Sq" ) { Double_t value = var_iter->second; this->setGenNtupleDoubleBranchValue(varName,value); } } } } } void LauTimeDepFitModel::propagateParUpdates() { // Update the complex mixing phase if (useSinCos_) { phiMixComplex_.setRealPart(cosPhiMix_.unblindValue()); phiMixComplex_.setImagPart(-1.0*sinPhiMix_.unblindValue()); } else { phiMixComplex_.setRealPart(TMath::Cos(-1.0*phiMix_.unblindValue())); phiMixComplex_.setImagPart(TMath::Sin(-1.0*phiMix_.unblindValue())); } // Update the total normalisation for the signal likelihood if (this->useDP() == kTRUE) { this->updateCoeffs(); sigModelB0bar_->updateCoeffs(coeffsB0bar_); sigModelB0_->updateCoeffs(coeffsB0_); this->calcInterTermNorm(); } // Update the signal events from the background numbers if not doing an extended fit // And update the tagging category fractions this->updateSigEvents(); } void LauTimeDepFitModel::updateSigEvents() { // The background parameters will have been set from Minuit. // We need to update the signal events using these. if (!this->doEMLFit()) { Double_t nTotEvts = this->eventsPerExpt(); Double_t signalEvents = nTotEvts; - // TODO loop over the background yields and subtract signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts); + for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { + LauAbsRValue* nBkgndEvents = (*iter); + if ( nBkgndEvents->isLValue() ) { + LauParameter* yield = dynamic_cast( nBkgndEvents ); + yield->range(-2.0*nTotEvts,2.0*nTotEvts); + } + } + + // Subtract background events (if any) from signal. + if (usingBkgnd_ == kTRUE) { + for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { + signalEvents -= (*iter)->value(); + } + } if ( ! signalEvents_->fixed() ) { signalEvents_->value(signalEvents); } } // tagging-category fractions for signal events flavTag_->setFirstTagCatFractions(); } void LauTimeDepFitModel::cacheInputFitVars() { // Fill the internal data trees of the signal and background models. // Note that we store the events of both charges in both the // negative and the positive models. It's only later, at the stage // when the likelihood is being calculated, that we separate them. LauFitDataTree* inputFitData = this->fitData(); evtCPEigenVals_.clear(); const Bool_t hasCPEV = ( (cpevVarName_ != "") && inputFitData->haveBranch( cpevVarName_ ) ); UInt_t nEvents = inputFitData->nEvents(); evtCPEigenVals_.reserve( nEvents ); LauFitData::const_iterator fitdata_iter; for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) { const LauFitData& dataValues = inputFitData->getData(iEvt); // if the CP-eigenvalue is in the data use those, otherwise use the default if ( hasCPEV ) { fitdata_iter = dataValues.find( cpevVarName_ ); const Int_t cpEV = static_cast( fitdata_iter->second ); if ( cpEV == 1 ) { cpEigenValue_ = CPEven; } else if ( cpEV == -1 ) { cpEigenValue_ = CPOdd; } else { std::cerr<<"WARNING in LauTimeDepFitModel::cacheInputFitVars : Unknown value: "<cacheInputFitVars(inputFitData); if (this->useDP() == kTRUE) { // DecayTime and SigmaDecayTime for (LauTagCatDtPdfMap::iterator dt_iter = signalDecayTimePdfs_.begin(); dt_iter != signalDecayTimePdfs_.end(); ++dt_iter) { (*dt_iter).second->setFlavourTagging(flavTag_); (*dt_iter).second->cacheInfo(*inputFitData); } } // ...and then the extra PDFs for (LauTagCatPdfListMap::iterator pdf_iter = sigExtraPdf_.begin(); pdf_iter != sigExtraPdf_.end(); ++pdf_iter) { this->cacheInfo(pdf_iter->second, *inputFitData); } + if(usingBkgnd_ == kTRUE){ + for (LauBkgndPdfsList::iterator iter = BkgndPdfs_.begin(); iter != BkgndPdfs_.end(); ++iter) { + this->cacheInfo((*iter), *inputFitData); + } + } if (this->useDP() == kTRUE) { sigModelB0bar_->fillDataTree(*inputFitData); sigModelB0_->fillDataTree(*inputFitData); + if (usingBkgnd_ == kTRUE) { + for (LauBkgndDPModelList::iterator iter = BkgndDPModels_.begin(); iter != BkgndDPModels_.end(); ++iter) { + (*iter)->fillDataTree(*inputFitData); + } + } } } Double_t LauTimeDepFitModel::getTotEvtLikelihood(const UInt_t iEvt) { // Find out whether the tag-side B was a B0 or a B0bar. curEvtTagFlv_ = flavTag_->getEvtTagFlvVals(iEvt); // Also get the tagging category. curEvtTagCat_ = flavTag_->getEvtTagCatVals(iEvt); // Get the CP eigenvalue of the current event cpEigenValue_ = evtCPEigenVals_[iEvt]; // Get the DP and DecayTime likelihood for signal (TODO and eventually backgrounds) this->getEvtDPDtLikelihood(iEvt); // Get the flavour tagging likelihood from eta PDFs (per tagging category - TODO backgrounds to come later) sigFlavTagLike_ = 1.0; //this->getEvtFlavTagLikelihood(iEvt); // Get the combined extra PDFs likelihood for signal (TODO and eventually backgrounds) this->getEvtExtraLikelihoods(iEvt); // Construct the total likelihood for signal, qqbar and BBbar backgrounds Double_t sigLike = sigDPLike_ * sigFlavTagLike_ * sigExtraLike_; //std::cout << "DP like = " << sigDPLike_ << std::endl; //std::cout << "flav tag like = " << sigFlavTagLike_ << std::endl; //std::cout << "extra like = " << sigExtraLike_ << std::endl; Double_t signalEvents = signalEvents_->unblindValue(); if (this->useDP() == kFALSE) { signalEvents *= 0.5 * (1.0 + curEvtTagFlv_ * signalAsym_->unblindValue()); } // TODO better to store this info for each event LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac(); const Double_t sigTagCatFrac = signalTagCatFrac[curEvtTagCat_].unblindValue(); //std::cout << "Signal tag cat frac = " << sigTagCatFrac << std::endl; // Construct the total event likelihood Double_t likelihood(sigLike*sigTagCatFrac); //std::cout << "Likelihoood = " << likelihood << std::endl; if ( ! signalEvents_->fixed() ) { likelihood *= signalEvents; } return likelihood; } Double_t LauTimeDepFitModel::getEventSum() const { Double_t eventSum(0.0); eventSum += signalEvents_->unblindValue(); return eventSum; } void LauTimeDepFitModel::getEvtDPDtLikelihood(const UInt_t iEvt) { // Function to return the signal and background likelihoods for the // Dalitz plot for the given event evtNo. - sigDPLike_ = 1.0; //There's always a likelihood term for signal, so we better not zero it. - if ( this->useDP() == kFALSE ) { + if ( ! this->useDP() ) { + // There's always going to be a term in the likelihood for the + // signal, so we'd better not zero it. + sigDPLike_ = 1.0; + + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { + if (usingBkgnd_ == kTRUE) { + bkgndDPLike_[bkgndID] = 1.0; + } else { + bkgndDPLike_[bkgndID] = 0.0; + } + } + return; } // Mistag probabilities. Defined as: omega = prob of the tagging B0 being reported as B0bar // Whether we want omega or omegaBar depends on q_tag, hence curEvtTagFlv_*... in the previous lines //Double_t misTagFrac = 0.5 * (1.0 - dilution_[curEvtTagCat_] - qDDo2); //Double_t misTagFracBar = 0.5 * (1.0 - dilution_[curEvtTagCat_] + qDDo2); // Calculate event quantities Double_t qD(0); Double_t qDDo2(0); if (flavTag_->getUsePerEvtMistag()){ //qDDo2 term accounted for automatically with per event information Double_t omega = flavTag_->getOmega(iEvt); qD = curEvtTagFlv_*(1.0-2.0*omega); } else { // TODO need to sort out references here LauTagCatParamMap dilution_ = flavTag_->getDilution(); LauTagCatParamMap deltaDilution_ = flavTag_->getDeltaDilution(); qD = curEvtTagFlv_*dilution_[curEvtTagCat_].unblindValue(); qDDo2 = curEvtTagFlv_*0.5*deltaDilution_[curEvtTagCat_].unblindValue(); } // Get the dynamics to calculate everything required for the likelihood calculation sigModelB0bar_->calcLikelihoodInfo(iEvt); sigModelB0_->calcLikelihoodInfo(iEvt); + // Background part + // TODO add them into the actual Likelihood calculations + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { + if (usingBkgnd_ == kTRUE) { + bkgndDPLike_[bkgndID] = BkgndDPModels_[bkgndID]->getLikelihood(iEvt); + } else { + bkgndDPLike_[bkgndID] = 0.0; + } + } + // Retrieve the amplitudes and efficiency from the dynamics const LauComplex& Abar = sigModelB0bar_->getEvtDPAmp(); const LauComplex& A = sigModelB0_->getEvtDPAmp(); Double_t eff = sigModelB0bar_->getEvtEff(); // Next calculate the DP terms Double_t aSqSum = A.abs2() + Abar.abs2(); Double_t aSqDif = A.abs2() - Abar.abs2(); LauComplex inter = Abar * A.conj() * phiMixComplex_; Double_t interTermIm = 2.0 * inter.im(); Double_t interTermRe = 2.0 * inter.re(); // First get all the decay time terms LauDecayTimePdf* decayTimePdf = signalDecayTimePdfs_[curEvtTagCat_]; decayTimePdf->calcLikelihoodInfo(iEvt); + // TODO Backgrounds + // Get the decay time acceptance Double_t dtEff = decayTimePdf->getEffiTerm(); Double_t dtEffN = decayTimePdf->getEffiTermNorm(); // First get all the decay time terms Double_t dtCos = decayTimePdf->getCosTerm(); Double_t dtSin = decayTimePdf->getSinTerm(); Double_t dtCosh = decayTimePdf->getCoshTerm(); Double_t dtSinh = decayTimePdf->getSinhTerm(); //std::cout << dtCos << " " << dtSin << " " << dtCosh << " " << dtSinh << std::endl; Double_t cosTerm = dtCos * qD * aSqDif; Double_t sinTerm = dtSin * qD * interTermIm; Double_t coshTerm = dtCosh * (1.0 + qDDo2) * aSqSum; Double_t sinhTerm = dtSinh * (1.0 + qDDo2) * interTermRe; if ( cpEigenValue_ == CPOdd ) { sinTerm *= -1.0; sinhTerm *= -1.0; } // ... to get the total and multiply by the efficiency (DP and decay time) Double_t ASq = coshTerm + cosTerm - sinTerm + sinhTerm; ASq *= eff; ASq *= dtEff; //std::cout<<"Cosh Cos Sin Sinh "<< coshTerm<< " " << cosTerm << " " << sinTerm << " " << sinhTerm << std::endl; //std::cout << "Total Amplitude : " << ASq << std::endl; // Calculate the DP and time normalisation Double_t normTermIndep = sigModelB0bar_->getDPNorm() + sigModelB0_->getDPNorm(); Double_t normTermCosh(0.0); Double_t norm(0.0); if (decayTimePdf->getFuncType() == LauDecayTimePdf::ExpTrig || decayTimePdf->getFuncType() == LauDecayTimePdf::SimFitNormBd || decayTimePdf->getFuncType() == LauDecayTimePdf::SimFitSigBd){ normTermCosh = decayTimePdf->getNormTermExp(); norm = normTermIndep*normTermCosh; //Test norm *= dtEffN; } if (decayTimePdf->getFuncType() == LauDecayTimePdf::ExpHypTrig){ normTermCosh = decayTimePdf->getNormTermCosh(); Double_t normTermDep = interTermReNorm_; Double_t normTermSinh = decayTimePdf->getNormTermSinh(); norm = normTermIndep*normTermCosh + normTermDep*normTermSinh; } // Calculate the normalised signal likelihood //std::cout << "ASq = " << ASq << std::endl; //std::cout << "norm = " << norm << std::endl; sigDPLike_ = ASq / norm; } void LauTimeDepFitModel::getEvtExtraLikelihoods(const UInt_t iEvt) { // Function to return the signal and background likelihoods for the // extra variables for the given event evtNo. sigExtraLike_ = 1.0; //There's always a likelihood term for signal, so we better not zero it. // First, those independent of the tagging of the event: // signal LauTagCatPdfListMap::iterator sig_iter = sigExtraPdf_.find(curEvtTagCat_); LauPdfList* pdfList = (sig_iter != sigExtraPdf_.end())? &(sig_iter->second) : 0; if (pdfList) { sigExtraLike_ = this->prodPdfValue( *pdfList, iEvt ); } + // Background + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { + if (usingBkgnd_) { + bkgndExtraLike_[bkgndID] = this->prodPdfValue( BkgndPdfs_[bkgndID], iEvt ); + } else { + bkgndExtraLike_[bkgndID] = 0.0; + } + } + } void LauTimeDepFitModel::getEvtFlavTagLikelihood(const UInt_t iEvt) { // Function to return the signal and background likelihoods for the // extra variables for the given event evtNo. sigFlavTagLike_ = 1.0; //There's always a likelihood term for signal, so we better not zero it. // First, those independent of the tagging of the event: // signal LauAbsPdf* pdf = sigFlavTagPdf_[curEvtTagCat_]; if (pdf) { pdf->calcLikelihoodInfo(iEvt); sigFlavTagLike_ = pdf->getLikelihood(); } if (sigFlavTagLike_<=0){ std::cout<<"INFO in LauTimeDepFitModel::getEvtFlavTagLikelihood : Event with 0 FlavTag Liklihood"<antiparticleCoeff()); coeffsB0_.push_back(coeffPars_[i]->particleCoeff()); } } void LauTimeDepFitModel::checkMixingPhase() { Double_t phase = phiMix_.value(); Double_t genPhase = phiMix_.genValue(); // Check now whether the phase lies in the right range (-pi to pi). Bool_t withinRange(kFALSE); while (withinRange == kFALSE) { if (phase > -LauConstants::pi && phase < LauConstants::pi) { withinRange = kTRUE; } else { // Not within the specified range if (phase > LauConstants::pi) { phase -= LauConstants::twoPi; } else if (phase < -LauConstants::pi) { phase += LauConstants::twoPi; } } } // A further problem can occur when the generated phase is close to -pi or pi. // The phase can wrap over to the other end of the scale - // this leads to artificially large pulls so we wrap it back. Double_t diff = phase - genPhase; if (diff > LauConstants::pi) { phase -= LauConstants::twoPi; } else if (diff < -LauConstants::pi) { phase += LauConstants::twoPi; } // finally store the new value in the parameter // and update the pull phiMix_.value(phase); phiMix_.updatePull(); } void LauTimeDepFitModel::embedSignal(Int_t tagCat, const TString& fileName, const TString& treeName, Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment) { if (signalTree_[tagCat]) { std::cerr<<"ERROR in LauTimeDepFitModel::embedSignal : Already embedding signal from file for tagging category "<findBranches(); if (!dataOK) { delete signalTree_[tagCat]; signalTree_[tagCat] = 0; std::cerr<<"ERROR in LauTimeDepFitModel::embedSignal : Problem creating data tree for embedding."<validBkgndClass( bkgndClass ) ) { + std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl; + std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; + return; + } + + UInt_t bkgndID = this->bkgndClassID( bkgndClass ); + + LauTagCatEmbDataMap bkgTree = bkgndTree_[bkgndID]; + + if (bkgTree[tagCat]) { + std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Already embedding background from a file." << std::endl; + return; + } + + bkgTree[tagCat] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment); + + Bool_t dataOK = bkgTree[tagCat]->findBranches(); + if (!dataOK) { + delete bkgTree[tagCat]; bkgTree[tagCat] = 0; + std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Problem creating data tree for embedding." << std::endl; + return; + } + + reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble; + + if (this->enableEmbedding() == kFALSE) { + this->enableEmbedding(kTRUE); + } +} + void LauTimeDepFitModel::setupSPlotNtupleBranches() { // add branches for storing the experiment number and the number of // the event within the current experiment this->addSPlotNtupleIntegerBranch("iExpt"); this->addSPlotNtupleIntegerBranch("iEvtWithinExpt"); // Store the efficiency of the event (for inclusive BF calculations). if (this->storeDPEff()) { this->addSPlotNtupleDoubleBranch("efficiency"); } // Store the total event likelihood for each species. this->addSPlotNtupleDoubleBranch("sigTotalLike"); + + if (usingBkgnd_) { + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + TString name( this->bkgndClassName(iBkgnd) ); + name += "TotalLike"; + this->addSPlotNtupleDoubleBranch(name); + } + } // Store the DP likelihoods if (this->useDP()) { this->addSPlotNtupleDoubleBranch("sigDPLike"); + if (usingBkgnd_) { + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + TString name( this->bkgndClassName(iBkgnd) ); + name += "DPLike"; + this->addSPlotNtupleDoubleBranch(name); + } + } } // Store the likelihoods for each extra PDF const LauPdfList* pdfList( &(sigExtraPdf_.begin()->second) ); this->addSPlotNtupleBranches(pdfList, "sig"); + if (usingBkgnd_) { + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + const TString& bkgndClass = this->bkgndClassName(iBkgnd); + const LauPdfList* pdfList2 = &(BkgndPdfs_[iBkgnd]); + this->addSPlotNtupleBranches(pdfList2, bkgndClass); + } + } } void LauTimeDepFitModel::addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix) { if (!extraPdfs) { return; } // Loop through each of the PDFs for (LauPdfList::const_iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) { // Count the number of input variables that are not // DP variables (used in the case where there is DP // dependence for e.g. MVA) UInt_t nVars(0); for ( std::vector::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nVars; } } if ( nVars == 1 ) { // If the PDF only has one variable then // simply add one branch for that variable TString varName = (*pdf_iter)->varName(); TString name(prefix); name += varName; name += "Like"; this->addSPlotNtupleDoubleBranch(name); } else if ( nVars == 2 ) { // If the PDF has two variables then we // need a branch for them both together and // branches for each TString allVars(""); for ( std::vector::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) { allVars += (*var_iter); TString name(prefix); name += (*var_iter); name += "Like"; this->addSPlotNtupleDoubleBranch(name); } TString name(prefix); name += allVars; name += "Like"; this->addSPlotNtupleDoubleBranch(name); } else { std::cerr<<"WARNING in LauTimeDepFitModel::addSPlotNtupleBranches : Can't yet deal with 3D PDFs."<begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) { // calculate the likelihood for this event (*pdf_iter)->calcLikelihoodInfo(iEvt); extraLike = (*pdf_iter)->getLikelihood(); totalLike *= extraLike; // Count the number of input variables that are not // DP variables (used in the case where there is DP // dependence for e.g. MVA) UInt_t nVars(0); for ( std::vector::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nVars; } } if ( nVars == 1 ) { // If the PDF only has one variable then // simply store the value for that variable TString varName = (*pdf_iter)->varName(); TString name(prefix); name += varName; name += "Like"; this->setSPlotNtupleDoubleBranchValue(name, extraLike); } else if ( nVars == 2 ) { // If the PDF has two variables then we // store the value for them both together // and for each on their own TString allVars(""); for ( std::vector::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) { allVars += (*var_iter); TString name(prefix); name += (*var_iter); name += "Like"; Double_t indivLike = (*pdf_iter)->getLikelihood( (*var_iter) ); this->setSPlotNtupleDoubleBranchValue(name, indivLike); } TString name(prefix); name += allVars; name += "Like"; this->setSPlotNtupleDoubleBranchValue(name, extraLike); } else { std::cerr<<"WARNING in LauAllFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs."<useDP()) { nameSet.insert("DP"); } LauPdfList pdfList( (sigExtraPdf_.begin()->second) ); for (LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) { // Loop over the variables involved in each PDF for ( std::vector::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) { // If they are not DP coordinates then add them if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { nameSet.insert( (*var_iter) ); } } } return nameSet; } LauSPlot::NumbMap LauTimeDepFitModel::freeSpeciesNames() const { LauSPlot::NumbMap numbMap; if (!signalEvents_->fixed() && this->doEMLFit()) { numbMap["sig"] = signalEvents_->genValue(); } + if ( usingBkgnd_ ) { + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + const TString& bkgndClass = this->bkgndClassName(iBkgnd); + const LauAbsRValue* par = bkgndEvents_[iBkgnd]; + if (!par->fixed()) { + numbMap[bkgndClass] = par->genValue(); + if ( ! par->isLValue() ) { + std::cerr << "WARNING in LauTimeDepFitModel::freeSpeciesNames : \"" << par->name() << "\" is a LauFormulaPar, which implies it is perhaps not entirely free to float in the fit, so the sWeight calculation may not be reliable" << std::endl; + } + } + } + } return numbMap; } LauSPlot::NumbMap LauTimeDepFitModel::fixdSpeciesNames() const { LauSPlot::NumbMap numbMap; if (signalEvents_->fixed() && this->doEMLFit()) { numbMap["sig"] = signalEvents_->genValue(); } + if ( usingBkgnd_ ) { + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + const TString& bkgndClass = this->bkgndClassName(iBkgnd); + const LauAbsRValue* par = bkgndEvents_[iBkgnd]; + if (par->fixed()) { + numbMap[bkgndClass] = par->genValue(); + } + } + } return numbMap; } LauSPlot::TwoDMap LauTimeDepFitModel::twodimPDFs() const { LauSPlot::TwoDMap twodimMap; const LauPdfList* pdfList = &(sigExtraPdf_.begin()->second); for (LauPdfList::const_iterator pdf_iter = pdfList->begin(); pdf_iter != pdfList->end(); ++pdf_iter) { // Count the number of input variables that are not DP variables UInt_t nVars(0); for ( std::vector::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nVars; } } if ( nVars == 2 ) { twodimMap.insert( std::make_pair( "sig", std::make_pair( (*pdf_iter)->varNames()[0], (*pdf_iter)->varNames()[1] ) ) ); } } + if (usingBkgnd_) { + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + const TString& bkgndClass = this->bkgndClassName(iBkgnd); + const LauPdfList& pdfList2 = BkgndPdfs_[iBkgnd]; + for (LauPdfList::const_iterator pdf_iter = pdfList2.begin(); pdf_iter != pdfList2.end(); ++pdf_iter) { + // Count the number of input variables that are not DP variables + UInt_t nVars(0); + std::vector varNames = (*pdf_iter)->varNames(); + for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { + if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { + ++nVars; + } + } + if ( nVars == 2 ) { + twodimMap.insert( std::make_pair( bkgndClass, std::make_pair( varNames[0], varNames[1] ) ) ); + } + } + } + } return twodimMap; } void LauTimeDepFitModel::storePerEvtLlhds() { std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Storing per-event likelihood values..."<fitData(); // if we've not been using the DP model then we need to cache all // the info here so that we can get the efficiency from it if (!this->useDP() && this->storeDPEff()) { sigModelB0bar_->initialise(coeffsB0bar_); sigModelB0_->initialise(coeffsB0_); sigModelB0bar_->fillDataTree(*inputFitData); sigModelB0_->fillDataTree(*inputFitData); } UInt_t evtsPerExpt(this->eventsPerExpt()); LauIsobarDynamics* sigModel(sigModelB0bar_); for (UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt) { // Find out whether we have B0bar or B0 curEvtTagFlv_ = flavTag_->getEvtTagFlvVals(iEvt); curEvtTagCat_ = flavTag_->getEvtTagCatVals(iEvt); curEvtMistag_ = flavTag_->getOmega(iEvt); LauTagCatPdfListMap::iterator sig_iter = sigExtraPdf_.find(curEvtTagCat_); LauPdfList* sigPdfs = (sig_iter != sigExtraPdf_.end())? &(sig_iter->second) : 0; // the DP information this->getEvtDPDtLikelihood(iEvt); if (this->storeDPEff()) { if (!this->useDP()) { sigModel->calcLikelihoodInfo(iEvt); } this->setSPlotNtupleDoubleBranchValue("efficiency",sigModel->getEvtEff()); } if (this->useDP()) { sigTotalLike_ = sigDPLike_; this->setSPlotNtupleDoubleBranchValue("sigDPLike",sigDPLike_); + if (usingBkgnd_) { + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + TString name = this->bkgndClassName(iBkgnd); + name += "DPLike"; + this->setSPlotNtupleDoubleBranchValue(name,bkgndDPLike_[iBkgnd]); + } + } } else { sigTotalLike_ = 1.0; } // the signal PDF values sigTotalLike_ *= this->setSPlotNtupleBranchValues(sigPdfs, "sig", iEvt); + // the background PDF values + LauBkgndPdfsList* bkgndPdfs(0); + if (usingBkgnd_) { + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + const TString& bkgndClass = this->bkgndClassName(iBkgnd); + LauPdfList& pdfs = (*bkgndPdfs)[iBkgnd]; + bkgndTotalLike_[iBkgnd] *= this->setSPlotNtupleBranchValues(&(pdfs), bkgndClass, iEvt); + } + } + // the total likelihoods this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_); + if (usingBkgnd_) { + const UInt_t nBkgnds = this->nBkgndClasses(); + for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { + TString name = this->bkgndClassName(iBkgnd); + name += "TotalLike"; + this->setSPlotNtupleDoubleBranchValue(name,bkgndTotalLike_[iBkgnd]); + } + } // fill the tree this->fillSPlotNtupleBranches(); } std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Finished storing per-event likelihood values."<validTagCat(tagCat)) { std::cerr<<"ERROR in LauTimeDepFitModel::setSignalFlavTagPdfs : Tagging category \""<