Page MenuHomeHEPForge

No OneTemporary

diff --git a/inc/LauFlavTag.hh b/inc/LauFlavTag.hh
index e76c7fc..bc34323 100644
--- a/inc/LauFlavTag.hh
+++ b/inc/LauFlavTag.hh
@@ -1,194 +1,197 @@
-
-// Copyright University of Warwick 2006 - 2014.
-// Distributed under the Boost Software License, Version 1.0.
-// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
-
-// Authors:
-// Thomas Latham
-// John Back
-// Paul Harrison
-
-/*! \file LauFlavTag.hh
- \brief File containing declaration of LauFlavTag class.
- */
-
-/*! \class LauFlavTag
- \brief Class for defining the flavour tagging approach.
-
- Define the flavour tagging categories and all associated parameters
- to be passed to the relevant fit models.
- */
-
-#ifndef LAU_FLAVTAG
-#define LAU_FLAVTAG
-
-#include <vector>
-#include <map>
-#include <set>
-#include <iostream>
-
-#include "TString.h"
-#include "TStopwatch.h"
-#include "TSystem.h"
-
-#include "LauConstants.hh"
-#include "LauParameter.hh"
-#include "LauFitDataTree.hh"
-#include "LauAbsFitModel.hh"
-
-class LauFlavTag {
-
- public:
- //! Constructor
- /*!
- \param [in] modelB0bar DP model for the antiparticle
- \param [in] modelB0 DP model for the particle
- \param [in] useUntaggedEvents should the untagged sample be used or excluded?
- \param [in] tagVarName the variable name in the data tree that specifies the event tag
- \param [in] tagCatVarName the variable name in the data tree that specifies the event tagging category
- */
- LauFlavTag(const std::vector<LauParameter*>& params, const Bool_t useUntaggedEvents = kTRUE, const TString& tagVarName = "tagFlv", const TString& tagCatVarName = "tagCat");
-
- //! Destructor
- virtual ~LauFlavTag();
-
- void initialise();
-
- //! Add a tagging category to the list of valid categories
- /*!
- NB category 0 is always valid and corresponds to untagged events.
- Whether untagged events are used in the fit or note is controlled by a constructor argument.
-
- \param [in] tagCat the tagging category ID
- */
- void addValidTagCat(const Int_t tagCat);
-
- //! Add several tagging categories to the list of valid categories
- /*!
- NB category 0 is always valid and corresponds to untagged events.
- Whether untagged events are used in the fit or note is controlled by a constructor argument.
-
- \param [in] tagCats the list of tagging category IDs
- */
- void addValidTagCats(const std::vector<Int_t>& tagCats);
-
- //! Check the validity of the given tagging category
- /*!
- \param [in] tagCat the tagging category ID
- */
- Bool_t validTagCat(const Int_t tagCat) const;
-
- //! Change the dilutions, delta dilutions and tagCatFrac for signal if needed
- /*!
- \param [in] tagCat the tagging category to adjust
- \param [in] tagCatFrac the tagging category fraction
- \param [in] dilution the tagging category average dilution = (1 - 2 * avg_mistag_fraction)
- \param [in] deltaDilution the tagging category dilution difference TODO - check sign convention
- \param [in] fixTCFrac whether to fix or float the tagging category fraction
- \param [in] usePerEvtMistag whether to use per event mistag information or not
- \param [in] mistagVarName the name of the branch in the tree containing per event mistag information
- */
- void setSignalTagCatPars(const Int_t tagCat, const Double_t tagCatFrac, const Double_t dilution, const Double_t deltaDilution,
- const Bool_t fixTCFrac = kTRUE, const Bool_t usePerEvtMistag = kFALSE, const TString& mistagVarName = "tagMistag");
-
- //! Read in the input fit data variables, e.g. m13Sq and m23Sq
- void cacheInputFitVars(LauFitDataTree* inputFitData);
-
- Bool_t getUsePerEvtMistag(){return usePerEvtMistag_;};
-
- typedef std::map< Int_t, LauParameter> LauTagCatParamMap;
-
- LauTagCatParamMap getDilution(){return dilution_;};
- LauTagCatParamMap getDeltaDilution(){return deltaDilution_;};
- LauTagCatParamMap& getSignalTagCatFrac(){return signalTagCatFrac_;};
-
- //! Calculates the fraction of the first tagging category based on the others
- /*!
- \param [in,out] theMap the container of tagcat fractions
- */
- void setFirstTagCatFrac(LauTagCatParamMap& theMap);
-
- //std::vector<Int_t> getEvtTagFlvVals(){return evtTagFlvVals_;};
- //std::vector<Int_t> getEvtTagCatVals(){return evtTagCatVals_;};
- //std::vector<Double_t> getEvtMistagVals(){return evtMistagVals_;};
- Int_t getEvtTagFlvVals(Int_t ievt){return evtTagFlvVals_[ievt];};
- Int_t getEvtTagCatVals(Int_t ievt){return evtTagCatVals_[ievt];};
- Double_t getEvtMistagVals(Int_t ievt){return evtMistagVals_[ievt];};
- TString getMistagVarName(){return mistagVarName_;};
-
- LauParameter* findParameter(const TString& parName);
- std::vector<LauParameter*> getCalibParameters();
-
- Double_t getPerEvtAvgMistag(){return perEvtAvgMistag_;};
- Double_t getOmega(Int_t ievt);
-
- protected:
- //! Check that the tagging category fractions are all present and sum to unity
- /*!
- \param [in] theMap the container of tagcat fractions
- */
- Bool_t checkTagCatFracMap(const LauTagCatParamMap& theMap) const;
-
- //! Check the signal tagging category fractions
- void checkSignalTagCatFractions();
-
- private:
- //! Whether or not to use untagged events
- const Bool_t useUntaggedEvents_;
-
- //! Signal tagging category fractions
- LauTagCatParamMap signalTagCatFrac_;
-
- //! Flavour tagging variable name
- TString tagVarName_;
-
- //! Tagging category variable name
- TString tagCatVarName_;
-
- //! Per event mistag variable name
- TString mistagVarName_;
-
- //! The allowed tagging categories
- std::set<Int_t> validTagCats_;
-
- //! 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_;
-
- //! Vector to store event flavour tags
- std::vector<Int_t> evtTagFlvVals_;
-
- //! Vector to store event tagging categories
- std::vector<Int_t> evtTagCatVals_;
-
- //! Vector to store event mistag values
- std::vector<Double_t> evtMistagVals_;
-
- //! Average dilutions
- LauTagCatParamMap dilution_;
-
- //! Dilution differences
- LauTagCatParamMap deltaDilution_;
-
- //! Use per event mistag
- Bool_t usePerEvtMistag_;
-
- //! Per-event average mistag value (eta hat)
- Double_t perEvtAvgMistag_;
-
- //! Calibration parameters
- LauParameter* calib_p0_;
- LauParameter* calib_p1_;
-
- //! Input parameters
- std::vector<LauParameter*> params_;
-
- ClassDef(LauFlavTag,0) // Flaour tagging set up
-};
-
-#endif
+
+// Copyright University of Warwick 2006 - 2014.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Authors:
+// Thomas Latham
+// John Back
+// Paul Harrison
+
+/*! \file LauFlavTag.hh
+ \brief File containing declaration of LauFlavTag class.
+ */
+
+/*! \class LauFlavTag
+ \brief Class for defining the flavour tagging approach.
+
+ Define the flavour tagging categories and all associated parameters
+ to be passed to the relevant fit models.
+ */
+
+#ifndef LAU_FLAVTAG
+#define LAU_FLAVTAG
+
+#include <vector>
+#include <map>
+#include <set>
+#include <iostream>
+
+#include "TString.h"
+#include "TStopwatch.h"
+#include "TSystem.h"
+
+#include "LauConstants.hh"
+#include "LauParameter.hh"
+#include "LauFitDataTree.hh"
+#include "LauAbsFitModel.hh"
+
+class LauFlavTag {
+
+ public:
+ //! Constructor
+ /*!
+ \param [in] modelB0bar DP model for the antiparticle
+ \param [in] modelB0 DP model for the particle
+ \param [in] useUntaggedEvents should the untagged sample be used or excluded?
+ \param [in] tagVarName the variable name in the data tree that specifies the event tag
+ \param [in] tagCatVarName the variable name in the data tree that specifies the event tagging category
+ */
+ LauFlavTag(const std::vector<LauParameter*>& params, const Bool_t useUntaggedEvents = kTRUE, const TString& tagVarName = "tagFlv", const TString& tagCatVarName = "tagCat");
+
+ //! Destructor
+ virtual ~LauFlavTag();
+
+ void initialise();
+
+ //! Add a tagging category to the list of valid categories
+ /*!
+ NB category 0 is always valid and corresponds to untagged events.
+ Whether untagged events are used in the fit or note is controlled by a constructor argument.
+
+ \param [in] tagCat the tagging category ID
+ */
+ void addValidTagCat(const Int_t tagCat);
+
+ //! Add several tagging categories to the list of valid categories
+ /*!
+ NB category 0 is always valid and corresponds to untagged events.
+ Whether untagged events are used in the fit or note is controlled by a constructor argument.
+
+ \param [in] tagCats the list of tagging category IDs
+ */
+ void addValidTagCats(const std::vector<Int_t>& tagCats);
+
+ //! Check the validity of the given tagging category
+ /*!
+ \param [in] tagCat the tagging category ID
+ */
+ Bool_t validTagCat(const Int_t tagCat) const;
+
+ //! Change the dilutions, delta dilutions and tagCatFrac for signal if needed
+ /*!
+ \param [in] tagCat the tagging category to adjust
+ \param [in] tagCatFrac the tagging category fraction
+ \param [in] dilution the tagging category average dilution = (1 - 2 * avg_mistag_fraction)
+ \param [in] deltaDilution the tagging category dilution difference TODO - check sign convention
+ \param [in] fixTCFrac whether to fix or float the tagging category fraction
+ \param [in] usePerEvtMistag whether to use per event mistag information or not
+ \param [in] mistagVarName the name of the branch in the tree containing per event mistag information
+ */
+ void setSignalTagCatPars(const Int_t tagCat, const Double_t tagCatFrac, const Double_t dilution, const Double_t deltaDilution,
+ const Bool_t fixTCFrac = kTRUE, const Bool_t usePerEvtMistag = kFALSE, const TString& mistagVarName = "tagMistag");
+
+ //! Read in the input fit data variables, e.g. m13Sq and m23Sq
+ void cacheInputFitVars(LauFitDataTree* inputFitData);
+
+ Bool_t getUsePerEvtMistag(){return usePerEvtMistag_;};
+
+ typedef std::map< Int_t, LauParameter> LauTagCatParamMap;
+
+ LauTagCatParamMap getDilution(){return dilution_;};
+ LauTagCatParamMap getDeltaDilution(){return deltaDilution_;};
+ LauTagCatParamMap& getSignalTagCatFrac(){return signalTagCatFrac_;};
+
+ //! Update the fraction of the first tagging category for all species
+ void setFirstTagCatFractions();
+
+ //std::vector<Int_t> getEvtTagFlvVals(){return evtTagFlvVals_;};
+ //std::vector<Int_t> getEvtTagCatVals(){return evtTagCatVals_;};
+ //std::vector<Double_t> getEvtMistagVals(){return evtMistagVals_;};
+ Int_t getEvtTagFlvVals(UInt_t iEvt){return evtTagFlvVals_[iEvt];};
+ Int_t getEvtTagCatVals(UInt_t iEvt){return evtTagCatVals_[iEvt];};
+ Double_t getEvtMistagVals(UInt_t iEvt){return evtMistagVals_[iEvt];};
+ const TString& getMistagVarName(){return mistagVarName_;};
+
+ LauParameter* findParameter(const TString& parName);
+ std::vector<LauParameter*> getCalibParameters();
+
+ Double_t getPerEvtAvgMistag() const {return perEvtAvgMistag_;};
+ Double_t getOmega(const UInt_t ievt) const;
+
+ protected:
+ //! Check the signal tagging category fractions
+ void checkSignalTagCatFractions();
+
+ //! Check that the background tagging category fractions are all present and sum to unity
+ /*!
+ \param [in] theMap the container of tagcat fractions
+ */
+ Bool_t checkTagCatFracMap(const LauTagCatParamMap& theMap) const;
+
+ //! Calculates the fraction of the first tagging category based on the others
+ /*!
+ \param [in] theMap the container of tagcat fractions
+ */
+ void setFirstTagCatFrac(LauTagCatParamMap& theMap);
+
+ private:
+ //! Whether or not to use untagged events
+ const Bool_t useUntaggedEvents_;
+
+ //! Signal tagging category fractions
+ LauTagCatParamMap signalTagCatFrac_;
+
+ //! Flavour tagging variable name
+ TString tagVarName_;
+
+ //! Tagging category variable name
+ TString tagCatVarName_;
+
+ //! Per event mistag variable name
+ TString mistagVarName_;
+
+ //! The allowed tagging categories
+ std::set<Int_t> validTagCats_;
+
+ //! 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_;
+
+ //! Vector to store event flavour tags
+ std::vector<Int_t> evtTagFlvVals_;
+
+ //! Vector to store event tagging categories
+ std::vector<Int_t> evtTagCatVals_;
+
+ //! Vector to store event mistag values
+ std::vector<Double_t> evtMistagVals_;
+
+ //! Average dilutions
+ LauTagCatParamMap dilution_;
+
+ //! Dilution differences
+ LauTagCatParamMap deltaDilution_;
+
+ //! Use per event mistag
+ Bool_t usePerEvtMistag_;
+
+ //! Per-event average mistag value (eta hat)
+ Double_t perEvtAvgMistag_;
+
+ //! Calibration parameters
+ LauParameter* calib_p0_;
+ LauParameter* calib_p1_;
+
+ //! Input parameters
+ std::vector<LauParameter*> params_;
+
+ ClassDef(LauFlavTag,0) // Flaour tagging set up
+};
+
+#endif
diff --git a/src/LauFlavTag.cc b/src/LauFlavTag.cc
index 234a8fe..7ef4107 100644
--- a/src/LauFlavTag.cc
+++ b/src/LauFlavTag.cc
@@ -1,344 +1,346 @@
-
-// Copyright University of Warwick 2006 - 2014.
-// Distributed under the Boost Software License, Version 1.0.
-// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
-
-// Authors:
-// Thomas Latham
-// John Back
-// Paul Harrison
-
-/*! \file LauFlavTag.cc
- \brief File containing implementation of LauFlavTag class.
- */
-
-#include <iostream>
-#include <iomanip>
-#include <fstream>
-#include <map>
-#include <vector>
-
-#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 "LauFlavTag.hh"
-
-ClassImp(LauFlavTag)
-
- LauFlavTag::LauFlavTag(const std::vector<LauParameter*>& params, const Bool_t useUntaggedEvents, const TString& tagVarName, const TString& tagCatVarName) :
- useUntaggedEvents_(useUntaggedEvents),
- signalTagCatFrac_(),
- tagVarName_(tagVarName),
- tagCatVarName_(tagCatVarName),
- mistagVarName_("tagMistag"),
- validTagCats_(),
- curEvtTagFlv_(0),
- curEvtTagCat_(0),
- curEvtMistag_(0.),
- evtTagFlvVals_(0),
- evtTagCatVals_(0),
- evtMistagVals_(0),
- dilution_(),
- deltaDilution_(),
- usePerEvtMistag_(kFALSE),
- perEvtAvgMistag_(),
- calib_p0_(),
- calib_p1_(),
- params_(params)
-{
- // Add the untagged category as a valid category
- this->addValidTagCat(0);
-
- // Set the fraction, average dilution and dilution difference for the untagged category
- this->setSignalTagCatPars(0, 1.0, 0.0, 0.0, kTRUE);
-}
-
-LauFlavTag::~LauFlavTag()
-{
-}
-
-void LauFlavTag::initialise()
-{
- this->checkSignalTagCatFractions();
-}
-
-void LauFlavTag::addValidTagCats(const std::vector<Int_t>& tagCats)
-{
- for (std::vector<Int_t>::const_iterator iter = tagCats.begin(); iter != tagCats.end(); ++iter) {
- this->addValidTagCat(*iter);
- }
-}
-
-void LauFlavTag::addValidTagCat(const Int_t tagCat)
-{
- validTagCats_.insert(tagCat);
-}
-
-void LauFlavTag::setSignalTagCatPars(const Int_t tagCat, const Double_t tagCatFrac, const Double_t dilution, const Double_t deltaDilution, const Bool_t fixTCFrac, const Bool_t usePerEvtMistag, const TString& mistagVarName)
-{
- if (!this->validTagCat(tagCat)) {
- std::cerr<<"ERROR in LauFlavTag::setSignalTagCatPars : Tagging category \""<<tagCat<<"\" not valid, not changing the parameters."<<std::endl;
- return;
- }
-
- TString tagCatFracName("signalTagCatFrac");
- tagCatFracName += tagCat;
-
- signalTagCatFrac_[tagCat].name(tagCatFracName);
- signalTagCatFrac_[tagCat].range(0.0,1.0);
- signalTagCatFrac_[tagCat].value(tagCatFrac); signalTagCatFrac_[tagCat].initValue(tagCatFrac); signalTagCatFrac_[tagCat].genValue(tagCatFrac);
- signalTagCatFrac_[tagCat].fixed(fixTCFrac);
-
- TString dilutionName("dilution");
- dilutionName += tagCat;
-
- dilution_[tagCat].name(dilutionName);
- dilution_[tagCat].range(0.0,1.0);
- dilution_[tagCat].value(dilution); dilution_[tagCat].initValue(dilution); dilution_[tagCat].genValue(dilution);
- dilution_[tagCat].fixed(kTRUE);
-
- TString deltaDilutionName("deltaDilution");
- deltaDilutionName += tagCat;
-
- deltaDilution_[tagCat].name(deltaDilutionName);
- deltaDilution_[tagCat].range(-2.0,2.0);
- deltaDilution_[tagCat].value(deltaDilution); deltaDilution_[tagCat].initValue(deltaDilution); deltaDilution_[tagCat].genValue(deltaDilution);
- deltaDilution_[tagCat].fixed(kTRUE);
-
- if (usePerEvtMistag){
- usePerEvtMistag_ = kTRUE;
- mistagVarName_ = mistagVarName;
- std::cout<<"INFO in LauFlavTag::setSignalTagCatPars : Using per event mistag information from branch "<<mistagVarName_<<" ."<<std::endl;
- }
- // We check they're set up correctly with checkSignalTagCatFractions()
- // only when the user has set them all up, see initialise()
-}
-
-void LauFlavTag::checkSignalTagCatFractions()
-{
- Double_t totalTaggedFrac(0.0);
- for (LauTagCatParamMap::const_iterator iter=signalTagCatFrac_.begin(); iter!=signalTagCatFrac_.end(); ++iter) {
- if (iter->first != 0) {
- const LauParameter& par = iter->second;
- totalTaggedFrac += par.value();
- }
- }
-
- if ( ((totalTaggedFrac < (1.0-1.0e-8))&&!useUntaggedEvents_) || (totalTaggedFrac > (1.0+1.0e-8)) ) {
- std::cerr<<"WARNING in LauFlavTag::checkSignalTagCatFractions : Tagging category fractions add up to "<<totalTaggedFrac<<", not 1.0; normalizing them."<<std::endl;
- for (LauTagCatParamMap::iterator iter=signalTagCatFrac_.begin(); iter!=signalTagCatFrac_.end(); ++iter) {
- LauParameter& par = iter->second;
- Double_t newVal = par.value() / totalTaggedFrac;
- par.value(newVal); par.initValue(newVal); par.genValue(newVal);
- }
- } else if (useUntaggedEvents_) {
- Double_t tagCatFrac = 1.0 - totalTaggedFrac;
- TString tagCatFracName("signalTagCatFrac0");
- signalTagCatFrac_[0].name(tagCatFracName);
- signalTagCatFrac_[0].range(0.0,1.0);
- signalTagCatFrac_[0].value(tagCatFrac); signalTagCatFrac_[0].initValue(tagCatFrac); signalTagCatFrac_[0].genValue(tagCatFrac);
- signalTagCatFrac_[0].fixed(kTRUE);
-
- TString dilutionName("dilution0");
- dilution_[0].name(dilutionName);
- dilution_[0].range(0.0,1.0);
- dilution_[0].value(0.0); dilution_[0].initValue(0.0); dilution_[0].genValue(0.0);
- dilution_[0].fixed(kTRUE);
-
- TString deltaDilutionName("deltaDilution0");
- deltaDilution_[0].name(deltaDilutionName);
- deltaDilution_[0].range(-2.0,2.0);
- deltaDilution_[0].value(0.0); deltaDilution_[0].initValue(0.0); deltaDilution_[0].genValue(0.0);
- deltaDilution_[0].fixed(kTRUE);
- }
-
- for (LauTagCatParamMap::const_iterator iter=dilution_.begin(); iter!=dilution_.end(); ++iter) {
- std::cout<<"INFO in LauFlavTag::checkSignalTagCatFractions : Setting dilution for tagging category "<<(*iter).first<<" to "<<(*iter).second<<std::endl;
- }
- for (LauTagCatParamMap::const_iterator iter=deltaDilution_.begin(); iter!=deltaDilution_.end(); ++iter) {
- std::cout<<"INFO in LauFlavTag::checkSignalTagCatFractions : Setting Delta(dilution) for tagging category "<<(*iter).first<<" to "<<(*iter).second<<std::endl;
- }
- for (LauTagCatParamMap::const_iterator iter=signalTagCatFrac_.begin(); iter!=signalTagCatFrac_.end(); ++iter) {
- std::cout<<"INFO in LauFlavTag::checkSignalTagCatFractions : Setting (signal) fraction for tagging category "<<(*iter).first<<" to "<<(*iter).second<<std::endl;
- }
-}
-
-void LauFlavTag::setFirstTagCatFrac(LauTagCatParamMap& theMap)
-{
- Double_t firstCatFrac = 1.0;
- Int_t firstCat(0);
- for (LauTagCatParamMap::iterator iter = theMap.begin(); iter != theMap.end(); ++iter) {
- if (iter == theMap.begin()) {
- firstCat = iter->first;
- continue;
- }
- LauParameter& par = iter->second;
- firstCatFrac -= par.unblindValue();
- }
- theMap[firstCat].value(firstCatFrac);
-}
-
-Bool_t LauFlavTag::validTagCat(Int_t tagCat) const
-{
- return (validTagCats_.find(tagCat) != validTagCats_.end());
-}
-
-Bool_t LauFlavTag::checkTagCatFracMap(const LauTagCatParamMap& theMap) const
-{
- // First check that there is an entry for each tagging category.
- // NB an entry won't have been added if it isn't a valid category
- // so don't need to check for that here.
- if (theMap.size() != signalTagCatFrac_.size()) {
- std::cerr<<"ERROR in LauFlavTag::checkTagCatFracMap : Not all tagging categories present."<<std::endl;
- return kFALSE;
- }
-
- // Now check that the fractions sum up to unity.
- Double_t tagCatFracSum(0.0);
- for (LauTagCatParamMap::const_iterator iter = theMap.begin(); iter != theMap.end(); ++iter) {
- const LauParameter& par = (*iter).second;
- tagCatFracSum += par.unblindValue();
- }
- if ((tagCatFracSum - 1.0) > 1E-10) {
- std::cerr<<"ERROR in LauFlavTag::checkTagCatFracMap : Tagging category event fractions do not sum to unity."<<std::endl;
- return kFALSE;
- }
-
- // If we've got to here then all is OK.
- return kTRUE;
-}
-
-void LauFlavTag::cacheInputFitVars(LauFitDataTree* inputFitData)
-{
- // Fill event by event data
- // Start by caching the tagging information
- evtTagCatVals_.clear();
- evtTagFlvVals_.clear();
- evtMistagVals_.clear();
- if ( ! inputFitData->haveBranch( tagCatVarName_ ) ) {
- std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << tagCatVarName_ << "\"." << std::endl;
- gSystem->Exit(EXIT_FAILURE);
- }
- if ( ! inputFitData->haveBranch( tagVarName_ ) ) {
- std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << tagVarName_ << "\"." << std::endl;
- gSystem->Exit(EXIT_FAILURE);
- }
- if ( ! inputFitData->haveBranch( mistagVarName_ ) ) {
- std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << mistagVarName_ << "\"." << std::endl;
- gSystem->Exit(EXIT_FAILURE);
- }
- UInt_t nEvents = inputFitData->nEvents();
- evtTagCatVals_.reserve( nEvents );
- evtTagFlvVals_.reserve( nEvents );
- evtMistagVals_.reserve( nEvents );
-
- //Running total to calc average mistag
- Double_t tot_mistag(0.);
-
- Int_t totalTagged = 0;
-
- LauFitData::const_iterator fitdata_iter;
- for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
- const LauFitData& dataValues = inputFitData->getData(iEvt);
-
- fitdata_iter = dataValues.find( tagCatVarName_ );
- curEvtTagCat_ = static_cast<Int_t>( fitdata_iter->second );
- if ( ! this->validTagCat( curEvtTagCat_ ) ) {
- std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging category " << curEvtTagCat_ << " for event " << iEvt << ", setting it to untagged" << std::endl;
- curEvtTagCat_ = 0;
- }
- evtTagCatVals_.push_back( curEvtTagCat_ );
-
- fitdata_iter = dataValues.find( tagVarName_ );
- curEvtTagFlv_ = static_cast<Int_t>( fitdata_iter->second );
- if ( TMath::Abs( curEvtTagFlv_ ) != 1 ) {
- if ( curEvtTagFlv_ > 0 ) {
- std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv_ << " for event " << iEvt << ", setting it to +1" << std::endl;
- curEvtTagFlv_ = +1;
- } else {
- std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv_ << " for event " << iEvt << ", setting it to -1" << std::endl;
- curEvtTagFlv_ = -1;
- }
- }
- evtTagFlvVals_.push_back( curEvtTagFlv_ );
-
- fitdata_iter = dataValues.find( mistagVarName_);
- curEvtMistag_ = static_cast<Double_t>( fitdata_iter->second );
- if (curEvtMistag_ > 0.5){
- std::cerr<<"WARNING in LauFlavTag::cacheInputFitVars : Mistag value "<<curEvtMistag_<<" is larger than 0.5, setting to 0.5"<<std::endl;
- curEvtMistag_ = 0.5;
- }
- if (curEvtMistag_ < 0.0){
- std::cerr<<"WARNING in LauFlavTag::cacheInputFitVars : Mistag value "<<curEvtMistag_<<" is less than 0.0, setting to 0.0"<<std::endl;
- curEvtMistag_ = 0.0;
- }
-
- // Only those with tags (cat != 0) contribute to the average
- if (curEvtTagCat_ != 0) {
- tot_mistag += curEvtMistag_;
- totalTagged += 1;
- }
- evtMistagVals_.push_back( curEvtMistag_ );
- }
- // Store average per event mistag for calibration
- perEvtAvgMistag_ = tot_mistag/totalTagged;
-
- // Check here that the tagging category fractions add up to 1, otherwise "normalise". Also set up the untagged cat.
- // NB this has to be done early in the initialization as other methods access the tagCats map.
- //this->checkSignalTagCatFractions();
-
-}
-
-LauParameter* LauFlavTag::findParameter(const TString& parName)
-{
- for ( std::vector<LauParameter*>::iterator iter = params_.begin(); iter != params_.end(); ++iter ) {
- if ((*iter)->name().Contains(parName)) {
- return (*iter);
- }
- }
- std::cerr << "ERROR in LauFlavTag::findParameter : Parameter \"" << parName << "\" not found." << std::endl;
- return 0;
-}
-
-Double_t LauFlavTag::getOmega(Int_t ievt)
-{
- return calib_p0_->value() + calib_p1_->value()*( evtMistagVals_[ievt] - perEvtAvgMistag_);
- // return calib_p0_->value() + calib_p1_->value()*( evtMistagVals_[ievt] );
-}
-
-std::vector<LauParameter*> LauFlavTag::getCalibParameters()
-{
- std::vector<LauParameter*> calibParams;
- calib_p0_ = this->findParameter("calib_p0");
- calib_p1_ = this->findParameter("calib_p1");
-
- // calib_p0_->value(perEvtAvgMistag_); // WOULD BE NICE TO DO THIS
-
- calibParams.push_back( calib_p0_ );
- calibParams.push_back( calib_p1_ );
- if (calibParams.size() != 2) {
- std::cout << "FATAL in LauFlavTag::setCalibParams : Expected two calibration parameters, received " << calibParams.size() << "." << std::endl;
- exit(1);
- }
- return calibParams;
-
-}
+
+// Copyright University of Warwick 2006 - 2014.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Authors:
+// Thomas Latham
+// John Back
+// Paul Harrison
+
+/*! \file LauFlavTag.cc
+ \brief File containing implementation of LauFlavTag class.
+ */
+
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <map>
+#include <vector>
+
+#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 "LauFlavTag.hh"
+
+ClassImp(LauFlavTag)
+
+ LauFlavTag::LauFlavTag(const std::vector<LauParameter*>& params, const Bool_t useUntaggedEvents, const TString& tagVarName, const TString& tagCatVarName) :
+ useUntaggedEvents_(useUntaggedEvents),
+ signalTagCatFrac_(),
+ tagVarName_(tagVarName),
+ tagCatVarName_(tagCatVarName),
+ mistagVarName_("tagMistag"),
+ validTagCats_(),
+ curEvtTagFlv_(0),
+ curEvtTagCat_(0),
+ curEvtMistag_(0.),
+ evtTagFlvVals_(0),
+ evtTagCatVals_(0),
+ evtMistagVals_(0),
+ dilution_(),
+ deltaDilution_(),
+ usePerEvtMistag_(kFALSE),
+ perEvtAvgMistag_(),
+ calib_p0_(),
+ calib_p1_(),
+ params_(params)
+{
+ // Add the untagged category as a valid category
+ this->addValidTagCat(0);
+
+ // Set the fraction, average dilution and dilution difference for the untagged category
+ this->setSignalTagCatPars(0, 1.0, 0.0, 0.0, kTRUE);
+}
+
+LauFlavTag::~LauFlavTag()
+{
+}
+
+void LauFlavTag::initialise()
+{
+ this->checkSignalTagCatFractions();
+}
+
+void LauFlavTag::addValidTagCats(const std::vector<Int_t>& tagCats)
+{
+ for (std::vector<Int_t>::const_iterator iter = tagCats.begin(); iter != tagCats.end(); ++iter) {
+ this->addValidTagCat(*iter);
+ }
+}
+
+void LauFlavTag::addValidTagCat(const Int_t tagCat)
+{
+ validTagCats_.insert(tagCat);
+}
+
+void LauFlavTag::setSignalTagCatPars(const Int_t tagCat, const Double_t tagCatFrac, const Double_t dilution, const Double_t deltaDilution, const Bool_t fixTCFrac, const Bool_t usePerEvtMistag, const TString& mistagVarName)
+{
+ if (!this->validTagCat(tagCat)) {
+ std::cerr<<"ERROR in LauFlavTag::setSignalTagCatPars : Tagging category \""<<tagCat<<"\" not valid, not changing the parameters."<<std::endl;
+ return;
+ }
+
+ TString tagCatFracName("signalTagCatFrac");
+ tagCatFracName += tagCat;
+
+ signalTagCatFrac_[tagCat].name(tagCatFracName);
+ signalTagCatFrac_[tagCat].range(0.0,1.0);
+ signalTagCatFrac_[tagCat].value(tagCatFrac); signalTagCatFrac_[tagCat].initValue(tagCatFrac); signalTagCatFrac_[tagCat].genValue(tagCatFrac);
+ signalTagCatFrac_[tagCat].fixed(fixTCFrac);
+
+ TString dilutionName("dilution");
+ dilutionName += tagCat;
+
+ dilution_[tagCat].name(dilutionName);
+ dilution_[tagCat].range(0.0,1.0);
+ dilution_[tagCat].value(dilution); dilution_[tagCat].initValue(dilution); dilution_[tagCat].genValue(dilution);
+ dilution_[tagCat].fixed(kTRUE);
+
+ TString deltaDilutionName("deltaDilution");
+ deltaDilutionName += tagCat;
+
+ deltaDilution_[tagCat].name(deltaDilutionName);
+ deltaDilution_[tagCat].range(-2.0,2.0);
+ deltaDilution_[tagCat].value(deltaDilution); deltaDilution_[tagCat].initValue(deltaDilution); deltaDilution_[tagCat].genValue(deltaDilution);
+ deltaDilution_[tagCat].fixed(kTRUE);
+
+ if (usePerEvtMistag){
+ usePerEvtMistag_ = kTRUE;
+ mistagVarName_ = mistagVarName;
+ std::cout<<"INFO in LauFlavTag::setSignalTagCatPars : Using per event mistag information from branch "<<mistagVarName_<<" ."<<std::endl;
+ }
+ // We check they're set up correctly with checkSignalTagCatFractions()
+ // only when the user has set them all up, see initialise()
+}
+
+void LauFlavTag::checkSignalTagCatFractions()
+{
+ Double_t totalTaggedFrac(0.0);
+ for (LauTagCatParamMap::const_iterator iter=signalTagCatFrac_.begin(); iter!=signalTagCatFrac_.end(); ++iter) {
+ if (iter->first != 0) {
+ const LauParameter& par = iter->second;
+ totalTaggedFrac += par.value();
+ }
+ }
+
+ if ( ((totalTaggedFrac < (1.0-1.0e-8))&&!useUntaggedEvents_) || (totalTaggedFrac > (1.0+1.0e-8)) ) {
+ std::cerr<<"WARNING in LauFlavTag::checkSignalTagCatFractions : Tagging category fractions add up to "<<totalTaggedFrac<<", not 1.0; normalizing them."<<std::endl;
+ for (LauTagCatParamMap::iterator iter=signalTagCatFrac_.begin(); iter!=signalTagCatFrac_.end(); ++iter) {
+ LauParameter& par = iter->second;
+ Double_t newVal = par.value() / totalTaggedFrac;
+ par.value(newVal); par.initValue(newVal); par.genValue(newVal);
+ }
+ } else if (useUntaggedEvents_) {
+ Double_t tagCatFrac = 1.0 - totalTaggedFrac;
+ TString tagCatFracName("signalTagCatFrac0");
+ signalTagCatFrac_[0].name(tagCatFracName);
+ signalTagCatFrac_[0].range(0.0,1.0);
+ signalTagCatFrac_[0].value(tagCatFrac); signalTagCatFrac_[0].initValue(tagCatFrac); signalTagCatFrac_[0].genValue(tagCatFrac);
+ signalTagCatFrac_[0].fixed(kTRUE);
+
+ TString dilutionName("dilution0");
+ dilution_[0].name(dilutionName);
+ dilution_[0].range(0.0,1.0);
+ dilution_[0].value(0.0); dilution_[0].initValue(0.0); dilution_[0].genValue(0.0);
+ dilution_[0].fixed(kTRUE);
+
+ TString deltaDilutionName("deltaDilution0");
+ deltaDilution_[0].name(deltaDilutionName);
+ deltaDilution_[0].range(-2.0,2.0);
+ deltaDilution_[0].value(0.0); deltaDilution_[0].initValue(0.0); deltaDilution_[0].genValue(0.0);
+ deltaDilution_[0].fixed(kTRUE);
+ }
+
+ for (LauTagCatParamMap::const_iterator iter=dilution_.begin(); iter!=dilution_.end(); ++iter) {
+ std::cout<<"INFO in LauFlavTag::checkSignalTagCatFractions : Setting dilution for tagging category "<<(*iter).first<<" to "<<(*iter).second<<std::endl;
+ }
+ for (LauTagCatParamMap::const_iterator iter=deltaDilution_.begin(); iter!=deltaDilution_.end(); ++iter) {
+ std::cout<<"INFO in LauFlavTag::checkSignalTagCatFractions : Setting Delta(dilution) for tagging category "<<(*iter).first<<" to "<<(*iter).second<<std::endl;
+ }
+ for (LauTagCatParamMap::const_iterator iter=signalTagCatFrac_.begin(); iter!=signalTagCatFrac_.end(); ++iter) {
+ std::cout<<"INFO in LauFlavTag::checkSignalTagCatFractions : Setting (signal) fraction for tagging category "<<(*iter).first<<" to "<<(*iter).second<<std::endl;
+ }
+}
+
+void LauFlavTag::setFirstTagCatFractions()
+{
+ this->setFirstTagCatFrac(signalTagCatFrac_);
+ // TODO - add background maps
+}
+
+void LauFlavTag::setFirstTagCatFrac(LauTagCatParamMap& theMap)
+{
+ Double_t firstCatFrac = 1.0;
+ Int_t firstCat(0);
+ for (LauTagCatParamMap::iterator iter = theMap.begin(); iter != theMap.end(); ++iter) {
+ if (iter == theMap.begin()) {
+ firstCat = iter->first;
+ continue;
+ }
+ LauParameter& par = iter->second;
+ firstCatFrac -= par.unblindValue();
+ }
+ theMap[firstCat].value(firstCatFrac);
+}
+
+Bool_t LauFlavTag::validTagCat(Int_t tagCat) const
+{
+ return (validTagCats_.find(tagCat) != validTagCats_.end());
+}
+
+Bool_t LauFlavTag::checkTagCatFracMap(const LauTagCatParamMap& theMap) const
+{
+ // TODO - this is for checking the the background maps are OK (so it is unused at the moment)
+
+ // First check that there is an entry for each tagging category.
+ // NB an entry won't have been added if it isn't a valid category
+ // so don't need to check for that here.
+ if (theMap.size() != signalTagCatFrac_.size()) {
+ std::cerr<<"ERROR in LauFlavTag::checkTagCatFracMap : Not all tagging categories present."<<std::endl;
+ return kFALSE;
+ }
+
+ // Now check that the fractions sum up to unity.
+ Double_t tagCatFracSum(0.0);
+ for (LauTagCatParamMap::const_iterator iter = theMap.begin(); iter != theMap.end(); ++iter) {
+ const LauParameter& par = (*iter).second;
+ tagCatFracSum += par.unblindValue();
+ }
+ if ((tagCatFracSum - 1.0) > 1E-10) {
+ std::cerr<<"ERROR in LauFlavTag::checkTagCatFracMap : Tagging category event fractions do not sum to unity."<<std::endl;
+ return kFALSE;
+ }
+
+ // If we've got to here then all is OK.
+ return kTRUE;
+}
+
+void LauFlavTag::cacheInputFitVars(LauFitDataTree* inputFitData)
+{
+ // Fill event by event data
+ // Start by caching the tagging information
+ evtTagCatVals_.clear();
+ evtTagFlvVals_.clear();
+ evtMistagVals_.clear();
+ if ( ! inputFitData->haveBranch( tagCatVarName_ ) ) {
+ std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << tagCatVarName_ << "\"." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ if ( ! inputFitData->haveBranch( tagVarName_ ) ) {
+ std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << tagVarName_ << "\"." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ if ( ! inputFitData->haveBranch( mistagVarName_ ) ) {
+ std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << mistagVarName_ << "\"." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ UInt_t nEvents = inputFitData->nEvents();
+ evtTagCatVals_.reserve( nEvents );
+ evtTagFlvVals_.reserve( nEvents );
+ evtMistagVals_.reserve( nEvents );
+
+ //Running total to calc average mistag
+ Double_t tot_mistag(0.);
+ Int_t totalTagged = 0;
+
+ LauFitData::const_iterator fitdata_iter;
+ for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
+ const LauFitData& dataValues = inputFitData->getData(iEvt);
+
+ fitdata_iter = dataValues.find( tagCatVarName_ );
+ curEvtTagCat_ = static_cast<Int_t>( fitdata_iter->second );
+ if ( ! this->validTagCat( curEvtTagCat_ ) ) {
+ std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging category " << curEvtTagCat_ << " for event " << iEvt << ", setting it to untagged" << std::endl;
+ curEvtTagCat_ = 0;
+ }
+ evtTagCatVals_.push_back( curEvtTagCat_ );
+
+ fitdata_iter = dataValues.find( tagVarName_ );
+ curEvtTagFlv_ = static_cast<Int_t>( fitdata_iter->second );
+ if ( TMath::Abs( curEvtTagFlv_ ) != 1 ) {
+ if ( curEvtTagFlv_ > 0 ) {
+ std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv_ << " for event " << iEvt << ", setting it to +1" << std::endl;
+ curEvtTagFlv_ = +1;
+ } else {
+ std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv_ << " for event " << iEvt << ", setting it to -1" << std::endl;
+ curEvtTagFlv_ = -1;
+ }
+ }
+ evtTagFlvVals_.push_back( curEvtTagFlv_ );
+
+ // TODO - surely this part should only be done if we have per-event mistagging
+ fitdata_iter = dataValues.find( mistagVarName_);
+ curEvtMistag_ = static_cast<Double_t>( fitdata_iter->second );
+ if (curEvtMistag_ > 0.5){
+ std::cerr<<"WARNING in LauFlavTag::cacheInputFitVars : Mistag value "<<curEvtMistag_<<" is larger than 0.5, setting to 0.5"<<std::endl;
+ curEvtMistag_ = 0.5;
+ }
+ if (curEvtMistag_ < 0.0){
+ std::cerr<<"WARNING in LauFlavTag::cacheInputFitVars : Mistag value "<<curEvtMistag_<<" is less than 0.0, setting to 0.0"<<std::endl;
+ curEvtMistag_ = 0.0;
+ }
+
+ // Only those with tags (cat != 0) contribute to the average
+ if (curEvtTagCat_ != 0) {
+ tot_mistag += curEvtMistag_;
+ totalTagged += 1;
+ }
+ evtMistagVals_.push_back( curEvtMistag_ );
+ }
+ // Store average per event mistag for calibration
+ perEvtAvgMistag_ = tot_mistag/totalTagged;
+}
+
+LauParameter* LauFlavTag::findParameter(const TString& parName)
+{
+ for ( std::vector<LauParameter*>::iterator iter = params_.begin(); iter != params_.end(); ++iter ) {
+ if ((*iter)->name().Contains(parName)) {
+ return (*iter);
+ }
+ }
+ std::cerr << "ERROR in LauFlavTag::findParameter : Parameter \"" << parName << "\" not found." << std::endl;
+ return 0;
+}
+
+Double_t LauFlavTag::getOmega(const UInt_t iEvt) const
+{
+ return calib_p0_->unblindValue() + calib_p1_->unblindValue() * (evtMistagVals_[iEvt] - perEvtAvgMistag_);
+}
+
+std::vector<LauParameter*> LauFlavTag::getCalibParameters()
+{
+ std::vector<LauParameter*> calibParams;
+ calib_p0_ = this->findParameter("calib_p0");
+ calib_p1_ = this->findParameter("calib_p1");
+
+ // calib_p0_->value(perEvtAvgMistag_); // WOULD BE NICE TO DO THIS
+
+ calibParams.push_back( calib_p0_ );
+ calibParams.push_back( calib_p1_ );
+ if (calibParams.size() != 2) {
+ std::cout << "FATAL in LauFlavTag::setCalibParams : Expected two calibration parameters, received " << calibParams.size() << "." << std::endl;
+ exit(1);
+ }
+ return calibParams;
+
+}
diff --git a/src/LauTimeDepFitModel.cc b/src/LauTimeDepFitModel.cc
index ab1e0ea..8104d0b 100644
--- a/src/LauTimeDepFitModel.cc
+++ b/src/LauTimeDepFitModel.cc
@@ -1,2321 +1,2330 @@
// Copyright University of Warwick 2006 - 2014.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file LauTimeDepFitModel.cc
\brief File containing implementation of LauTimeDepFitModel class.
*/
#include <iostream>
#include <iomanip>
#include <fstream>
#include <map>
#include <vector>
#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"
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),
- flavTag_(flavTag),
+ flavTag_(flavTag),
nSigComp_(0),
nSigDPPar_(0),
nDecayTimePar_(0),
nExtraPdfPar_(0),
nNormPar_(0),
- nCalibPar_(0),
+ nCalibPar_(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), -3.0, 3.0, kFALSE),
cosPhiMix_("cosPhiMix", TMath::Cos(2.0*LauConstants::beta), -3.0, 3.0, kFALSE),
useSinCos_(kFALSE),
phiMixComplex_(TMath::Cos(-2.0*LauConstants::beta),TMath::Sin(-2.0*LauConstants::beta)),
signalDecayTimePdfs_(),
curEvtDecayTime_(0.0),
curEvtDecayTimeErr_(0.0),
sigExtraPdf_(),
sigFlavTagPdf_(),
bkgdFlavTagPdf_(),
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 );
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(LauParameter* 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;
}
bkgndEvents_[bkgndID] = nBkgndEvents;
bkgndEvents_[bkgndID]->name( nBkgndEvents->name()+"Events" );
Double_t value = nBkgndEvents->value();
bkgndEvents_[bkgndID]->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
bkgndAsym_[bkgndID] = new LauParameter(nBkgndEvents->name()+"Asym",0.0,-1.0,1.0,kTRUE);
}
void LauTimeDepFitModel::setNBkgndEvents(LauParameter* nBkgndEvents, LauParameter* 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] = nBkgndEvents;
bkgndEvents_[bkgndID]->name( nBkgndEvents->name()+"Events" );
Double_t value = nBkgndEvents->value();
bkgndEvents_[bkgndID]->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
bkgndAsym_[bkgndID] = bkgndAsym;
bkgndAsym_[bkgndID]->name( nBkgndEvents->name()+"Asym" );
bkgndAsym_[bkgndID]->range(-1.0,1.0);
}
void LauTimeDepFitModel::setSignalDtPdf(Int_t tagCat, LauDecayTimePdf* pdf)
{
if (!flavTag_->validTagCat(tagCat)) {
std::cerr<<"ERROR in LauTimeDepFitModel::setSignalDtPdf : Tagging category \""<<tagCat<<"\" not valid, not adding PDF."<<std::endl;
return;
}
if (pdf==0) {
std::cerr<<"ERROR in LauTimeDepFitModel::setSignalDtPdf : The PDF pointer is null, not adding it."<<std::endl;
return;
}
signalDecayTimePdfs_[tagCat] = pdf;
}
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 \""<<tagCat<<"\" not valid, not adding models."<<std::endl;
return;
}
if (pdf==0) {
std::cerr<<"ERROR in LauTimeDepFitModel::setSignalPdfs : The PDF pointer is null."<<std::endl;
return;
}
sigExtraPdf_[tagCat].push_back(pdf);
}
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."<<std::endl;
gSystem->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."<<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 calibration parameters
- this->setCalibParams();
+ // Handle flavour-tagging calibration parameters
+ this->setCalibParams();
// Check that we have the expected number of fit variables
const LauParameterPList& fitVars = this->fitPars();
if (fitVars.size() != (nSigDPPar_ + nDecayTimePar_ + nExtraPdfPar_ + nNormPar_ + nCalibPar_)) {
std::cerr<<"ERROR in LauTimeDepFitModel::initialise : Number of fit parameters not of expected size."<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
this->setExtraNtupleVars();
-
- //Flavour tagging
- flavTag_->initialise();
-
}
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"<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (sigModelB0_ == 0) {
std::cerr<<"ERROR in LauTimeDepFitModel::initialiseDPModels : B0 signal DP model doesn't exist"<<std::endl;
gSystem->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"<<std::endl;
sigModelB0bar_->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();
}
void LauTimeDepFitModel::calcInterferenceTermIntegrals()
{
const std::vector<LauDPPartialIntegralInfo*>& integralInfoListB0bar = sigModelB0bar_->getIntegralInfos();
const std::vector<LauDPPartialIntegralInfo*>& 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<Double_t>& fNormB0bar = sigModelB0bar_->getFNorm();
const std::vector<Double_t>& 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 \""<<compName<<"\"."<<std::endl;
return;
}
std::cerr<<"WARNING in LauTimeDepFitModel::setAmpCoeffSet : B0bar signal DP model doesn't contain component \""<<compName<<"\" but does contain the conjugate \""<<conjName<<"\", restting name to use the conjugate."<<std::endl;
TString tmp = compName;
compName = conjName;
conjName = tmp;
coeffSet->name( compName );
}
if ( conjugate ) {
if ( ! sigModelB0_->hasResonance(conjName) ) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<<conjName<<"\"."<<std::endl;
return;
}
} else {
if ( ! sigModelB0_->hasResonance(compName) ) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<<compName<<"\"."<<std::endl;
return;
}
}
// Do we already have it in our list of names?
for (std::vector<LauAbsCoeffSet*>::const_iterator iter=coeffPars_.begin(); iter!=coeffPars_.end(); ++iter) {
if ((*iter)->name() == compName) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : Have already set coefficients for \""<<compName<<"\"."<<std::endl;
return;
}
}
coeffSet->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 \""<<compName<<"\" to the fit model."<<std::endl;
}
void LauTimeDepFitModel::calcAsymmetries(Bool_t initValues)
{
// Calculate the CP asymmetries
// Also calculate the fit fraction asymmetries
for (UInt_t i = 0; i < nSigComp_; i++) {
acp_[i] = coeffPars_[i]->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_);
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_);
}
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..."<<std::endl;
fitVars.push_back(signalEvents_);
++nNormPar_;
} else {
std::cout<<"INFO in LauTimeDepFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..."<<std::endl;
}
// if not using the DP in the model we need an explicit signal asymmetry parameter
if (this->useDP() == kFALSE) {
fitVars.push_back(signalAsym_);
++nNormPar_;
}
- LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac();
+ // 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_;
}
}
void LauTimeDepFitModel::setCalibParams()
{
- std::vector<LauParameter *> calibParams = flavTag_->getCalibParameters();
+ std::vector<LauParameter *> calibParams = flavTag_->getCalibParameters();
- // If calib params aren't fixed, add them to the vector of floating parameters
+ // If calib params aren't fixed, add them to the vector of floating parameters
- LauParameter * p0 = calibParams[0];
- LauParameter * p1 = calibParams[1];
+ LauParameter * p0 = calibParams[0];
+ LauParameter * p1 = calibParams[1];
- if (!p0->fixed()) {
- LauParameterPList& fitVars = this->fitPars();
- nCalibPar_ += 1;
- fitVars.push_back(p0);
- }
+ if (!p0->fixed()) {
+ LauParameterPList& fitVars = this->fitPars();
+ nCalibPar_ += 1;
+ fitVars.push_back(p0);
+ }
- if (!p1->fixed()) {
- LauParameterPList& fitVars = this->fitPars();
- nCalibPar_ += 1;
- fitVars.push_back(p1);
- }
+ if (!p1->fixed()) {
+ LauParameterPList& fitVars = this->fitPars();
+ nCalibPar_ += 1;
+ fitVars.push_back(p1);
+ }
}
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: "<<fitFracB0bar_.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0bar_[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0bar_[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j = i; j < nSigComp_; j++) {
TString name = fitFracB0bar_[i][j].name();
name.Insert( name.Index("FitFrac"), "B0bar" );
fitFracB0bar_[i][j].name(name);
extraVars.push_back(fitFracB0bar_[i][j]);
}
}
fitFracB0_ = sigModelB0_->getFitFractions();
if (fitFracB0_.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0_.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0_[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0_[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j = i; j < nSigComp_; j++) {
TString name = fitFracB0_[i][j].name();
name.Insert( name.Index("FitFrac"), "B0" );
fitFracB0_[i][j].name(name);
extraVars.push_back(fitFracB0_[i][j]);
}
}
// Calculate the ACPs and FitFrac asymmetries
this->calcAsymmetries(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::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();
}
if (useSinCos_) {
cosPhiMix_.updatePull();
sinPhiMix_.updatePull();
} else {
this->checkMixingPhase();
}
// Update the pulls on all the extra PDFs' parameters
for (LauTagCatPdfListMap::iterator iter = sigExtraPdf_.begin(); iter != sigExtraPdf_.end(); ++iter) {
this->updateFitParameters(iter->second);
}
- LauTagCatParamMap signalTagCatFrac = flavTag_->getSignalTagCatFrac();
+ 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: "<<fitFracB0bar.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0bar[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0bar[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
LauParArray fitFracB0 = sigModelB0_->getFitFractions();
if (fitFracB0.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j(i); j<nSigComp_; ++j) {
fitFracB0bar_[i][j].value(fitFracB0bar[i][j].value());
fitFracB0_[i][j].value(fitFracB0[i][j].value());
}
}
meanEffB0bar_.value(sigModelB0bar_->getMeanEff().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); i<nSigComp_; ++i) {
for (UInt_t j(i); j<nSigComp_; ++j) {
extraVars.push_back(fitFracB0bar_[i][j]);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j(i); j<nSigComp_; ++j) {
extraVars.push_back(fitFracB0_[i][j]);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
extraVars.push_back(fitFracAsymm_[i]);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
extraVars.push_back(acp_[i]);
}
extraVars.push_back(meanEffB0bar_);
extraVars.push_back(meanEffB0_);
extraVars.push_back(DPRateB0bar_);
extraVars.push_back(DPRateB0_);
this->printFitFractions(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 "<<i<<" ("<<compName<<") = "<<fitFracB0bar_[i][i]<<std::endl;
}
output<<"B0bar overall DP rate (integral of matrix element squared) = "<<DPRateB0bar_<<std::endl;
output<<"B0bar average efficiency weighted by whole DP dynamics = "<<meanEffB0bar_<<std::endl;
// Then for the B0 sample
for (UInt_t i = 0; i < nSigComp_; i++) {
const TString compName(coeffPars_[i]->name());
const TString conjName(sigModelB0bar_->getConjResName(compName));
output<<"B0 FitFraction for component "<<i<<" ("<<conjName<<") = "<<fitFracB0_[i][i]<<std::endl;
}
output<<"B0 overall DP rate (integral of matrix element squared) = "<<DPRateB0_<<std::endl;
output<<"B0 average efficiency weighted by whole DP dynamics = "<<meanEffB0_<<std::endl;
}
void LauTimeDepFitModel::printAsymmetries(std::ostream& output)
{
for (UInt_t i = 0; i < nSigComp_; i++) {
const TString compName(coeffPars_[i]->name());
output<<"Fit Fraction asymmetry for component "<<i<<" ("<<compName<<") = "<<fitFracAsymm_[i]<<std::endl;
}
for (UInt_t i = 0; i < nSigComp_; i++) {
const TString compName(coeffPars_[i]->name());
output<<"ACP for component "<<i<<" ("<<compName<<") = "<<acp_[i].value()<<" +- "<<acp_[i].error()<<std::endl;
}
}
void LauTimeDepFitModel::writeOutTable(const TString& outputFile)
{
// Write out the results of the fit to a tex-readable table
std::ofstream fout(outputFile);
LauPrint print;
std::cout<<"INFO in LauTimeDepFitModel::writeOutTable : Writing out results of the fit to the tex file "<<outputFile<<std::endl;
if (this->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"<<std::endl;
fout<<"\\end{tabular}"<<std::endl<<std::endl;
// print the fit fractions in another
fout<<"\\begin{tabular}{|l|c|c|c|c|}"<<std::endl;
fout<<"\\hline"<<std::endl;
fout<<"Component & \\Bzb\\ Fit Fraction & \\Bz\\ Fit Fraction & Fit Fraction Asymmetry & $A_{\\CP}$ \\\\"<<std::endl;
fout<<"\\hline"<<std::endl;
Double_t fitFracSumB0bar(0.0);
Double_t fitFracSumB0(0.0);
for (UInt_t i = 0; i < nSigComp_; i++) {
TString resName = coeffPars_[i]->name();
resName = resName.ReplaceAll("_", "\\_");
fout<<resName<<" & $";
Double_t fitFracB0bar = fitFracB0bar_[i][i].value();
fitFracSumB0bar += fitFracB0bar;
print.printFormat(fout, fitFracB0bar);
fout << "$ & $" << std::endl;
Double_t fitFracB0 = fitFracB0_[i][i].value();
fitFracSumB0 += fitFracB0;
print.printFormat(fout, fitFracB0);
fout << "$ & $" << std::endl;
Double_t fitFracAsymm = fitFracAsymm_[i].value();
print.printFormat(fout, fitFracAsymm);
fout << "$ & $" << std::endl;
Double_t acp = acp_[i].value();
Double_t acpErr = acp_[i].error();
print.printFormat(fout, acp);
fout<<" \\pm ";
print.printFormat(fout, acpErr);
fout<<"$ \\\\"<<std::endl;
}
fout<<"\\hline"<<std::endl;
// Also print out sum of fit fractions
fout << "Fit Fraction Sum = & $";
print.printFormat(fout, fitFracSumB0bar);
fout << "$ & $";
print.printFormat(fout, fitFracSumB0);
fout << "$ & & \\\\" << std::endl;
fout << "\\hline \n\\hline" << std::endl;
fout << "DP rate = & $";
print.printFormat(fout, DPRateB0bar_.value());
fout << "$ & $";
print.printFormat(fout, DPRateB0_.value());
fout << "$ & & \\\\" << std::endl;
fout << "$< \\varepsilon > =$ & $";
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|}"<<std::endl;
fout<<"\\hline"<<std::endl;
fout<<"\\Extra Signal PDFs' Parameters: & \\\\"<<std::endl;
for (LauTagCatPdfListMap::const_iterator iter = sigExtraPdf_.begin(); iter != sigExtraPdf_.end(); ++iter) {
this->printFitParameters(iter->second, fout);
}
fout<<"\\hline \n\\end{tabular}"<<std::endl<<std::endl;
}
}
void LauTimeDepFitModel::checkInitFitParams()
{
// Update the number of signal events to be total-sum(background events)
this->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)..."<<std::endl;
for (UInt_t i = 0; i < nSigComp_; i++) {
coeffPars_[i]->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();
+ 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:"<<std::endl;
std::cout<<" : Signal asymmetry = "<<sigAsym<<" and number of signal events = "<<signalEvents_->genValue()<<std::endl;
return nEvtsGen;
}
Bool_t LauTimeDepFitModel::genExpt()
{
// Routine to generate toy Monte Carlo events according to the various models we have defined.
// Determine the number of events to generate for each hypothesis
LauGenInfo nEvts = this->eventsToGenerate();
Bool_t genOK(kTRUE);
Int_t evtNum(0);
// Loop over the hypotheses and generate the appropriate number of
// events for each one
for (LauGenInfo::const_iterator iter = nEvts.begin(); iter != nEvts.end(); ++iter) {
// find the category of events (e.g. signal)
const TString& evtCategory(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); iEvt<nEvtsGen; ++iEvt) {
this->setGenNtupleDoubleBranchValue( "evtWeight", evtWeight );
if (evtCategory == "signal") {
this->setGenNtupleIntegerBranchValue("genSig",1);
// 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;
}
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
}
// 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_);
// 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 "<<iEvt<<" out of "<<nEvtsGen<<" "<<evtCategory<<" events."<<std::endl;}
} // end of loop over events for given species, tagCat and tagFlv
if (!genOK) {
break;
}
} // end of loop over tagCats for the given species and tagFlv
if (!genOK) {
break;
}
} //end of loop over species and tagFlv.
if (this->useDP() && genOK) {
sigModelB0bar_->checkToyMC(kTRUE);
sigModelB0_->checkToyMC(kTRUE);
std::cout<<"aSqMaxSet = "<<aSqMaxSet_<<" and aSqMaxVar = "<<aSqMaxVar_<<std::endl;
// Get the fit fractions if they're to be written into the latex table
if (this->writeLatexTable()) {
LauParArray fitFracB0bar = sigModelB0bar_->getFitFractions();
if (fitFracB0bar.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0bar.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0bar[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0bar[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
LauParArray fitFracB0 = sigModelB0_->getFitFractions();
if (fitFracB0.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j(i); j<nSigComp_; ++j) {
fitFracB0bar_[i][j].value(fitFracB0bar[i][j].value());
fitFracB0_[i][j].value(fitFracB0[i][j].value());
}
}
meanEffB0bar_.value(sigModelB0bar_->getMeanEff().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();
}
}
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;
// 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; // B0 tag
} else {
curEvtTagFlv_ = -1; // B0bar tag
}
// Generate the mistag
// - For tagged events use random for now... to be fixed
// - Mistag = 0.5 for untagged events (matches LHCb data but doesn't matter)
if (curEvtTagCat_==0){
curEvtMistag_ = 0.5;
} else {
curEvtMistag_ = LauRandom::randomFun()->Gaus(0.45, 0.1);
while (curEvtMistag_ > 0.5) {
curEvtMistag_ = LauRandom::randomFun()->Gaus(0.45, 0.1);
}
}
// Calculate event quantities that depend only on the tagCat and tagFlv
Double_t qD(0.);
Double_t qDDo2(0.);
if(flavTag_->getUsePerEvtMistag()){
qD = curEvtTagFlv_*(1-2*curEvtMistag_);
}else{
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_);
// ...and check that the calculation went ok, otherwise loop again
if (decayTimePdf->state() != LauDecayTimePdf::Good) {
std::cout<<"decayTimePdf state is bad"<<std::endl;
++nGenLoop_;
continue;
}
// First get all the decay time terms
//Double_t dtExp = decayTimePdf->getExpTerm();
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;
// 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;
//std::cout << "dtCos * qD * aSqDif (dtSin * interTermIm) " << dtCos << " " << qD << " " << aSqDif << " " << dtSin << " " << interTermIm << std::endl;
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;
//std::cout<<"Cosh Cos Sin Sinh "<< coshTerm<< " " << cosTerm << " " << sinTerm << " " << sinhTerm << std::endl;
//std::cout << "Total Amplitude : " << ASq << std::endl;
//ASq /= decayTimePdf->getNormTerm();
ASq *= eff;
//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_<<std::endl;
} else if (aSqMaxVar_ > aSqMaxSet_) {
aSqMaxSet_ = 1.01 * aSqMaxVar_;
genOK = kFALSE;
std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Found a larger ASq value: setting aSqMaxSet_ to "<<aSqMaxSet_<<std::endl;
}
if (genOK) {
//Some variables, like Fisher or NN, might use m13Sq and m23Sq from the kinematics
//kinematicsB0bar_ is up to date, update kinematicsB0_
kinematicsB0_->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."<<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<TString>::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
- if (!this->doEMLFit()) {
- this->updateSigEvents();
- }
+ // 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.
- Double_t nTotEvts = this->eventsPerExpt();
- Double_t signalEvents = nTotEvts;
-
- LauTagCatParamMap signalTagCatFrac = flavTag_->getSignalTagCatFrac();
+ 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);
+ if ( ! signalEvents_->fixed() ) {
+ signalEvents_->value(signalEvents);
+ }
+ }
// tagging-category fractions for signal events
- flavTag_->setFirstTagCatFrac(signalTagCatFrac);
-
- signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts);
-
- if ( ! signalEvents_->fixed() ) {
- signalEvents_->value(signalEvents);
- }
+ 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();
- flavTag_->cacheInputFitVars(inputFitData);
-
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<Int_t>( fitdata_iter->second );
if ( cpEV == 1 ) {
cpEigenValue_ = CPEven;
} else if ( cpEV == -1 ) {
cpEigenValue_ = CPOdd;
} else {
std::cerr<<"WARNING in LauTimeDepFitModel::cacheInputFitVars : Unknown value: "<<cpEV<<" for CP eigenvalue, setting to CP-even"<<std::endl;
cpEigenValue_ = CPEven;
}
}
evtCPEigenVals_.push_back( cpEigenValue_ );
}
// We'll cache the DP amplitudes at the end because we'll
// append some points that the other PDFs won't deal with.
+ // Flavour tagging information
+ flavTag_->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->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 (this->useDP() == kTRUE) {
sigModelB0bar_->fillDataTree(*inputFitData);
sigModelB0_->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)
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());
}
- LauTagCatParamMap signalTagCatFrac = flavTag_->getSignalTagCatFrac();
+ // 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*signalTagCatFrac[curEvtTagCat_].unblindValue());
+ 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 ) {
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-2*omega);
+ Double_t omega = flavTag_->getOmega(iEvt);
+ qD = curEvtTagFlv_*(1.0-2.0*omega);
} else {
- LauTagCatParamMap dilution_ = flavTag_->getDilution();
- LauTagCatParamMap deltaDilution_ = flavTag_->getDeltaDilution();
+ // 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);
// 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);
// 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();
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;
ASq *= eff;
// 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){
- normTermCosh = decayTimePdf->getNormTermCosh();
+ normTermCosh = decayTimePdf->getNormTermExp();
norm = normTermIndep*normTermCosh;
}
if (decayTimePdf->getFuncType() == LauDecayTimePdf::ExpHypTrig){
- normTermCosh = decayTimePdf->getNormTermExp();
+ 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 );
}
}
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();
}
// TODO Add in the background components too
}
void LauTimeDepFitModel::updateCoeffs()
{
coeffsB0bar_.clear(); coeffsB0_.clear();
coeffsB0bar_.reserve(nSigComp_); coeffsB0_.reserve(nSigComp_);
for (UInt_t i = 0; i < nSigComp_; ++i) {
coeffsB0bar_.push_back(coeffPars_[i]->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 "<<tagCat<<"."<<std::endl;
return;
}
if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
std::cerr<<"WARNING in LauTimeDepFitModel::embedSignal : Conflicting options provided, will not reuse events at all."<<std::endl;
reuseEventsWithinExperiment = kFALSE;
}
signalTree_[tagCat] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
Bool_t dataOK = signalTree_[tagCat]->findBranches();
if (!dataOK) {
delete signalTree_[tagCat]; signalTree_[tagCat] = 0;
std::cerr<<"ERROR in LauTimeDepFitModel::embedSignal : Problem creating data tree for embedding."<<std::endl;
return;
}
reuseSignal_ = reuseEventsWithinEnsemble;
}
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");
// Store the DP likelihoods
if (this->useDP()) {
this->addSPlotNtupleDoubleBranch("sigDPLike");
}
// Store the likelihoods for each extra PDF
const LauPdfList* pdfList( &(sigExtraPdf_.begin()->second) );
this->addSPlotNtupleBranches(pdfList, "sig");
}
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<TString>::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<TString>::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."<<std::endl;
}
}
}
Double_t LauTimeDepFitModel::setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt)
{
// Store the PDF value for each variable in the list
Double_t totalLike(1.0);
Double_t extraLike(0.0);
if ( !extraPdfs ) {
return totalLike;
}
for (LauPdfList::iterator pdf_iter = extraPdfs->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<TString>::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<TString>::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."<<std::endl;
}
}
return totalLike;
}
LauSPlot::NameSet LauTimeDepFitModel::variableNames() const
{
LauSPlot::NameSet nameSet;
if (this->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<TString>::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();
}
return numbMap;
}
LauSPlot::NumbMap LauTimeDepFitModel::fixdSpeciesNames() const
{
LauSPlot::NumbMap numbMap;
if (signalEvents_->fixed() && this->doEMLFit()) {
numbMap["sig"] = signalEvents_->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<TString>::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] ) ) );
}
}
return twodimMap;
}
void LauTimeDepFitModel::storePerEvtLlhds()
{
std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Storing per-event likelihood values..."<<std::endl;
LauFitDataTree* inputFitData = this->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);
+ 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_);
} else {
sigTotalLike_ = 1.0;
}
// the signal PDF values
sigTotalLike_ *= this->setSPlotNtupleBranchValues(sigPdfs, "sig", iEvt);
// the total likelihoods
this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_);
// fill the tree
this->fillSPlotNtupleBranches();
}
std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Finished storing per-event likelihood values."<<std::endl;
}
void LauTimeDepFitModel::weightEvents( const TString& /*dataFileName*/, const TString& /*dataTreeName*/ )
{
std::cerr << "ERROR in LauTimeDepFitModel::weightEvents : Method not available for this fit model." << std::endl;
return;
}
void LauTimeDepFitModel::savePDFPlots(const TString& /*label*/)
{
}
void LauTimeDepFitModel::savePDFPlotsWave(const TString& /*label*/, const Int_t& /*spin*/)
{
}
void LauTimeDepFitModel::setSignalFlavTagPdfs(const Int_t tagCat, LauAbsPdf* pdf)
{
if (!flavTag_->validTagCat(tagCat)) {
std::cerr<<"ERROR in LauTimeDepFitModel::setSignalFlavTagPdfs : Tagging category \""<<tagCat<<"\" not valid, not adding models."<<std::endl;
return;
}
if (pdf==0) {
std::cerr<<"ERROR in LauTimeDepFitModel::setSignalFlavTagPdfs : The PDF pointer is null."<<std::endl;
return;
}
sigFlavTagPdf_[tagCat] = pdf;
}
void LauTimeDepFitModel::setBkgdFlavTagPdfs(const TString name, LauAbsPdf* pdf)
{
// TODO - when backgrounds are added - implement a check that "name" really is the name of a background category
if (pdf==0) {
std::cerr<<"ERROR in LauTimeDepFitModel::setBkgdFlavTagPdfs : The PDF pointer is null."<<std::endl;
return;
}
bkgdFlavTagPdf_[name] = pdf;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 2:30 PM (1 d, 10 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804745
Default Alt Text
(124 KB)

Event Timeline