Page MenuHomeHEPForge
Files F11222124

No OneTemporary
Actions

View Options

diff --git a/inc/LauTimeDepFitModel.hh b/inc/LauTimeDepFitModel.hh
index 266ca97..dbe9c96 100644
--- a/inc/LauTimeDepFitModel.hh
+++ b/inc/LauTimeDepFitModel.hh
@@ -1,742 +1,751 @@
/*
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 <vector>
#include <map>
#include <set>
#include <iostream>
#include "TString.h"
#include "TStopwatch.h"
#include "TSystem.h"
#include "LauAbsFitModel.hh"
#include "LauComplex.hh"
#include "LauConstants.hh"
#include "LauEmbeddedData.hh"
#include "LauFlavTag.hh"
#include "LauParameter.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 */
QFS = 0, /*!< Quasi Flavour Specific 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);
//TODO what is this asymmetry?
//! 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);
//TODO what is this asymmetry?
//! Set the background event yield and asymmetry
/*!
The name of the parameters must be that of the corresponding background category (so that they can be correctly assigned)
\param [in] nBkgndEvents contains the background yield and option to fix it
\param [in] bkgndAsym contains the background asymmetry and option to fix it
*/
virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym);
//! Set the background DP models (null pointer for BbarModel implies same model for both)
/*!
\param [in] bkgndClass the name of the background class
\param [in] BModel the DP model of the background for B (particle) decays
\param [in] BbarModel the DP model of the background for Bbar (anti-particle) decays
*/
void setBkgndDPModels(const TString& bkgndClass, LauAbsBkgndDPModel* BModel, LauAbsBkgndDPModel* BbarModel);
//! 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 cpevVarName 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 signal decay time PDF
/*!
\param [in] pdf the signal decay time PDF
*/
void setSignalDtPdf(LauDecayTimePdf* pdf);
//! Set the decay time PDF for a given background class
/*!
\param [in] bkgndClass the name of the background class
\param [in] pdf the background decay time PDF
*/
void setBkgndDtPdf(const TString& bkgndClass, LauDecayTimePdf* pdf);
//! Set the signal PDF for a given variable
/*!
\param [in] pdf the PDF to be added to the signal model
*/
void setSignalPdfs(LauAbsPdf* pdf);
//! Set the background PDF for a given variable and background class
/*!
\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);
//! Embed full simulation events for the signal, rather than generating toy from the PDFs
/*!
\param [in] fileName the name of the file containing the events
\param [in] treeName the name of the tree
\param [in] reuseEventsWithinEnsemble sample with replacement but only replace events once each experiment has been generated
\param [in] reuseEventsWithinExperiment sample with immediate replacement
\param [in] useReweighting perform an accept/reject routine using the configured signal amplitude model based on the MC-truth DP coordinate
*/
void embedSignal(const TString& fileName, const TString& treeName,
const Bool_t reuseEventsWithinEnsemble, const Bool_t reuseEventsWithinExperiment = kFALSE, const Bool_t useReweighting = kFALSE);
//! Embed full simulation events for the signal, rather than generating toy from the PDFs
/*!
\param [in] bkgndClass the name of the background class
\param [in] fileName the name of the file containing the events
\param [in] treeName the name of the tree
\param [in] reuseEventsWithinEnsemble sample with replacement but only replace events once each experiment has been generated
\param [in] reuseEventsWithinExperiment sample with immediate replacement
\param [in] useReweighting perform an accept/reject routine using the configured signal amplitude model based on the MC-truth DP coordinate
*/
void embedBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName,
const Bool_t reuseEventsWithinEnsemble, const Bool_t reuseEventsWithinExperiment = kFALSE, const Bool_t useReweighting = 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);
//! Setup blinding of the value of the mixing phase
/*!
\param [in] blindingString the string used to seed the blinding transformation for the phiMix parameter
\param [in] width the width of the Gaussian from which the offset should be sampled
\param [in] flipSign activate possible random sign flip
*/
void blindPhiMix(const TString& blindingString, const Double_t width = 1.0, const Bool_t flipSign = kTRUE);
//! Setup blinding of the value of the mixing phase
/*!
\param [in] sinBlindingString the string used to seed the blinding transformation for the sin(phiMix) parameter
\param [in] cosBlindingString the string used to seed the blinding transformation for the cos(phiMix) parameter
\param [in] width the width of the Gaussian from which the offset should be sampled
\param [in] flipSign activate possible random sign flip
*/
void blindPhiMix(const TString& sinBlindingString, const TString& cosBlindingString, const Double_t width = 0.5, const Bool_t flipSign = kTRUE);
+ //! Returns whether the mixing phase parameter is blinded
+ Bool_t phiMixBlinded() const { return phiMix_.blind(); }
+
+ //! Returns whether the sine of the mixing phase parameter is blinded
+ Bool_t sinPhiMixBlinded() const { return sinPhiMix_.blind(); }
+
+ //! Returns whether the cosine of the mixing phase parameter is blinded
+ Bool_t cosPhiMixBlinded() const { return cosPhiMix_.blind(); }
+
//! 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);
//! 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 tagging efficiency parameters
void setTagEffParams();
//! 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
/*!
\param [in] AProd is the production asymmetry definied as (N_BqBar-N_Bq)/(N_BqBar+N_Bq)
\param [in] AProdFix whether the production asymmetry should be fixed in the fit
*/
void setAsymmetries(const Double_t AProd, const Bool_t AProdFix);
//! Set production and detections asymmetries
/*!
\param [in] bkgndClass the name of the background class
\param [in] AProd is the production asymmetry definied as (N_BqBar-N_Bq)/(N_BqBar+N_Bq)
\param [in] AProdFix whether the production asymmetry should be fixed in the fit
*/
void setBkgndAsymmetries(const TString& bkgndClass, const Double_t AProd, const Bool_t AProdFix);
//! 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();
protected:
//! container to hold information on the number of events (and their weight) to generate for each category
typedef std::map<TString,std::pair<Int_t,Double_t>> LauGenInfo;
//! container of background DP models
typedef std::vector<LauAbsBkgndDPModel*> LauBkgndDPModelList;
//! container of background PDFs
typedef std::vector<LauPdfPList> LauBkgndPdfsList;
//! container of background yields or asymmetries
typedef std::vector<LauAbsRValue*> LauBkgndYieldList;
//! container of boolean toggles for each background category
typedef std::vector<Bool_t> LauBkgndReuseEventsList;
//! Determine the number of events to generate for each event category
LauGenInfo eventsToGenerate();
//! Generate signal event
/*!
\return whether the generation was successful (kTRUE) or not (kFALSE)
*/
Bool_t generateSignalEvent();
//! Generate background event
/*!
\param [in] bgID ID number of the background class
\return whether the generation was successful (kTRUE) or not (kFALSE)
*/
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(LauPdfPList& 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 LauPdfPList& 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 prefix for the branch names
\param [in] iEvt the event number
*/
Double_t setSPlotNtupleBranchValues(LauPdfPList& extraPdfs, const TString& prefix, const 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; }
//! Calculate the component integrals of the interference term
void calcInterferenceTermIntegrals();
//! Calculate the total integral of the interference term
void calcInterferenceTermNorm();
private:
//! Keep access to the base class methods that we're further overloading
using LauAbsFitModel::addFitParameters;
//! Add the parameters from a decay time PDF into the fit
/*!
\param [in] decayTimePdf the container of PDFs
*/
UInt_t addFitParameters(LauDecayTimePdf* decayTimePdf);
//! Add the parameters from each decay time PDF into the fit
/*!
\param [in] decayTimePdfList the container of PDFs
*/
UInt_t addFitParameters(std::vector<LauDecayTimePdf*>& decayTimePdfList);
//! 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 background Dalitz plot models for particles
LauBkgndDPModelList BkgndDPModelsB_;
//! The background Dalitz plot models for anti-particles
LauBkgndDPModelList BkgndDPModelsBbar_;
//! The background PDFs
LauBkgndPdfsList BkgndPdfs_;
//! Background boolean
Bool_t usingBkgnd_;
//! LauFlavTag object for flavour tagging
LauFlavTag* flavTag_;
//! Flavour tag for current event
std::vector<LauFlavTag::Flavour> curEvtTagFlv_;
//! Per event mistag for current event
std::vector<Double_t> curEvtMistag_;
//! Per event transformed mistag for current event
std::vector<Double_t> curEvtMistagPrime_;
//! True tag flavour (i.e. flavour at production) for the current event (only relevant for toy generation)
LauFlavTag::Flavour curEvtTrueTagFlv_;
//! Flavour at decay for the current event (only relevant for QFS)
LauFlavTag::Flavour curEvtDecayFlv_;
//! 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 tagging efficneyc parameters
UInt_t nTagEffPar_;
//! 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<LauComplex> coeffsB0bar_;
//! The complex coefficients for particle
std::vector<LauComplex> coeffsB0_;
//! Magnitudes and Phases or Real and Imaginary Parts
std::vector<LauAbsCoeffSet*> 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<LauComplex> > 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<LauParameter> fitFracAsymm_;
//! A_CP parameter
std::vector<LauParameter> 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<CPEigenvalue> 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_;
//! Signal Decay time PDFs (one per tagging category)
LauDecayTimePdf* signalDecayTimePdf_;
//! Background types
std::vector<LauFlavTag::BkgndType> BkgndTypes_;
//! Background Decay time PDFs (one per tagging category)
std::vector<LauDecayTimePdf*> BkgndDecayTimePdfs_;
//! Decay time for the current event
Double_t curEvtDecayTime_;
//! Decay time error for the current event
Double_t curEvtDecayTimeErr_;
//! PDFs for other variables
LauPdfPList sigExtraPdf_;
//! Production asymmetry between B0 and B0bar
LauParameter AProd_;
//! Production asymmetry between B0 and B0bar for bkgnds
std::vector<LauParameter*> AProdBkgnd_;
// 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
LauEmbeddedData* signalTree_;
//! The background event tree for embedding fully simulated events
std::vector<LauEmbeddedData*> bkgndTree_;
//! Boolean to control reuse of embedded signal events
Bool_t reuseSignal_;
//! Boolean to use reweighting when embedding
Bool_t useReweighting_;
//! 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_;
//! Total signal likelihood
Double_t sigTotalLike_;
//! Background DP likelihood value(s)
std::vector<Double_t> bkgndDPLike_;
//! Background likelihood value(s) from extra PDFs
std::vector<Double_t> bkgndExtraLike_;
//! Total background likelihood(s)
std::vector<Double_t> bkgndTotalLike_;
ClassDef(LauTimeDepFitModel,0) // Time-dependent neutral model
};
#endif

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 11:27 AM (14 h, 39 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111430
Default Alt Text
(25 KB)

Event Timeline