diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.cc b/examples/ProgOpts/Test_Dpipi_ProgOpts.cc index 22f1508..adca8b8 100644 --- a/examples/ProgOpts/Test_Dpipi_ProgOpts.cc +++ b/examples/ProgOpts/Test_Dpipi_ProgOpts.cc @@ -1,151 +1,152 @@ #include #include "boost/program_options.hpp" #include "Test_Dpipi_ProgOpts.hh" namespace po = boost::program_options; void validate( boost::any& v, const std::vector& values, Command* /*target_type*/, int) { // Make sure no previous assignment to 'v' has been made po::validators::check_first_occurrence(v); // Extract the first string from 'values'. If there is more than // one string, it's an error, and exception will be thrown. const std::string& s { po::validators::get_single_string(values) }; if ( s == "gen" ) { v = boost::any( Command::Generate ); } else if ( s == "fit" ) { v = boost::any( Command::Fit ); } else if ( s == "simfit" ) { v = boost::any( Command::SimFit ); } else { throw po::validation_error(po::validation_error::invalid_option_value, "command", s, 3); } } void validate( boost::any& v, const std::vector& values, LauTimeDepFitModel::CPEigenvalue* /*target_type*/, int) { // Make sure no previous assignment to 'v' has been made po::validators::check_first_occurrence(v); // Extract the first string from 'values'. If there is more than // one string, it's an error, and exception will be thrown. const std::string& s { po::validators::get_single_string(values) }; if ( s == "QFS" ) { v = boost::any( LauTimeDepFitModel::CPEigenvalue::QFS ); } else if ( s == "CPEven" ) { v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPEven ); } else if ( s == "CPOdd" ) { v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPOdd ); } else { throw po::validation_error(po::validation_error::invalid_option_value); } } namespace LauDecayTime { void validate( boost::any& v, const std::vector& values, EfficiencyMethod* /*target_type*/, int) { // Make sure no previous assignment to 'v' has been made po::validators::check_first_occurrence(v); // Extract the first string from 'values'. If there is more than // one string, it's an error, and exception will be thrown. const std::string& s { po::validators::get_single_string(values) }; if ( s == "uniform" || s == "flat" ) { v = boost::any( EfficiencyMethod::Uniform ); } else if ( s == "binned" || s == "hist" ) { v = boost::any( EfficiencyMethod::Binned ); } else if ( s == "spline" ) { v = boost::any( EfficiencyMethod::Spline ); } else { throw po::validation_error(po::validation_error::invalid_option_value); } } }; TestDpipi_ProgramSettings::TestDpipi_ProgramSettings(const int argc, const char ** argv) { po::options_description main_desc{"Main options"}; main_desc.add_options() ("command", po::value(&command)->required(), "main command: gen, fit, or simfit") ; po::positional_options_description p; p.add("command", 1); po::options_description common_desc{"Common options"}; common_desc.add_options() ("help", "produce help message") ("dtype", po::value(&dType)->default_value(LauTimeDepFitModel::CPEigenvalue::QFS,"QFS"), "type of D decay: QFS, CPOdd, or CPEven") ("dta-model", po::value(&timeEffModel)->default_value(LauDecayTime::EfficiencyMethod::Uniform,"uniform"), "decay-time acceptance model: uniform/flat, binned/hist, spline") ("dtr", po::value(&timeResolution)->default_value(kTRUE), "enable/disable decay-time resolution") ("dtr-perevent", po::value(&perEventTimeErr)->default_value(kFALSE), "enable/disable use of per-event decay-time error (requires decay-time resolution to be enabled to take effect)") ("seed", po::value(&RNGseed)->default_value(0), "set the seed for the RNG; if not set, the time is used to generate a seed") ("dir", po::value(&directory)->default_value("",""), "set the directory used to find the nTuples within the root file. Defaults to no directory") ; po::options_description gen_desc{"Generation options"}; gen_desc.add_options() ("firstExptGen", po::value(&firstExptGen)->default_value(0), "first experiment to generate") ("nExptGen", po::value(&nExptGen)->default_value(1), "number of experiments to generate") ; po::options_description fit_desc{"Fitting options"}; fit_desc.add_options() ("firstExpt", po::value(&firstExptFit)->default_value(0), "first experiment to fit") ("nExpt", po::value(&nExptFit)->default_value(1), "number of experiments to fit") ("iFit", po::value(&iFit)->default_value(0), "fit ID - used to distinguish multiple fits to a given dataset (used to disentangle multiple minima)") ("fixTau", po::value(&fixLifetime)->default_value(kTRUE), "enable/disable floating of B lifetime parameter") ("fixDm", po::value(&fixDeltaM)->default_value(kTRUE), "enable/disable floating of B mixing frequency parameter") ("fixPhiMix", po::value(&fixPhiMix)->default_value(kFALSE), "enable/disable floating of B mixing phase parameter(s)") ("fixSplineKnots", po::value(&fixSplineKnots)->default_value(kFALSE), "enable/disable floating of the decay-time-acceptance spline") ("useSinCos", po::value(&useSinCos)->default_value(kTRUE), "enable/disable using sine and cosine of B mixing phase as the floating parameters rather than the phase itself") ("useAveDeltaCalibVals", po::value(&useAveDeltaCalibVals)->default_value(kTRUE), "enable/disable fitting calib values as their average +/- delta values rather than having separate values for B and Bbar") ("addTaggers", po::value(&addTaggers)->default_value(kTRUE), "add/omit taggers") ("floatCalibPars", po::value(&floatCalibPars)->default_value(kTRUE), "enable/disable floating of the FT calibration parameters") + ("floatYields", po::value(&floatYields)->default_value(kFALSE), "enable/disable floating of the yields of each component - if enabled/disabled implies inclusion/exclusion of non-DP PDFs in/from the likelihood") ("fitBack", po::value(&fitBack)->default_value(kFALSE), "enable/disable refitting some generated data from a model") ; po::options_description simfit_desc{"Simultaneous fitting options"}; simfit_desc.add_options() ("port", po::value(&port)->default_value(0), "port number on which sim fit coordinator is listening") ; po::options_description all_desc; all_desc.add(main_desc).add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc); po::variables_map vm; try { po::store(po::command_line_parser(argc, argv). options(all_desc). positional(p). run(), vm); po::notify(vm); } catch ( boost::wrapexcept& e ) { std::cout << argv[0] << " requires a command, one of 'gen', 'fit', or 'simfit'\n" << "Append --help to those commands to see help on the related options" << std::endl; parsedOK = kFALSE; return; } catch ( po::validation_error& e ) { std::cerr << e.what() << std::endl; parsedOK = kFALSE; return; } if ( vm.count("help") ) { po::options_description help_desc; help_desc.add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc); std::cout << help_desc << std::endl; helpRequested = kTRUE; return; } } diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.hh b/examples/ProgOpts/Test_Dpipi_ProgOpts.hh index 138d19c..e5b4d9e 100644 --- a/examples/ProgOpts/Test_Dpipi_ProgOpts.hh +++ b/examples/ProgOpts/Test_Dpipi_ProgOpts.hh @@ -1,45 +1,46 @@ #ifndef TEST_DPIPI_PROGOPTS #define TEST_DPIPI_PROGOPTS #include "Rtypes.h" #include "LauDecayTime.hh" #include "LauTimeDepFitModel.hh" #include "Command.hh" #include struct TestDpipi_ProgramSettings { TestDpipi_ProgramSettings(const int argc, const char ** argv); Bool_t parsedOK{kTRUE}; Bool_t helpRequested{kFALSE}; Command command{Command::Generate}; LauTimeDepFitModel::CPEigenvalue dType{LauTimeDepFitModel::CPEigenvalue::QFS}; LauDecayTime::EfficiencyMethod timeEffModel{LauDecayTime::EfficiencyMethod::Uniform}; UInt_t firstExptGen{0}; UInt_t nExptGen{1}; UInt_t firstExptFit{0}; UInt_t nExptFit{1}; UInt_t iFit{0}; UInt_t port{0}; UInt_t RNGseed{0}; Bool_t timeResolution{kTRUE}; Bool_t perEventTimeErr{kFALSE}; Bool_t fixLifetime{kTRUE}; Bool_t fixDeltaM{kTRUE}; Bool_t fixPhiMix{kFALSE}; Bool_t fixSplineKnots{kFALSE}; Bool_t useSinCos{kTRUE}; Bool_t useAveDeltaCalibVals{kTRUE}; Bool_t addTaggers{kTRUE}; Bool_t floatCalibPars{kTRUE}; + Bool_t floatYields{kFALSE}; std::string directory{""}; Bool_t fitBack{kFALSE}; }; #endif diff --git a/inc/LauFlavTag.hh b/inc/LauFlavTag.hh index f0bb278..5e03dac 100644 --- a/inc/LauFlavTag.hh +++ b/inc/LauFlavTag.hh @@ -1,536 +1,539 @@ /* Copyright 2017 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file 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 #include #include "TString.h" #include "LauParameter.hh" class LauAbsPdf; class LauFitDataTree; class LauFlavTag final { public: //! Define sign convention for B and Bbar flavours enum Flavour { Bbar = -1, //< Bbar flavour Unknown = 0, //< Unknown flavour B = +1 //< B flavour }; //! Define different types of background to control the behaviour for each source // TODO Might want to move this somewhere more general later enum class BkgndType { Combinatorial, //< combinatorial background FlavourSpecific, //< flavour-specific decay (i.e. the flavour of the parent can be determined from the decay products), e.g. B0 -> K+ pi- pi0 SelfConjugate, //< decays where both B and Bbar can decay to a single self-conjugate final state, e.g. B0 -> pi+ pi- K_S0 NonSelfConjugate //< decays where both B and Bbar can decay to either of two final states that are conjugates of each other, e.g. B0 -> K+ pi- K_S0 and B0 -> K- pi+ K_S0 }; //! Constructor /*! \param [in] useAveDelta use average and delta variables for tagging calibration and efficiency \param [in] useEtaPrime use eta prime rather the eta as the mistag throughout \param [in] bkgndInfo map containing names and types of the background sources (if applicable) */ LauFlavTag(const Bool_t useAveDelta = kFALSE, const Bool_t useEtaPrime = kFALSE, const std::map bkgndInfo={}); //! Initialise // TODO is this needed? Commented for the moment (here and where called in LauTimeDepFitModel) //void initialise(); // TODO - need to decide which functions need to be public (interface) and which should be private (implementation details) // - improve/extend Doxygen comments //! Change the dilutions, delta dilutions and tagCatFrac for signal if needed /*! \param [in] name the name of the tagger \param [in] tagVarName the tagging variable name of the tagger in the ntuple \param [in] mistagVarName the associated mistag variable name of the same tagger in the ntuple \param [in] etapdf the mistag distribution for the tagger \param [in] tagEff tagging efficiency - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag \param [in] calib_p0 calibration parameter p0 - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag \param [in] calib_p1 calibration parameter p1 - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag */ void addTagger(const TString& name, const TString& tagVarName, const TString& mistagVarName, LauAbsPdf* etapdf, const std::pair tagEff, const std::pair calib_p0, const std::pair calib_p1); //! Change the dilutions, delta dilutions and tagCatFrac for signal if needed /*! \param [in] name the name of the tagger \param [in] tagVarName the tagging variable name of the tagger in the ntuple \param [in] mistagVarName the associated mistag variable name of the same tagger in the ntuple \param [in] etapdf the mistag distribution for the tagger \param [in] tagEff tagging efficiency histograms - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag \param [in] calib_p0 calibration parameter p0 - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag \param [in] calib_p1 calibration parameter p1 - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag */ void addTagger(const TString& name, const TString& tagVarName, const TString& mistagVarName, LauAbsPdf* etapdf, const std::pair tagEff, const std::pair calib_p0, const std::pair calib_p1); //! Read in the input fit data variables /*! \param [in] inputFitData the data source \param [in] decayTimeVarName the name of the decay time variable within the input data (used if the tagging efficiencies depend on decay time) */ void cacheInputFitVars(LauFitDataTree* inputFitData, const TString& decayTimeVarName=""); //! Generate values for the signal tag decisions and mis-tag probabilities /*! \param [in] trueTagFlv the true tag flavour \param [in] curEvtDecayTime the generated decay time value (used if the tagging efficiencies depend on decay time) */ void generateEventInfo(const Flavour trueTagFlv, const Double_t curEvtDecayTime); //! Generate values for the background tag decisions and mis-tag probabilities /*! \param [in] bkgndID the background category ID for which to generate \param [in] trueTagFlv the true tag flavour \param [in] trueDecayFlv the true decay flavour (if known) \param [in] curEvtDecayTime the generated decay time value (used if the tagging efficiencies depend on decay time) */ void generateBkgndEventInfo(const std::size_t bkgndID, const Flavour trueTagFlv, const Flavour trueDecayFlv, const Double_t curEvtDecayTime); //! Retrieve the cached info for a given event /*! \param [in] iEvt the event to retrieve */ void updateEventInfo(const std::size_t iEvt); //! Retrieve the name of the true tag variable const TString& getTrueTagVarName() const {return trueTagVarName_;}; //! Retrieve the name of the decay flavour variable const TString& getDecayFlvVarName() const {return decayFlvVarName_;}; //! Retrieve the names of the tag decision variables for each tagger const std::vector& getTagVarNames() const {return tagVarNames_;}; //! Retrieve the names of the mis-tag probability variables for each tagger const std::vector& getMistagVarNames() const {return mistagVarNames_;}; //! Retrieve the current value of the true tag variable Flavour getCurEvtTrueTagFlv() const {return curEvtTrueTagFlv_;}; //! Retrieve the current value of the decay flavour variable Flavour getCurEvtDecayFlv() const {return curEvtDecayFlv_;}; //! Retrieve the current values of the tag decision variables for each tagger const std::vector& getCurEvtTagFlv() const {return curEvtTagFlv_;}; //! Retrieve the current values of the mis-tag probability variables for each tagger const std::vector& getCurEvtMistag() const {return curEvtMistag_;}; //! Retrieve the number of taggers std::size_t getNTaggers() const {return tagVarNames_.size();} //! Get vector of calibration p0 for B0 parameters for each tagger std::vector getCalibP0B0(){return calib_p0_B0_;}; //! Get vector of calibration p0 for B0bar parameters for each tagger std::vector getCalibP0B0bar(){return calib_p0_B0bar_;}; //! Get vector of calibration p1 for B0 parameters for each tagger std::vector getCalibP1B0(){return calib_p1_B0_;}; //! Get vector of calibration p1 for B0bar parameters for each tagger std::vector getCalibP1B0bar(){return calib_p1_B0bar_;}; //! Get vector of calibration p0 average parameters for each tagging category std::vector getCalibP0Ave(){return calib_p0_ave_;}; //! Get vector of calibration p0 difference parameters for each tagging category std::vector getCalibP0Delta(){return calib_p0_delta_;}; //! Get vector of calibration p1 average parameters for each tagging category std::vector getCalibP1Ave(){return calib_p1_ave_;}; //! Get vector of calibration p1 difference parameters for each tagging category std::vector getCalibP1Delta(){return calib_p1_delta_;}; //! Get vector of tagging efficiency for B0 parameters for each tagging category std::vector getTagEffB0(){return tagEff_B0_;}; //! Get vector of tagging efficiency for B0bar parameters for each tagging category std::vector getTagEffB0bar(){return tagEff_B0bar_;}; //! Get vector of tagging efficiency average parameters for each tagging category std::vector getTagEffAve(){return tagEff_ave_;}; //! Get vector of tagging efficiency difference parameters for each tagging category std::vector getTagEffDelta(){return tagEff_delta_;}; //! Get 2D vector of background tagging efficiency for B0 parameters for each tagger (inner vec) and background source (outer vec) auto getTagEffBkgndB0(){return tagEffBkgnd_B0_;}; //! Get 2D vector of background tagging efficiency for B0bar parameters for each tagger (inner vec) and background source (outer vec) auto getTagEffBkgndB0bar(){return tagEffBkgnd_B0bar_;}; //! Get 2D vector of background tagging efficiency average parameters for each tagger (inner vec) and background source (outer vec) auto getTagEffBkgndAve(){return tagEffBkgnd_ave_;}; //! Get 2D vector of background tagging efficiency difference parameters for each tagger (inner vec) and background source (outer vec) auto getTagEffBkgndDelta(){return tagEffBkgnd_delta_;}; //! Set background parameters for a given background and given tagger /*! Set the eta (mis-tag probability) distribution and tagging efficiency constants for B and Bbar \param [in] bkgndName background category name \param [in] taggerName name of the tagger \param [in] etaPdf the eta PDF itself \param [in] tagEff the tagging efficiency parameters */ void setBkgndParams(const TString& bkgndName, const TString& taggerName, LauAbsPdf* etaPdf, std::pair tagEff); //! Set background parameters for a given background and given tagger /*! Set the eta (mis-tag probability) distribution and tagging efficiencies (as a function of decay time) for B and Bbar \param [in] bkgndName background category name \param [in] taggerName name of the tagger \param [in] etaPdf the eta PDF itself \param [in] tagEff the tagging efficiency histograms */ void setBkgndParams(const TString& bkgndName, const TString& taggerName, LauAbsPdf* etaPdf, std::pair tagEff); //! Set background parameters for a given background and given tagger /*! Set the eta (mis-tag probability) distribution and tagging efficiency constants for B and Bbar This version provides different parameters for each decay flavour (used for Combinatorial background in QFS modes) \param [in] bkgndName background category name \param [in] taggerName name of the tagger \param [in] etaPdfB the eta PDF for decay flavour = B \param [in] tagEffB the tagging efficiency parameters for decay flavour = B \param [in] etaPdfBbar the eta PDF for decay flavour = Bbar \param [in] tagEffBbar the tagging efficiency parameters for decay flavour = Bbar */ void setBkgndParams(const TString& bkgndName, const TString& taggerName, LauAbsPdf* etaPdfB, std::pair tagEffB, LauAbsPdf* etaPdfBbar, std::pair tagEffBbar); //! Set background parameters for a given background and given tagger /*! Set the eta (mis-tag probability) distribution and tagging efficiencies (as a function of decay time) for B and Bbar This version provides different parameters for each decay flavour (used for Combinatorial background in QFS modes) \param [in] bkgndName background category name \param [in] taggerName name of the tagger \param [in] etaPdfB the eta PDF for decay flavour = B \param [in] tagEffB the tagging efficiency histograms for decay flavour = B \param [in] etaPdfBbar the eta PDF for decay flavour = Bbar \param [in] tagEffBbar the tagging efficiency histograms for decay flavour = Bbar */ void setBkgndParams(const TString& bkgndName, const TString& taggerName, LauAbsPdf* etaPdfB, std::pair tagEffB, LauAbsPdf* etaPdfBbar, std::pair tagEffBbar); //! Returns little omega (calibrated mistag) /*! \param [in] taggerID index of the tagger in the taggers vectors \param [in] flag choose to calculate omega or omegabar (corrsonding to B or Bbar) */ Double_t getLittleOmega(const std::size_t taggerID, const Flavour flag) const; //! Capital Omega for signal decays /*! \param [in] taggerID index of the tagger in the taggers vectors \param [in] flag choose to calculate Omega or Omegabar (corrsonding to B or Bbar) */ Double_t getCapitalOmega(const std::size_t taggerID, const Flavour flag) const; //! Returns little omega (calibrated mistag) for backgrounds /*! \param [in] bkgndID index of the background vector \param [in] taggerID index of the tagger in the taggers vectors \param [in] flag choose to calculate omega or omegabar (corrsonding to B or Bbar) \param [in] decayFlv the decay flavour of the B */ Double_t getLittleOmegaBkgnd(const std::size_t bkgndID, const std::size_t taggerID, const Flavour flag, const Flavour decayFlv) const; //! Capital Omega for backgrounds /*! \param [in] bkgndID index of the background vector \param [in] taggerID index of the tagger in the taggers vectors \param [in] flag choose to calculate Omega or Omegabar (corrsonding to B or Bbar) \param [in] decayFlv the decay flavour of the B */ Double_t getCapitalOmegaBkgnd(const std::size_t bkgndID, const std::size_t taggerID, const Flavour flag, const Flavour decayFlv) const; //! Get the generated value of signal eta for the given tagger /*! \param [in] taggerID index of the tagger in the taggers vectors */ Double_t getEtaGen(const std::size_t taggerID); //! Get the generated value of background eta for the given tagger and background category /*! \param [in] taggerID index of the tagger in the taggers vectors \param [in] bkgndID index of the background in the backgrounds vectors */ Double_t getEtaGenBkgnd(const std::size_t taggerID, const std::size_t bkgndID, const Flavour decayFlv); //! Return the Boolean controlling if we use the alternative tagging calibration and efficiency parameters Bool_t getUseAveDelta() const {return useAveDelta_;}; //! Specify the name of the true tag variable /*! \param [in] trueTagVarName the name of the true tag variable */ void setTrueTagVarName(TString trueTagVarName); //! Specify the name of the decay flavour variable /*! \param [in] decayFlvVarName the name of the decay flavour variable */ void setDecayFlvVarName(TString decayFlvVarName); //! Gaussian constraints for P0 parameters for a given tagger /*! \param [in] taggerName name of the tagger \param [in] constraint1 the (mean, sigma) for the particle or average parameter \param [in] constraint2 the (mean, sigma) for the antiparticle or delta parameter */ void addP0GaussianConstraints(const TString& taggerName, const std::pair constraint1, const std::pair constraint2); //! Gaussian constraints for P1 parameters for a given tagger /*! \param [in] taggerName name of the tagger \param [in] constraint1 the (mean, sigma) for the particle or average parameter \param [in] constraint2 the (mean, sigma) for the antiparticle or delta parameter */ void addP1GaussianConstraints(const TString& taggerName, const std::pair constraint1, const std::pair constraint2); //! Gaussian constraints for tagging efficiency parameters for a given tagger /*! \param [in] taggerName name of the tagger \param [in] constraint1 the (mean, sigma) for the particle or average parameter \param [in] constraint2 the (mean, sigma) for the antiparticle or delta parameter */ void addTagEffGaussianConstraints(const TString& taggerName, const std::pair constraint1, const std::pair constraint2); //! Retrieve the names of the background categories std::vector getBkgndNames(); //! Retrieve the types of the background categories std::vector getBkgndTypes(){return bkgndTypes_;} //! Float the P0 calibration parameters for B tags /*! \param [in] taggerName name of the tagger, defaults to empty to float for all taggers */ void floatCalibParP0B0(const TString& taggerName = ""); //! Float the P1 calibration parameters for B tags /*! \param [in] taggerName name of the tagger, defaults to empty to float for all taggers */ void floatCalibParP1B0(const TString& taggerName = ""); //! Float the P0 calibration parameters for Bbar tags /*! \param [in] taggerName name of the tagger, defaults to empty to float for all taggers */ void floatCalibParP0B0bar(const TString& taggerName = ""); //! Float the P1 calibration parameters for Bbar tags /*! \param [in] taggerName name of the tagger, defaults to empty to float for all taggers */ void floatCalibParP1B0bar(const TString& taggerName = ""); //! Float the P0 average calibration parameters /*! \param [in] taggerName name of the tagger, defaults to empty to float for all taggers */ void floatCalibParP0Ave(const TString& taggerName = ""); //! Float the P0 delta calibration parameters /*! \param [in] taggerName name of the tagger, defaults to empty to float for all taggers */ void floatCalibParP0Delta(const TString& taggerName = ""); //! Float the P1 average calibration parameters /*! \param [in] taggerName name of the tagger, defaults to empty to float for all taggers */ void floatCalibParP1Ave(const TString& taggerName = ""); //! Float the P1 delta calibration parameters /*! \param [in] taggerName name of the tagger, defaults to empty to float for all taggers */ void floatCalibParP1Delta(const TString& taggerName = ""); //! Float all calibration parameters void floatAllCalibPars(); + //! Update pulls on all floating parameters + void updatePulls(); + private: //! Internal function to extend vectors void extendVectors(const TString& tagVarName, const TString& mistagVarName); //! Internal function to setup Calib parameters void setupCalibParams(const TString& taggerName, const std::pair calib_p0, const std::pair calib_p1); //! Internal function to find parameters in per-tagger vectors LauParameter* findParameter( const TString& taggerName, std::vector& parameters ); //! Flag to use alternative calibration parameters const Bool_t useAveDelta_; //! Flag to use eta prime not eta for the mistag const Bool_t useEtaPrime_; //! Map to link tagger names to their index in all the vectors std::map taggerIndex_; //! Map to link background names to their index in all the vectors std::map bkgndIndex_; //! Flavour tagging variable name std::vector tagVarNames_; //! Per event mistag variable name std::vector mistagVarNames_; //! True tag variable name for normalisation decays TString trueTagVarName_; //! Decay flavour tag variable name for normalisation decays TString decayFlvVarName_; //! Vector of background types std::vector bkgndTypes_; //! Vector of flags indicating if the background parameters depend on the decay flavour std::vector bkgndDecayFlvDep_; //! Vector of flavour tags for each event std::vector> evtTagFlv_; //! Flavour tag for current event std::vector curEvtTagFlv_; //! Vector of mistags for each event std::vector> evtMistag_; //! Per event mistag for current event std::vector curEvtMistag_; //! Vector of true tags for each event std::vector evtTrueTagFlv_; //! Vector of decay tags for each event std::vector evtDecayFlv_; //! True tag from normalisation mode for current event Flavour curEvtTrueTagFlv_{Flavour::Unknown}; //! True tag from normalisation mode for current event Flavour curEvtDecayFlv_{Flavour::Unknown}; //! Average background mistag value (eta hat) std::vector>> avgMistagBkgnd_; //! Per-event average mistag value (eta hat) std::vector perEvtAvgMistag_; //! Decay time values for each event std::vector evtDecayTime_; //! Decay time value of the current event Double_t curEvtDecayTime_; //! Calibration parameters p0 for B0 std::vector calib_p0_B0_; //! Calibration parameters p0 for B0bar std::vector calib_p0_B0bar_; //! Calibration parameters p1 for B0 std::vector calib_p1_B0_; //! Calibration parameters p1 for B0bar std::vector calib_p1_B0bar_; //! Alternative calibration parameters p0 average std::vector calib_p0_ave_; //! Alternative calibration parameters p0 difference std::vector calib_p0_delta_; //! Alternative calibration parameters p1 average std::vector calib_p1_ave_; //! Alternative calibration parameters p1 difference std::vector calib_p1_delta_; //! Tagging efficiency parameters for B0 std::vector tagEff_B0_; //! Tagging efficiency parameters for B0bar std::vector tagEff_B0bar_; //! Tagging efficiency parameters average of B0 and B0bar std::vector tagEff_ave_; //! Tagging efficiency parameters difference between B0 and B0bar std::vector tagEff_delta_; //! Tagging efficiency histograms for B0 std::vector tagEff_hist_B0_; //! Tagging efficiency histograms for B0bar std::vector tagEff_hist_B0bar_; //! Tagging efficiency histograms average of B0 and B0bar std::vector tagEff_hist_ave_; //! Tagging efficiency histograms difference between B0 and B0bar std::vector tagEff_hist_delta_; //! Background tagging efficiency parameters for B0 std::vector>> tagEffBkgnd_B0_; //! Background tagging efficiency parameters for B0bar std::vector>> tagEffBkgnd_B0bar_; //! Background tagging efficiency parameters average of B0 and B0bar std::vector>> tagEffBkgnd_ave_; //! Background tagging efficiency parameters difference between B0 and B0bar std::vector>> tagEffBkgnd_delta_; //! Background tagging efficiency histograms for B0 std::vector>> tagEffBkgnd_hist_B0_; //! Background tagging efficiency histograms for B0bar std::vector>> tagEffBkgnd_hist_B0bar_; //! Background tagging efficiency histograms average of B0 and B0bar std::vector>> tagEffBkgnd_hist_ave_; //! Background tagging efficiency histograms difference between B0 and B0bar std::vector>> tagEffBkgnd_hist_delta_; //! Eta PDFs std::vector etaPdfs_; //! Eta PDFs for backgrounds per tagger (inner vec) and per background source (outer vec) std::vector>> etaBkgndPdfs_; ClassDef(LauFlavTag,0) // Flavour tagging set up }; //! output stream operator /*! \param [in,out] stream the output stream to which the background type is to be printed \param [in] bkgndType the type to be written out return the output stream */ std::ostream& operator<<( std::ostream& stream, const LauFlavTag::BkgndType bkgndType ); #endif diff --git a/inc/LauTimeDepFitModel.hh b/inc/LauTimeDepFitModel.hh index 9f42bde..d2f70d3 100644 --- a/inc/LauTimeDepFitModel.hh +++ b/inc/LauTimeDepFitModel.hh @@ -1,723 +1,725 @@ /* Copyright 2006 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauTimeDepFitModel.hh \brief File containing declaration of LauTimeDepFitModel class. */ /*! \class LauTimeDepFitModel \brief Class for defining a time-dependent fit model. LauTimeDepFitModel is a class that allows the user to define a three-body Dalitz plot according to the isobar model, i.e. defining a set of resonances that have complex amplitudes that can interfere with each other. It extends the LauSimpleFitModel and LauCPFitModel models in that it allows the fitting of time-dependent particle/antiparticle decays to flavour-conjugate Dalitz plots, including their interference through mixing. */ #ifndef LAU_TIMEDEP_FIT_MODEL #define LAU_TIMEDEP_FIT_MODEL #include #include #include #include #include "TString.h" #include "TStopwatch.h" #include "TSystem.h" #include "LauAbsFitModel.hh" #include "LauComplex.hh" #include "LauConstants.hh" #include "LauEmbeddedData.hh" #include "LauFlavTag.hh" #include "LauParameter.hh" class LauAbsBkgndDPModel; class LauAbsCoeffSet; class LauAbsPdf; class LauDecayTimePdf; class LauIsobarDynamics; class LauKinematics; class LauScfMap; class LauTimeDepFitModel : public LauAbsFitModel { public: //! Possible CP eigenvalues (the intrinsic CP of the final state particles) enum CPEigenvalue { CPOdd = -1, /*!< CP odd final state */ QFS = 0, /*!< Quasi Flavour Specific final state */ CPEven = 1 /*!< CP even final state */ }; //! Constructor /*! \param [in] modelB0bar DP model for the antiparticle \param [in] modelB0 DP model for the particle \param [in] flavTag flavour tagging information */ LauTimeDepFitModel(LauIsobarDynamics* modelB0bar, LauIsobarDynamics* modelB0, LauFlavTag* flavTag); //! Destructor virtual ~LauTimeDepFitModel(); //! Set the signal event yield /*! \param [in] nSigEvents contains the signal yield and option to fix it */ virtual void setNSigEvents(LauParameter* nSigEvents); + //TODO what is this asymmetry? //! Set the signal event yield and asymmetry /*! \param [in] nSigEvents contains the signal yield and option to fix it \param [in] sigAsym contains the signal asymmetry and option to fix it */ virtual void setNSigEvents(LauParameter* nSigEvents, LauParameter* sigAsym); //! Set the number of background events /*! The name of the parameter must be that of the corresponding background category (so that it can be correctly assigned) \param [in] nBkgndEvents contains the name, yield and option to fix the yield of the background */ virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents); + //TODO what is this asymmetry? //! Set the background event yield and asymmetry /*! The name of the parameters must be that of the corresponding background category (so that they can be correctly assigned) \param [in] nBkgndEvents contains the background yield and option to fix it \param [in] bkgndAsym contains the background asymmetry and option to fix it */ virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym); //! Set the background DP models (null pointer for BbarModel implies same model for both) /*! \param [in] bkgndClass the name of the background class \param [in] BModel the DP model of the background for B (particle) decays \param [in] BbarModel the DP model of the background for Bbar (anti-particle) decays */ void setBkgndDPModels(const TString& bkgndClass, LauAbsBkgndDPModel* BModel, LauAbsBkgndDPModel* BbarModel); //! Switch on/off storage of amplitude info in generated ntuple void storeGenAmpInfo(Bool_t storeInfo) { storeGenAmpInfo_ = storeInfo; } //! Set CP eigenvalue /*! The CP eigenvalue can be supplied on an event-by-event basis, e.g. if the data contains daughters that are D0 mesons that can decay to either K+ K- (CP even) or KS pi0 (CP odd). This method allows you to set the default value that should be used if the data does not contain this information as well as the name of the variable in the data that will specify this information. If completely unspecified all events will be assumed to be CP even. \param defaultCPEV the default for the eigenvalue \param cpevVarName the variable name in the data tree that specifies the CP eigenvalue */ inline void setCPEigenvalue( const CPEigenvalue defaultCPEV, const TString& cpevVarName = "" ) { cpEigenValue_ = defaultCPEV; cpevVarName_ = cpevVarName; } //! Set the DP amplitude coefficients /*! \param [in] coeffSet the set of coefficients */ void setAmpCoeffSet(LauAbsCoeffSet* coeffSet); //! Set the signal decay time PDF /*! \param [in] pdf the signal decay time PDF */ void setSignalDtPdf(LauDecayTimePdf* pdf); //! Set the decay time PDF for a given background class /*! \param [in] bkgndClass the name of the background class \param [in] pdf the background decay time PDF */ void setBkgndDtPdf(const TString& bkgndClass, LauDecayTimePdf* pdf); //! Set the signal PDF for a given variable /*! \param [in] pdf the PDF to be added to the signal model */ void setSignalPdfs(LauAbsPdf* pdf); //! Set the background PDF for a given variable and background class /*! \param [in] bkgndClass the name of the background class \param [in] pdf the PDF to be added to the background model */ void setBkgndPdf(const TString& bkgndClass, LauAbsPdf* pdf); //! Embed full simulation events for the signal, rather than generating toy from the PDFs /*! \param [in] fileName the name of the file containing the events \param [in] treeName the name of the tree \param [in] reuseEventsWithinEnsemble sample with replacement but only replace events once each experiment has been generated \param [in] reuseEventsWithinExperiment sample with immediate replacement \param [in] useReweighting perform an accept/reject routine using the configured signal amplitude model based on the MC-truth DP coordinate */ void embedSignal(const TString& fileName, const TString& treeName, const Bool_t reuseEventsWithinEnsemble, const Bool_t reuseEventsWithinExperiment = kFALSE, const Bool_t useReweighting = kFALSE); //! Embed full simulation events for the signal, rather than generating toy from the PDFs /*! \param [in] bkgndClass the name of the background class \param [in] fileName the name of the file containing the events \param [in] treeName the name of the tree \param [in] reuseEventsWithinEnsemble sample with replacement but only replace events once each experiment has been generated \param [in] reuseEventsWithinExperiment sample with immediate replacement \param [in] useReweighting perform an accept/reject routine using the configured signal amplitude model based on the MC-truth DP coordinate */ void embedBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName, const Bool_t reuseEventsWithinEnsemble, const Bool_t reuseEventsWithinExperiment = kFALSE, const Bool_t useReweighting = kFALSE); //! Set the value of the mixing phase /*! \param [in] phiMix the value of the mixing phase \param [in] fixPhiMix whether the value should be fixed or floated \param [in] useSinCos whether to use the sine and cosine as separate parameters or to just use the mixing phase itself */ void setPhiMix(const Double_t phiMix, const Bool_t fixPhiMix, const Bool_t useSinCos = kFALSE); //! Initialise the fit virtual void initialise(); //! Initialise the signal DP model virtual void initialiseDPModels(); //! Recalculate Normalization the signal DP models virtual void recalculateNormalisation(); //! Update the coefficients virtual void updateCoeffs(); // Toy MC generation and fitting overloaded functions virtual Bool_t genExpt(); //! Set the maximum value of A squared to be used in the accept/reject /*! \param [in] value the new value */ inline void setASqMaxValue(const Double_t value) {aSqMaxSet_ = value;} //! Weight events based on the DP model /*! \param [in] dataFileName the name of the data file \param [in] dataTreeName the name of the data tree */ virtual void weightEvents( const TString& dataFileName, const TString& dataTreeName ); //! Calculate things that depend on the fit parameters after they have been updated by Minuit virtual void propagateParUpdates(); //! Read in the input fit data variables, e.g. m13Sq and m23Sq virtual void cacheInputFitVars(); //! Check the initial fit parameters virtual void checkInitFitParams(); //! Get the fit results and store them /*! \param [in] tablePrefixName prefix for the name of the output file */ virtual void finaliseFitResults(const TString& tablePrefixName); //! Save the pdf Plots for all the resonances of experiment number fitExp /*! TODO - not working in this model!! \param [in] label prefix for the file name to be saved */ virtual void savePDFPlots(const TString& label); //! Save the pdf Plots for the sum of ressonances correspondint to "sin" of experiment number fitExp /*! TODO - not working in this model!! \param [in] label prefix for the file name to be saved \param [in] spin spin of the wave to be saved */ virtual void savePDFPlotsWave(const TString& label, const Int_t& spin); //! Print the fit fractions, total DP rate and mean efficiency /*! \param [out] output the stream to which to print */ virtual void printFitFractions(std::ostream& output); //! Print the asymmetries /*! \param [out] output the stream to which to print */ virtual void printAsymmetries(std::ostream& output); //! Write the fit results in latex table format /*! \param [in] outputFile the name of the output file */ virtual void writeOutTable(const TString& outputFile); //! Store the per event likelihood values virtual void storePerEvtLlhds(); // Methods to do with calculating the likelihood functions // and manipulating the fitting parameters. //! Get the total likelihood for each event /*! \param [in] iEvt the event number */ virtual Double_t getTotEvtLikelihood(const UInt_t iEvt); //! Calculate the signal and background likelihoods for the DP for a given event /*! \param [in] iEvt the event number */ virtual void getEvtDPDtLikelihood(const UInt_t iEvt); //! Determine the signal and background likelihood for the extra variables for a given event /*! \param [in] iEvt the event number */ virtual void getEvtExtraLikelihoods(const UInt_t iEvt); //! Get the total number of events /*! \return the total number of events */ virtual Double_t getEventSum() const; //! Set the fit parameters for the DP model void setSignalDPParameters(); //! Set the fit parameters for the decay time PDFs void setDecayTimeParameters(); //! Set the fit parameters for the extra PDFs void setExtraPdfParameters(); //! Set the initial yields void setFitNEvents(); //! Set the calibration parameters void setCalibParams(); //! Set the tagging efficiency parameters void setTagEffParams(); //! Set the asymmetry parameters void setAsymParams(); //! Set the tagging asymmetry parameters void setFlavTagAsymParams(); //! Set-up other parameters that are derived from the fit results, e.g. fit fractions void setExtraNtupleVars(); //! Set production and detections asymmetries /*! \param [in] AProd is the production asymmetry definied as (N_BqBar-N_Bq)/(N_BqBar+N_Bq) \param [in] AProdFix whether the production asymmetry should be fixed in the fit */ void setAsymmetries(const Double_t AProd, const Bool_t AProdFix); //! Set production and detections asymmetries /*! \param [in] bkgndClass the name of the background class \param [in] AProd is the production asymmetry definied as (N_BqBar-N_Bq)/(N_BqBar+N_Bq) \param [in] AProdFix whether the production asymmetry should be fixed in the fit */ void setBkgndAsymmetries(const TString& bkgndClass, const Double_t AProd, const Bool_t AProdFix); //! Randomise the initial fit parameters void randomiseInitFitPars(); //! Method to set up the storage for background-related quantities called by setBkgndClassNames virtual void setupBkgndVectors(); //! Calculate the CP asymmetries /*! \param [in] initValues is this before or after the fit */ void calcAsymmetries(const Bool_t initValues = kFALSE); //! Finalise value of mixing phase void checkMixingPhase(); protected: //! container to hold information on the number of events (and their weight) to generate for each category typedef std::map> LauGenInfo; //! container of background DP models typedef std::vector LauBkgndDPModelList; //! container of background PDFs typedef std::vector LauBkgndPdfsList; //! container of background yields or asymmetries typedef std::vector LauBkgndYieldList; //! container of boolean toggles for each background category typedef std::vector LauBkgndReuseEventsList; //! Determine the number of events to generate for each event category LauGenInfo eventsToGenerate(); //! Generate signal event /*! \return whether the generation was successful (kTRUE) or not (kFALSE) */ Bool_t generateSignalEvent(); //! Generate background event /*! \param [in] bgID ID number of the background class \return whether the generation was successful (kTRUE) or not (kFALSE) */ Bool_t generateBkgndEvent(UInt_t bgID); //! Setup the required ntuple branches void setupGenNtupleBranches(); //! Store all of the DP and decay time information void setDPDtBranchValues(); //! Generate from the extra PDFs /*! \param [in] extraPdfs the list of extra PDFs \param [in] embeddedData the embedded data sample */ void generateExtraPdfValues(LauPdfPList& extraPdfs, LauEmbeddedData* embeddedData); //! Add sPlot branches for the extra PDFs /*! \param [in] extraPdfs the list of extra PDFs \param [in] prefix the list of prefixes for the branch names */ void addSPlotNtupleBranches(const LauPdfPList& extraPdfs, const TString& prefix); //! Set the branches for the sPlot ntuple with extra PDFs /*! \param [in] extraPdfs the list of extra PDFs \param [in] prefix the prefix for the branch names \param [in] iEvt the event number */ Double_t setSPlotNtupleBranchValues(LauPdfPList& extraPdfs, const TString& prefix, const UInt_t iEvt); //! Update the signal events after Minuit sets background parameters void updateSigEvents(); //! Add branches to store experiment number and the event number within the experiment virtual void setupSPlotNtupleBranches(); //! Returns the names of all variables in the fit virtual LauSPlot::NameSet variableNames() const; //! Returns the names and yields of species that are free in the fit virtual LauSPlot::NumbMap freeSpeciesNames() const; //! Returns the names and yields of species that are fixed in the fit virtual LauSPlot::NumbMap fixdSpeciesNames() const; //! Returns the species and variables for all 2D PDFs in the fit virtual LauSPlot::TwoDMap twodimPDFs() const; //! Check if the signal is split into well-reconstructed and mis-reconstructed types virtual Bool_t splitSignal() const { return kFALSE; } //! Check if the mis-reconstructed signal is to be smeared in the DP virtual Bool_t scfDPSmear() const { return kFALSE; } //! Calculate the component integrals of the interference term void calcInterferenceTermIntegrals(); //! Calculate the total integral of the interference term void calcInterferenceTermNorm(); private: //! Keep access to the base class methods that we're further overloading using LauAbsFitModel::addFitParameters; //! Add the parameters from a decay time PDF into the fit /*! \param [in] decayTimePdf the container of PDFs */ UInt_t addFitParameters(LauDecayTimePdf* decayTimePdf); //! Add the parameters from each decay time PDF into the fit /*! \param [in] decayTimePdfList the container of PDFs */ UInt_t addFitParameters(std::vector& decayTimePdfList); //! Dalitz plot PDF for the antiparticle LauIsobarDynamics* sigModelB0bar_; //! Dalitz plot PDF for the particle LauIsobarDynamics* sigModelB0_; //! Kinematics object for antiparticle LauKinematics* kinematicsB0bar_; //! Kinematics object for particle LauKinematics* kinematicsB0_; //! The background Dalitz plot models for particles LauBkgndDPModelList BkgndDPModelsB_; //! The background Dalitz plot models for anti-particles LauBkgndDPModelList BkgndDPModelsBbar_; //! The background PDFs LauBkgndPdfsList BkgndPdfs_; //! Background boolean Bool_t usingBkgnd_; //! LauFlavTag object for flavour tagging LauFlavTag* flavTag_; //! Flavour tag for current event std::vector curEvtTagFlv_; //! Per event mistag for current event std::vector curEvtMistag_; //! Per event transformed mistag for current event std::vector curEvtMistagPrime_; //! True tag flavour (i.e. flavour at production) for the current event (only relevant for toy generation) LauFlavTag::Flavour curEvtTrueTagFlv_; //! Flavour at decay for the current event (only relevant for QFS) LauFlavTag::Flavour curEvtDecayFlv_; //! Number of signal components UInt_t nSigComp_; //! Number of signal DP parameters UInt_t nSigDPPar_; //! Number of decay time PDF parameters UInt_t nDecayTimePar_; //! Number of extra PDF parameters UInt_t nExtraPdfPar_; //! Number of normalisation parameters (yields, asymmetries) UInt_t nNormPar_; //! Number of calibration parameters (p0, p1) UInt_t nCalibPar_; //! Number of tagging efficneyc parameters UInt_t nTagEffPar_; //! Number of efficiency parameters (p0, p1) UInt_t nEffiPar_; //! Number of asymmetry parameters UInt_t nAsymPar_; //! Number of tagging asymmetry parameters UInt_t nTagAsymPar_; //! The complex coefficients for antiparticle std::vector coeffsB0bar_; //! The complex coefficients for particle std::vector coeffsB0_; //! Magnitudes and Phases or Real and Imaginary Parts std::vector coeffPars_; //! The integrals of the efficiency corrected amplitude cross terms for each pair of amplitude components /*! Calculated as the sum of A* x Abar x efficiency */ std::vector< std::vector > fifjEffSum_; //! The normalisation for the term 2.0*Re(A*Abar*phiMix) Double_t interTermReNorm_; //! The normalisation for the term 2.0*Im(A*Abar*phiMix) Double_t interTermImNorm_; //! The antiparticle fit fractions LauParArray fitFracB0bar_; //! The particle fit fractions LauParArray fitFracB0_; //! The fit fraction asymmetries std::vector fitFracAsymm_; //! A_CP parameter std::vector acp_; //! The mean efficiency for the antiparticle LauParameter meanEffB0bar_; //! The mean efficiency for the particle LauParameter meanEffB0_; //! The average DP rate for the antiparticle LauParameter DPRateB0bar_; //! The average DP rate for the particle LauParameter DPRateB0_; //! Signal yields LauParameter* signalEvents_; //! Signal asymmetry LauParameter* signalAsym_; //! Background yield(s) LauBkgndYieldList bkgndEvents_; //! Background asymmetries(s) LauBkgndYieldList bkgndAsym_; //! CP eigenvalue variable name TString cpevVarName_; //! CP eigenvalue for current event CPEigenvalue cpEigenValue_; //! Vector to store event CP eigenvalues std::vector evtCPEigenVals_; //! The mass difference between the neutral mass eigenstates LauParameter deltaM_; //! The width difference between the neutral mass eigenstates LauParameter deltaGamma_; //! The lifetime LauParameter tau_; //! The mixing phase LauParameter phiMix_; //! The sine of the mixing phase LauParameter sinPhiMix_; //! The cosine of the mixing phase LauParameter cosPhiMix_; //! Flag whether to use the sine and cosine of the mixing phase or the phase itself as the fit parameters Bool_t useSinCos_; //! e^{-i*phiMix} LauComplex phiMixComplex_; //! Signal Decay time PDFs (one per tagging category) LauDecayTimePdf* signalDecayTimePdf_; //! Background types std::vector BkgndTypes_; //! Background Decay time PDFs (one per tagging category) std::vector BkgndDecayTimePdfs_; //! Decay time for the current event Double_t curEvtDecayTime_; //! Decay time error for the current event Double_t curEvtDecayTimeErr_; //! PDFs for other variables LauPdfPList sigExtraPdf_; //! Production asymmetry between B0 and B0bar LauParameter AProd_; //! Production asymmetry between B0 and B0bar for bkgnds std::vector AProdBkgnd_; // Toy generation stuff //! The maximum allowed number of attempts when generating an event Int_t iterationsMax_; //! The number of unsucessful attempts to generate an event so far Int_t nGenLoop_; //! The value of A squared for the current event Double_t ASq_; //! The maximum value of A squared that has been seen so far while generating Double_t aSqMaxVar_; //! The maximum allowed value of A squared Double_t aSqMaxSet_; //! Flag for storage of amplitude info in generated ntuple Bool_t storeGenAmpInfo_; //! The signal event tree for embedding fully simulated events LauEmbeddedData* signalTree_; //! The background event tree for embedding fully simulated events std::vector bkgndTree_; //! Boolean to control reuse of embedded signal events Bool_t reuseSignal_; //! Boolean to use reweighting when embedding Bool_t useReweighting_; //! Vector of booleans to reuse background events LauBkgndReuseEventsList reuseBkgnd_; // Likelihood values //! Signal DP likelihood value Double_t sigDPLike_; //! Signal likelihood from extra PDFs Double_t sigExtraLike_; //! Total signal likelihood Double_t sigTotalLike_; //! Background DP likelihood value(s) std::vector bkgndDPLike_; //! Background likelihood value(s) from extra PDFs std::vector bkgndExtraLike_; //! Total background likelihood(s) std::vector bkgndTotalLike_; ClassDef(LauTimeDepFitModel,0) // Time-dependent neutral model }; #endif diff --git a/src/LauDecayTimePdf.cc b/src/LauDecayTimePdf.cc index f7dc0dd..bbd9896 100644 --- a/src/LauDecayTimePdf.cc +++ b/src/LauDecayTimePdf.cc @@ -1,614 +1,614 @@ /* Copyright 2006 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauDecayTimePdf.cc \brief File containing implementation of LauDecayTimePdf class. */ #include #include #include #include #include #include #include "TMath.h" #include "TRandom.h" #include "TSystem.h" #include "TH1.h" #include "RooMath.h" #include "Lau1DCubicSpline.hh" #include "Lau1DHistPdf.hh" #include "LauConstants.hh" #include "LauComplex.hh" #include "LauDecayTimePdf.hh" #include "LauFitDataTree.hh" #include "LauParameter.hh" #include "LauParamFixed.hh" #include "LauRandom.hh" #include "LauNonSmearedDecayTimeCalculator.hh" #include "LauNonSmearedUniformEfficiencyDecayTimeIntegrator.hh" #include "LauNonSmearedBinnedEfficiencyDecayTimeIntegrator.hh" #include "LauNonSmearedSplineEfficiencyDecayTimeIntegrator.hh" #include "LauSmearedDecayTimeCalculator.hh" #include "LauSmearedUniformEfficiencyDecayTimeIntegrator.hh" #include "LauSmearedBinnedEfficiencyDecayTimeIntegrator.hh" #include "LauSmearedSplineEfficiencyDecayTimeIntegrator.hh" ClassImp(LauDecayTimePdf) LauDecayTimePdf::LauDecayTimePdf(const TString& theVarName, const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const TH1* dtHist, const LauDecayTime::TimeMeasurementMethod method) : type_{LauDecayTime::FuncType::Hist}, method_{method}, varName_{theVarName}, minAbscissa_{minAbscissaVal}, maxAbscissa_{maxAbscissaVal}, pdfHist_{ (dtHist==nullptr) ? nullptr : std::make_unique( varName_, dtHist, minAbscissa_, maxAbscissa_ ) } // TODO - check this, once we've consolidated all the data members { } LauDecayTimePdf::LauDecayTimePdf(const TString& theVarName, const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const TString& theVarErrName, const Double_t minAbscissaErr, const Double_t maxAbscissaErr, const TH1* dtHist, const TH1* dtErrHist, const LauDecayTime::TimeMeasurementMethod method) : type_{LauDecayTime::FuncType::Hist}, method_{method}, varName_{theVarName}, varErrName_{theVarErrName}, minAbscissa_{minAbscissaVal}, maxAbscissa_{maxAbscissaVal}, minAbscissaError_{minAbscissaErr}, maxAbscissaError_{maxAbscissaErr}, pdfHist_{ (dtHist==nullptr) ? nullptr : std::make_unique( varName_, dtHist, minAbscissa_, maxAbscissa_ ) }, errHist_{ (dtErrHist==nullptr) ? nullptr : std::make_unique( varErrName_, dtErrHist, minAbscissaError_, maxAbscissaError_ ) } // TODO - check this, once we've consolidated all the data members { } LauDecayTimePdf::LauDecayTimePdf(const TString& theVarName, const Double_t minAbscissaVal, const Double_t maxAbscissaVal, std::unique_ptr physicsModel, const LauDecayTime::TimeMeasurementMethod method) : type_{physicsModel->getFunctionType()}, method_{method}, varName_{theVarName}, minAbscissa_{minAbscissaVal}, maxAbscissa_{maxAbscissaVal}, physicsModel_{std::move(physicsModel)} // TODO - check this, once we've consolidated all the data members { auto effModel = std::make_unique(1.0); calculator_ = std::make_unique( minAbscissa_, maxAbscissa_, *physicsModel_ ); integrator_ = std::make_unique( minAbscissa_, maxAbscissa_, *physicsModel_, *effModel ); efficiencyModel_ = std::move( effModel ); } LauDecayTimePdf::LauDecayTimePdf(const TString& theVarName, const Double_t minAbscissaVal, const Double_t maxAbscissaVal, std::unique_ptr physicsModel, std::unique_ptr resolutionModel, const LauDecayTime::TimeMeasurementMethod method) : type_{physicsModel->getFunctionType()}, method_{method}, varName_{theVarName}, minAbscissa_{minAbscissaVal}, maxAbscissa_{maxAbscissaVal}, physicsModel_{std::move(physicsModel)}, resolutionModel_{std::move(resolutionModel)}, smear_{true} // TODO - check this, once we've consolidated all the data members { auto effModel = std::make_unique(1.0); calculator_ = std::make_unique( minAbscissa_, maxAbscissa_, *physicsModel_, *resolutionModel_ ); integrator_ = std::make_unique( minAbscissa_, maxAbscissa_, *physicsModel_, *effModel, *resolutionModel_ ); efficiencyModel_ = std::move( effModel ); } LauDecayTimePdf::LauDecayTimePdf(const TString& theVarName, const Double_t minAbscissaVal, const Double_t maxAbscissaVal, const TString& theVarErrName, const Double_t minAbscissaErr, const Double_t maxAbscissaErr, std::unique_ptr physicsModel, std::unique_ptr resolutionModel, const TH1* dtErrHist, const LauDecayTime::TimeMeasurementMethod method) : type_{physicsModel->getFunctionType()}, method_{method}, varName_{theVarName}, varErrName_{theVarErrName}, minAbscissa_{minAbscissaVal}, maxAbscissa_{maxAbscissaVal}, minAbscissaError_{minAbscissaErr}, maxAbscissaError_{maxAbscissaErr}, physicsModel_{std::move(physicsModel)}, resolutionModel_{std::move(resolutionModel)}, smear_{true}, scaleWithPerEventError_{true}, errHist_{ (dtErrHist==nullptr) ? nullptr : std::make_unique( varErrName_, dtErrHist, minAbscissaError_, maxAbscissaError_ ) } // TODO - check this, once we've consolidated all the data members { auto effModel = std::make_unique(1.0); calculator_ = std::make_unique( minAbscissa_, maxAbscissa_, *physicsModel_, *resolutionModel_ ); integrator_ = std::make_unique( minAbscissa_, maxAbscissa_, *physicsModel_, *effModel, *resolutionModel_ ); efficiencyModel_ = std::move( effModel ); } void LauDecayTimePdf::initialise() { if (type_ == LauDecayTime::FuncType::Hist) { if ( not pdfHist_ ) { std::cerr << "ERROR in LauDecayTimePdf::initialise : Hist type should have a histogram PDF supplied" << std::endl; gSystem->Exit(EXIT_FAILURE); } return; } if ( smear_ and scaleWithPerEventError_ ) { if ( not errHist_ ) { std::cerr << "ERROR in LauDecayTimePdf::initialise : scaling with per-event error but no error distribution supplied" << std::endl; gSystem->Exit(EXIT_FAILURE); } } // TODO other consistency checks? // Initialise the physics model physicsModel_->initialise(); // Initialise the resolution model if ( resolutionModel_ ) { resolutionModel_->initialise(); } // Initialise the efficiency model efficiencyModel_->initialise(); // Get all the parameters and consolidate them so we can pass them to the fit model std::vector physicsPars { physicsModel_->getParameters() }; physicsParFloating_ = this->anyParFloating( physicsPars ); this->addParams( physicsPars ); if ( resolutionModel_ ) { std::vector resoPars { resolutionModel_->getParameters() }; resoParFloating_ = this->anyParFloating( resoPars ); this->addParams( resoPars ); } std::vector effiPars { efficiencyModel_->getParameters() }; effiParFloating_ = this->anyParFloating( effiPars ); this->addParams( effiPars ); anythingFloating_ = physicsParFloating_ or resoParFloating_ or effiParFloating_; } void LauDecayTimePdf::cacheInfo(const LauFitDataTree& inputData) { const std::size_t nEvents { inputData.nEvents() }; // If we're a histogram form then there's not so much to do if ( type_ == LauDecayTime::FuncType::Hist ) { // Pass the data to the decay-time PDF for caching pdfHist_->cacheInfo(inputData); // Pass the data to the decay-time error PDF for caching if ( errHist_ ) { errHist_->cacheInfo(inputData); } // Make cache of effiTerms // Efficiency will always be 1 by definition effiTerms_.assign(nEvents,1.0); return; } // Otherwise... // Check that the input data contains the decay time variable bool hasBranch { inputData.haveBranch( varName_ ) }; if (!hasBranch) { std::cerr << "ERROR in LauDecayTimePdf::cacheInfo : Input data does not contain variable \"" << varName_ << "\"" << std::endl; gSystem->Exit(EXIT_FAILURE); } // If we're scaling by the per-event error, also check that the input data contains the decay time error variable const bool needPerEventNormTerms { smear_ and scaleWithPerEventError_ }; if ( needPerEventNormTerms ) { hasBranch = inputData.haveBranch( varErrName_ ); if (!hasBranch) { std::cerr << "ERROR in LauDecayTimePdf::cacheInfo : Input data does not contain variable \"" << varErrName_ << "\"" << std::endl; gSystem->Exit(EXIT_FAILURE); } // Pass the data to the decay-time error PDF for caching errHist_->cacheInfo(inputData); } // Clear the vectors and reserve enough space in the caches of the terms abscissas_.clear(); abscissas_.resize(nEvents); abscissaErrors_.clear(); abscissaErrors_.resize(nEvents); termsStore_.clear(); termsStore_.resize(nEvents); effiTerms_.clear(); effiTerms_.resize(nEvents); // Correctly size the normalisation cache elements // depending on whether we're doing per-event resolution normTermsStore_.clear(); if ( needPerEventNormTerms ) { normTermsStore_.resize(nEvents); } else { normTermsStore_.resize(1); } // Determine the abscissa and abscissa error values for each event for ( std::size_t iEvt {0}; iEvt < nEvents; iEvt++ ) { const LauFitData& dataValues { inputData.getData(iEvt) }; const Double_t abscissa { dataValues.at( varName_ ) }; if ( abscissa > maxAbscissa_ or abscissa < minAbscissa_ ) { std::cerr << "ERROR in LauDecayTimePdf::cacheInfo : Given value of the decay time: " << abscissa << " is outside allowed range: [" << minAbscissa_ << "," << maxAbscissa_ << "]" << std::endl; gSystem->Exit(EXIT_FAILURE); } abscissas_[iEvt] = abscissa; const Double_t abscissaErr { needPerEventNormTerms ? dataValues.at( varErrName_ ) : 0.0 }; if ( errHist_ and ( abscissaErr > maxAbscissaError_ or abscissaErr < minAbscissaError_ ) ) { std::cerr << "ERROR in LauDecayTimePdf::cacheInfo : Given value of the decay-time error: " << abscissaErr << " is outside allowed range: [" << minAbscissaError_ << "," << maxAbscissaError_ << "]." << std::endl; gSystem->Exit(EXIT_FAILURE); } abscissaErrors_[iEvt] = abscissaErr; } // Cache the numerator terms calculator_->cacheInfo( abscissas_, abscissaErrors_ ); // Cache the normalisations integrator_->cacheInfo( abscissaErrors_ ); // Force storing of all info this first time around physicsParChanged_ = true; resoParChanged_ = true; effiParChanged_ = true; this->updateCache(); } bool LauDecayTimePdf::anyParFloating( const std::vector& pars ) const { LauParamFixed isFixed; return not std::all_of( pars.begin(), pars.end(), isFixed ); } void LauDecayTimePdf::addParams( std::vector& pars ) { for ( auto& par : pars ) { params_.push_back( par ); } } void LauDecayTimePdf::updateCache() { // Calculate the values of all terms for each event const std::size_t nEvents { abscissas_.size() }; // If any of the physics or resolution parameters have changed we need // to update all numerator terms if ( physicsParChanged_ or resoParChanged_ ) { for (std::size_t iEvt {0}; iEvt < nEvents; iEvt++) { termsStore_[iEvt] = calculator_->getTerms( iEvt ); } } // If any of the efficiency parameters have changed we need to // recalculate the efficiency if ( effiParChanged_ ) { for (std::size_t iEvt {0}; iEvt < nEvents; iEvt++) { effiTerms_[iEvt] = this->calcEffiTerm( abscissas_[iEvt] ); } } // Calculate the normalisation terms // If we're not doing per-event scaling, // we only need to calculate the normalisations once const std::array nNormEntries { 1, nEvents }; const bool needPerEventNormTerms { smear_ and scaleWithPerEventError_ }; for ( std::size_t iEvt{0}; iEvt < nNormEntries[needPerEventNormTerms]; ++iEvt ) { normTermsStore_[iEvt] = integrator_->getNormTerms( iEvt ); } } void LauDecayTimePdf::calcLikelihoodInfo(const std::size_t iEvt) { // Extract all the terms and their normalisations if (type_ == LauDecayTime::FuncType::Hist) { pdfHist_->calcLikelihoodInfo(iEvt); histTerm_ = pdfHist_->getLikelihood(); } else { terms_ = termsStore_[iEvt]; const bool needPerEventNormTerms { smear_ and scaleWithPerEventError_ }; const std::size_t normTermElement { needPerEventNormTerms * iEvt }; normTerms_ = normTermsStore_[normTermElement]; } // Extract the decay time error PDF value if ( errHist_ ) { errHist_->calcLikelihoodInfo(iEvt); errTerm_ = errHist_->getLikelihood(); } else { errTerm_ = 1.0; } // Extract the decay time efficiency effiTerm_ = effiTerms_[iEvt]; } void LauDecayTimePdf::calcLikelihoodInfo( const Double_t abscissa ) { // Check whether any of the gaussians should be scaled - if any of them should we need the per-event error if ( smear_ and scaleWithPerEventError_ ) { std::cerr << "ERROR in LauDecayTimePdf::calcLikelihoodInfo : Per-event error on decay time not provided, cannot calculate anything." << std::endl; return; } this->calcLikelihoodInfo(abscissa, 0.0); } Double_t LauDecayTimePdf::calcEffiTerm( const Double_t abscissa ) const { Double_t effiTerm { efficiencyModel_->getEfficiency( abscissa ) }; // TODO print warning messages for these, but not every time // - do we even want the upper limit imposed? if ( effiTerm > 1.0 ) { effiTerm = 1.0; } // else if ( effiTerm < 0.0 ) { effiTerm = 0.0; } else if ( effiTerm <= 0.0 ) { effiTerm = 1e-20; } return effiTerm; } void LauDecayTimePdf::calcLikelihoodInfo(const Double_t abscissa, const Double_t abscissaErr) { // Check that the decay time and the decay time error are in valid ranges if ( abscissa > maxAbscissa_ or abscissa < minAbscissa_ ) { std::cerr << "ERROR in LauDecayTimePdf::calcLikelihoodInfo : Given value of the decay time: " << abscissa << " is outside allowed range: [" << minAbscissa_ << "," << maxAbscissa_ << "]" << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( errHist_ and ( abscissaErr > maxAbscissaError_ or abscissaErr < minAbscissaError_ ) ) { std::cerr << "ERROR in LauDecayTimePdf::calcLikelihoodInfo : Given value of the decay-time error: " << abscissaErr << " is outside allowed range: [" << minAbscissaError_ << "," << maxAbscissaError_ << "]." << std::endl; gSystem->Exit(EXIT_FAILURE); } // Calculate the decay time error PDF value if ( errHist_ ) { const std::vector absErrVec {abscissaErr}; errHist_->calcLikelihoodInfo(absErrVec); errTerm_ = errHist_->getLikelihood(); } else { errTerm_ = 1.0; } // For the histogram PDF just calculate that term and return if ( type_ == LauDecayTime::FuncType::Hist ){ pdfHist_->calcLikelihoodInfo(abscissa); histTerm_ = pdfHist_->getLikelihood(); effiTerm_ = 1.0; return; } // Determine the decay time efficiency effiTerm_ = this->calcEffiTerm( abscissa ); // Determine the other terms terms_ = calculator_->calcTerms( abscissa, {1.0, abscissaErr} ); // TODO - do we need to do this? //if ( smear_ and scaleWithPerEventError_ ) { //normTerms_ = integrator_->calcNormTerms( {1.0, abscissaErr} ); //} } Double_t LauDecayTimePdf::generateError( bool forceNew ) { if ( errHist_ && ( forceNew or not abscissaErrorGenerated_ ) ) { LauFitData errData { errHist_->generate(nullptr) }; abscissaError_ = errData.at( varErrName_ ); abscissaErrorGenerated_ = true; } return abscissaError_; } Double_t LauDecayTimePdf::generate( const LauKinematics* kinematics ) { // generateError SHOULD have been called before this // function but will call it here just to make sure // (has no effect if has already been called) abscissaError_ = this->generateError(); Double_t abscissa{0.0}; switch ( type_ ) { case LauDecayTime::FuncType::ExpTrig : case LauDecayTime::FuncType::ExpHypTrig : std::cerr << "ERROR in LauDecayTimePdf::generate : direct generation does not make sense for ExpTrig and ExpHypTrig types" << std::endl; return abscissa; case LauDecayTime::FuncType::Hist : { const LauFitData genAbscissa { pdfHist_ -> generate( kinematics ) }; abscissa = genAbscissa.at(this->varName()); break; } case LauDecayTime::FuncType::Exp : { do { abscissa = LauRandom::randomFun()->Uniform(minAbscissa_,maxAbscissa_); this->calcLikelihoodInfo( abscissa, abscissaError_ ); const Double_t pdfVal { this->getExpTerm() * this->getEffiTerm() }; const Double_t maxVal { this->getMaxHeight(kinematics) }; if ( pdfVal > maxVal ) { std::cerr << "WARNING in LauDecayTimePdf::generate : PDF value = " << pdfVal << " is larger than the maximum PDF height " << maxVal << "\n"; std::cerr << " : This occurs for the abscissa = " << abscissa << "\n"; std::cerr << " : This really should not happen!!" << std::endl; } if ( LauRandom::randomFun()->Rndm() <= pdfVal/maxVal ) { break; } } while ( true ); break; } case LauDecayTime::FuncType::Delta : { // TODO std::cerr << "WARNING in LauDecayTimePdf::generate : generation of Delta case not currently implemented" << std::endl; break; } case LauDecayTime::FuncType::DeltaExp : { // TODO std::cerr << "WARNING in LauDecayTimePdf::generate : generation of DeltaExp case not currently implemented" << std::endl; break; } } // mark that we need a new error to be generated next time abscissaErrorGenerated_ = false; return abscissa; } Double_t LauDecayTimePdf::getMaxHeight(const LauKinematics* /*kinematics*/) { if ( heightUpToDate_ ) { return maxHeight_; } // TODO this is brute force - can we do this in a more refined way? // Scan in small increments across the space to find the maximum height const std::size_t nPoints { 1000 }; const Double_t range { maxAbscissa_ - minAbscissa_ }; const Double_t delta { range / nPoints }; maxHeight_ = 0.0; for ( Double_t point {minAbscissa_}; point <= maxAbscissa_; point += delta ) { this->calcLikelihoodInfo(point, abscissaError_); Double_t heightAtPoint { 0.0 }; if ( type_ == LauDecayTime::FuncType::Exp ) { heightAtPoint = this->getExpTerm(); } else { // TODO - implement the Delta and ExpDelta cases std::cerr << "WARNING in LauDecayTimePdf::getMaxHeight : only Exp case currently implemented" << std::endl; } heightAtPoint *= this->getEffiTerm(); if ( heightAtPoint > maxHeight_ ) { maxHeight_ = heightAtPoint; } } // Mutliply by 120% to be on the safe side maxHeight_ *= 1.2; // Mark the height as being up to date // (unless we're scaling by per-event error) heightUpToDate_ = not ( smear_ and scaleWithPerEventError_ ); return maxHeight_; } void LauDecayTimePdf::setBinnedEfficiency( std::unique_ptr effModel ) { if ( not effModel ) { std::cerr << "WARNING in LauDecayTimePdf::setBinnedEfficiency : supplied binned efficiency model pointer is null" << std::endl; return; } // Create the approptiate integrator // NB need to use effModel here since it needs to be of concrete type if ( smear_ ) { integrator_ = std::make_unique( minAbscissa_, maxAbscissa_, *physicsModel_, *effModel, *resolutionModel_ ); } else { integrator_ = std::make_unique( minAbscissa_, maxAbscissa_, *physicsModel_, *effModel ); } // Store the efficiency model (as a pointer to base) efficiencyModel_ = std::move(effModel); } void LauDecayTimePdf::updatePulls() { - for ( std::vector::iterator iter = params_.begin(); iter != params_.end(); ++iter ) { - std::vector params = (*iter)->getPars(); - for (std::vector::iterator params_iter = params.begin(); params_iter != params.end(); ++params_iter ) { - if (!(*iter)->fixed()) { - (*params_iter)->updatePull(); + for ( LauAbsRValue* rValPtr : params_ ) { + std::vector params = rValPtr->getPars(); + for ( LauParameter* param : params ) { + if ( not param->fixed() ) { + param->updatePull(); } } } } void LauDecayTimePdf::propagateParUpdates() { // If none of the parameters are floating there's nothing to do if ( not anythingFloating_ ) { return; } physicsModel_->propagateParUpdates(); physicsParChanged_ = physicsModel_->anythingChanged(); if ( resolutionModel_ ) { resolutionModel_->propagateParUpdates(); resoParChanged_ = resolutionModel_->anythingChanged(); } efficiencyModel_->propagateParUpdates(); effiParChanged_ = efficiencyModel_->anythingChanged(); // If nothing has changed, there's nothing to do anythingChanged_ = physicsParChanged_ or resoParChanged_ or effiParChanged_; if ( not anythingChanged_ ) { return; } calculator_->propagateParUpdates(); integrator_->propagateParUpdates(); // Otherwise we need to update the cache this->updateCache(); } diff --git a/src/LauFlavTag.cc b/src/LauFlavTag.cc index 0dab8de..b49649c 100644 --- a/src/LauFlavTag.cc +++ b/src/LauFlavTag.cc @@ -1,1428 +1,1465 @@ /* Copyright 2017 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauFlavTag.cc \brief File containing implementation of LauFlavTag class. */ #include #include #include #include "TMath.h" #include "TString.h" #include "TSystem.h" #include "Lau1DHistPdf.hh" #include "LauAbsPdf.hh" #include "LauFlavTag.hh" #include "LauRandom.hh" ClassImp(LauFlavTag) std::ostream& operator<<( std::ostream& stream, const LauFlavTag::BkgndType bkgndType ) { switch ( bkgndType ) { case LauFlavTag::BkgndType::Combinatorial : stream << "Combinatorial"; break; case LauFlavTag::BkgndType::FlavourSpecific : stream << "FlavourSpecific"; break; case LauFlavTag::BkgndType::SelfConjugate : stream << "SelfConjugate"; break; case LauFlavTag::BkgndType::NonSelfConjugate : stream << "NonSelfConjugate"; break; } return stream; } LauFlavTag::LauFlavTag(const Bool_t useAveDelta, const Bool_t useEtaPrime, const std::map bkgndInfo) : useAveDelta_(useAveDelta), useEtaPrime_(useEtaPrime) { // Put map values into vectors bkgndTypes_.reserve( bkgndInfo.size() ); bkgndDecayFlvDep_.reserve( bkgndInfo.size() ); for (const auto& [ bkgndName, bkgndType ] : bkgndInfo){ const std::size_t bkgndPos { bkgndIndex_.size() }; bkgndIndex_.insert(std::make_pair(bkgndName, bkgndPos)); bkgndTypes_.push_back(bkgndType); bkgndDecayFlvDep_.push_back(false); std::cout<< "INFO in LauFlavTag::LauFlavTag : adding background " << bkgndName << " of type " << bkgndType < LauFlavTag::getBkgndNames() { std::vector names; names.reserve( bkgndIndex_.size() ); for ( auto& [ name, _ ] : bkgndIndex_ ) { names.push_back( name ); } return names; } void LauFlavTag::addTagger(const TString& name, const TString& tagVarName, const TString& mistagVarName, LauAbsPdf* etapdf, const std::pair tagEff, const std::pair calib_p0, const std::pair calib_p1) { // Check that we don't already have a tagger with the same name if ( taggerIndex_.find(name) != taggerIndex_.end() ) { std::cerr << "ERROR in LauFlavTag::addTagger : tagger called " << name << " already added" << std::endl; gSystem->Exit(EXIT_FAILURE); } // Check that the PDF pointer is valid if ( not etapdf ) { std::cerr << "ERROR in LauFlavTag::addTagger : Eta PDF pointer is NULL" << std::endl; gSystem->Exit(EXIT_FAILURE); } // Find how many taggers have already been added const std::size_t index { tagVarNames_.size() }; // Update map to relate tagger name and index in the vectors taggerIndex_[name] = index; // Extend vectors this->extendVectors(tagVarName, mistagVarName); etaPdfs_.push_back(etapdf); Lau1DHistPdf* etahistpdf = dynamic_cast(etapdf); if (etahistpdf){ perEvtAvgMistag_.push_back(etahistpdf->getMean()); } else { std::cerr << "WARNING in LauFlavTag::addTagger : Couldn't determine average eta value from PDF. Setting it to 0.4." << std::endl; perEvtAvgMistag_.push_back(0.4); } //Calib parameters this->setupCalibParams(name,calib_p0,calib_p1); //Tagging efficiencies if (!useAveDelta_){ TString tagEff_b0Name("tagEff_b0_"+name); TString tagEff_b0barName("tagEff_b0bar_"+name); LauParameter* tageffb0 = new LauParameter(tagEff_b0Name,tagEff.first,0.0,1.0,kTRUE); tagEff_B0_.push_back(tageffb0); tagEff_B0_[index]->initValue(tagEff.first); tagEff_B0_[index]->genValue(tagEff.first); tagEff_B0_[index]->fixed(kTRUE); if (tagEff.second==-1.0){ tagEff_B0bar_.push_back(tagEff_B0_[index]->createClone(tagEff_b0barName)); } else { LauParameter* tageffb0bar = new LauParameter(tagEff_b0barName,tagEff.second,0.0,1.0,kTRUE); tagEff_B0bar_.push_back(tageffb0bar); tagEff_B0bar_[index]->initValue(tagEff.second); tagEff_B0bar_[index]->genValue(tagEff.second); tagEff_B0bar_[index]->fixed(kTRUE); } } else { //Use average and delta variables TString tagEff_aveName("tagEff_ave_"+name); TString tagEff_deltaName("tagEff_delta_"+name); LauParameter* tageffave = new LauParameter(tagEff_aveName,tagEff.first,0.0,1.0,kTRUE); tagEff_ave_.push_back(tageffave); tagEff_ave_[index]->initValue(tagEff.first); tagEff_ave_[index]->genValue(tagEff.first); tagEff_ave_[index]->fixed(kTRUE); LauParameter* tageffdelta = new LauParameter(tagEff_deltaName,tagEff.second,-1.0,1.0,kTRUE); tagEff_delta_.push_back(tageffdelta); tagEff_delta_[index]->initValue(tagEff.second); tagEff_delta_[index]->genValue(tagEff.second); tagEff_delta_[index]->fixed(kTRUE); } std::cout<<"INFO in LauFlavTag::addTagger : Added tagger with name "<< name << std::endl; } void LauFlavTag::addTagger(const TString& name, const TString& tagVarName, const TString& mistagVarName, LauAbsPdf* etapdf, const std::pair tagEff, const std::pair calib_p0, const std::pair calib_p1) { // Check that we don't already have a tagger with the same name if ( taggerIndex_.find(name) != taggerIndex_.end() ) { std::cerr << "ERROR in LauFlavTag::addTagger : tagger called " << name << " already added" << std::endl; gSystem->Exit(EXIT_FAILURE); } // Check that the PDF pointer is valid if ( not etapdf ) { std::cerr << "ERROR in LauFlavTag::addTagger : Eta PDF pointer is NULL" << std::endl; gSystem->Exit(EXIT_FAILURE); } // Find how many taggers have already been added const std::size_t index { tagVarNames_.size() }; // Update map to relate tagger name and index in the vectors taggerIndex_[name] = index; // Extend vectors this->extendVectors(tagVarName, mistagVarName); etaPdfs_.push_back(etapdf); Lau1DHistPdf* etahistpdf = dynamic_cast(etapdf); if (etahistpdf){ perEvtAvgMistag_.push_back(etahistpdf->getMean()); } else { std::cerr << "WARNING in LauFlavTag::addTagger : Couldn't determine average eta value from PDF. Setting it to 0.4." << std::endl; perEvtAvgMistag_.push_back(0.4); } //Calib parameters this->setupCalibParams(name,calib_p0,calib_p1); //Tagging efficiencies if (!useAveDelta_){ tagEff_hist_B0_.push_back(tagEff.first); tagEff_B0_.push_back(nullptr); if (tagEff.second==nullptr){ tagEff_hist_B0bar_.push_back(tagEff.first); tagEff_B0bar_.push_back(nullptr); } else { tagEff_hist_B0bar_.push_back(tagEff.second); tagEff_B0bar_.push_back(nullptr); } } else { //Use average and delta variables tagEff_hist_ave_.push_back(tagEff.first); tagEff_hist_delta_.push_back(tagEff.second); tagEff_ave_.push_back(nullptr); tagEff_delta_.push_back(nullptr); } std::cout<<"INFO in LauFlavTag::addTagger : Added tagger with name "<< name << std::endl; } void LauFlavTag::extendVectors(const TString& tagVarName, const TString& mistagVarName) { tagVarNames_.push_back(tagVarName); mistagVarNames_.push_back(mistagVarName); curEvtTagFlv_.push_back(Flavour::Unknown); curEvtMistag_.push_back(Flavour::Unknown); if ( not bkgndIndex_.empty() ){ if (!useAveDelta_){ //Loop over the outer vector and extend the inner vectors for (auto& innerVec : tagEffBkgnd_B0_){ innerVec.push_back({nullptr,nullptr,nullptr}); } for (auto& innerVec : tagEffBkgnd_hist_B0_){ innerVec.push_back({nullptr,nullptr,nullptr}); } for (auto& innerVec : tagEffBkgnd_B0bar_){ innerVec.push_back({nullptr,nullptr,nullptr}); } for (auto& innerVec : tagEffBkgnd_hist_B0bar_){ innerVec.push_back({nullptr,nullptr,nullptr}); } } else { //Loop over the outer vector and extend the inner vectors for (auto& innerVec : tagEffBkgnd_ave_){ innerVec.push_back({nullptr,nullptr,nullptr}); } for (auto& innerVec : tagEffBkgnd_hist_ave_){ innerVec.push_back({nullptr,nullptr,nullptr}); } for (auto& innerVec : tagEffBkgnd_delta_){ innerVec.push_back({nullptr,nullptr,nullptr}); } for (auto& innerVec : tagEffBkgnd_hist_delta_){ innerVec.push_back({nullptr,nullptr,nullptr}); } } for (auto& innerVec : etaBkgndPdfs_){ innerVec.push_back({nullptr,nullptr,nullptr}); } for (auto& innerVec : avgMistagBkgnd_){ innerVec.push_back({0.0,0.0,0.0}); } } } void LauFlavTag::setupCalibParams(const TString& taggerName, const std::pair calib_p0, const std::pair calib_p1) { const std::size_t taggerID { taggerIndex_.at( taggerName ) }; if (!useAveDelta_){ TString calib_p0_b0Name("calib_p0_b0_"+taggerName); TString calib_p0_b0barName("calib_p0_b0bar_"+taggerName); TString calib_p1_b0Name("calib_p1_b0_"+taggerName); TString calib_p1_b0barName("calib_p1_b0bar_"+taggerName); LauParameter* calibp0b0 = new LauParameter(calib_p0_b0Name,calib_p0.first,-10.0,10.0,kTRUE); calib_p0_B0_.push_back(calibp0b0); calib_p0_B0_[taggerID]->initValue(calib_p0.first); calib_p0_B0_[taggerID]->genValue(calib_p0.first); calib_p0_B0_[taggerID]->fixed(kTRUE); LauParameter* calibp1b0 = new LauParameter(calib_p1_b0Name,calib_p1.first,0.0,1.5,kTRUE); calib_p1_B0_.push_back(calibp1b0); calib_p1_B0_[taggerID]->initValue(calib_p1.first); calib_p1_B0_[taggerID]->genValue(calib_p1.first); calib_p1_B0_[taggerID]->fixed(kTRUE); if (calib_p0.second==-1.0 && calib_p1.second==-1.0){ calib_p0_B0bar_.push_back(calib_p0_B0_[taggerID]->createClone(calib_p0_b0barName)); calib_p1_B0bar_.push_back(calib_p1_B0_[taggerID]->createClone(calib_p1_b0barName)); } else { LauParameter* calibp0b0bar = new LauParameter(calib_p0_b0barName,calib_p0.second,-10.0,10.0,kTRUE); calib_p0_B0bar_.push_back(calibp0b0bar); calib_p0_B0bar_[taggerID]->initValue(calib_p0.second); calib_p0_B0bar_[taggerID]->genValue(calib_p0.second); calib_p0_B0bar_[taggerID]->fixed(kTRUE); LauParameter* calibp1b0bar = new LauParameter(calib_p1_b0barName,calib_p1.second,0.0,1.5,kTRUE); calib_p1_B0bar_.push_back(calibp1b0bar); calib_p1_B0bar_[taggerID]->initValue(calib_p1.second); calib_p1_B0bar_[taggerID]->genValue(calib_p1.second); calib_p1_B0bar_[taggerID]->fixed(kTRUE); } } else { //Use average and delta variables TString calib_p0_aveName("calib_p0_ave_"+taggerName); TString calib_p0_deltaName("calib_p0_delta_"+taggerName); TString calib_p1_aveName("calib_p1_ave_"+taggerName); TString calib_p1_deltaName("calib_p1_delta_"+taggerName); LauParameter* calibp0ave = new LauParameter(calib_p0_aveName,calib_p0.first,-10.0,10.0,kTRUE); calib_p0_ave_.push_back(calibp0ave); calib_p0_ave_[taggerID]->initValue(calib_p0.first); calib_p0_ave_[taggerID]->genValue(calib_p0.first); calib_p0_ave_[taggerID]->fixed(kTRUE); LauParameter* calibp1ave = new LauParameter(calib_p1_aveName,calib_p1.first,0.0,1.5,kTRUE); calib_p1_ave_.push_back(calibp1ave); calib_p1_ave_[taggerID]->initValue(calib_p1.first); calib_p1_ave_[taggerID]->genValue(calib_p1.first); calib_p1_ave_[taggerID]->fixed(kTRUE); LauParameter* calibp0delta = new LauParameter(calib_p0_deltaName,calib_p0.second,-10.0,10.0,kTRUE); calib_p0_delta_.push_back(calibp0delta); calib_p0_delta_[taggerID]->initValue(calib_p0.second); calib_p0_delta_[taggerID]->genValue(calib_p0.second); calib_p0_delta_[taggerID]->fixed(kTRUE); LauParameter* calibp1delta = new LauParameter(calib_p1_deltaName,calib_p1.second,-10.0,10.0,kTRUE); calib_p1_delta_.push_back(calibp1delta); calib_p1_delta_[taggerID]->initValue(calib_p1.second); calib_p1_delta_[taggerID]->genValue(calib_p1.second); calib_p1_delta_[taggerID]->fixed(kTRUE); } } void LauFlavTag::cacheInputFitVars(LauFitDataTree* inputFitData, const TString& decayTimeVarName) { evtTagFlv_.clear(); evtMistag_.clear(); evtTrueTagFlv_.clear(); evtDecayFlv_.clear(); evtDecayTime_.clear(); // Loop over the taggers to check the branches for (std::size_t i=0; i < tagVarNames_.size(); ++i){ if ( not inputFitData->haveBranch( tagVarNames_[i] ) ) { std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << tagVarNames_[i] << "\"." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( not inputFitData->haveBranch( mistagVarNames_[i] ) ) { std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << mistagVarNames_[i] << "\"." << std::endl; gSystem->Exit(EXIT_FAILURE); } } if ( trueTagVarName_ != "" and not inputFitData->haveBranch( trueTagVarName_ ) ) { std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << trueTagVarName_ << "\"." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( decayFlvVarName_ != "" and not inputFitData->haveBranch( decayFlvVarName_ ) ) { std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << decayFlvVarName_ << "\"." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( decayTimeVarName != "" and not inputFitData->haveBranch( decayTimeVarName ) ) { std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << decayTimeVarName << "\"." << std::endl; gSystem->Exit(EXIT_FAILURE); } const std::size_t nEvents { inputFitData->nEvents() }; evtTagFlv_.reserve( nEvents ); evtMistag_.reserve( nEvents ); evtTrueTagFlv_.reserve( nEvents ); evtDecayFlv_.reserve( nEvents ); evtDecayTime_.reserve( nEvents ); LauFitData::const_iterator fitdata_iter; for (std::size_t iEvt = 0; iEvt < nEvents; iEvt++) { const LauFitData& dataValues = inputFitData->getData(iEvt); // For untagged events see if we have a truth tag for normalisation modes Int_t curEvtTrueTagFlv { ( trueTagVarName_ != "" ) ? static_cast( dataValues.at( trueTagVarName_ ) ) : 0 }; if ( curEvtTrueTagFlv > 1 ) { std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid true tag value " << curEvtTrueTagFlv << " for event " << iEvt << ", setting it to +1" << std::endl; curEvtTrueTagFlv = +1; } else if ( curEvtTrueTagFlv < -1 ){ std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid true tag value " << curEvtTrueTagFlv << " for event " << iEvt << ", setting it to -1" << std::endl; curEvtTrueTagFlv = -1; } curEvtTrueTagFlv_ = static_cast(curEvtTrueTagFlv); evtTrueTagFlv_.push_back(curEvtTrueTagFlv_); // Flavour at decay // TODO put this in a try catch block current error message is unhelpful if this throws Int_t curEvtDecayFlv { ( decayFlvVarName_ != "" ) ? static_cast( dataValues.at( decayFlvVarName_ ) ) : 0 }; if ( curEvtDecayFlv > 1 ) { std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid decay flavour value " << curEvtDecayFlv << " for event " << iEvt << ", setting it to +1" << std::endl; curEvtDecayFlv = +1; } else if ( curEvtDecayFlv < -1 ){ std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid decay flavour value " << curEvtDecayFlv << " for event " << iEvt << ", setting it to -1" << std::endl; curEvtDecayFlv = -1; } curEvtDecayFlv_ = static_cast(curEvtDecayFlv); evtDecayFlv_.push_back(curEvtDecayFlv_); for (std::size_t i=0; i < tagVarNames_.size(); ++i){ Int_t curEvtTagFlv { static_cast( dataValues.at( tagVarNames_[i] ) ) }; if ( curEvtTagFlv > 1 ) { std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv << " for event " << iEvt << ", setting it to +1" << std::endl; curEvtTagFlv = +1; } else if ( curEvtTagFlv < -1 ) { std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv << " for event " << iEvt << ", setting it to -1" << std::endl; curEvtTagFlv = -1; } curEvtTagFlv_[i] = static_cast( curEvtTagFlv ); curEvtMistag_[i] = static_cast( dataValues.at( mistagVarNames_[i] ) ); // Calibrated mistag > 0.5 is just a tag flip - handled automatically in getCapitalOmega function if (curEvtMistag_[i] > 0.5){ std::cerr<<"WARNING in LauFlavTag::cacheInputFitVars : Mistag value "<( dataValues.at( decayTimeVarName ) ); evtDecayTime_.push_back(curEvtDecayTime_); } } void LauFlavTag::updateEventInfo(const std::size_t iEvt) { //Assign current event variables curEvtTagFlv_ = evtTagFlv_[iEvt]; curEvtMistag_ = evtMistag_[iEvt]; curEvtTrueTagFlv_ = evtTrueTagFlv_[iEvt]; curEvtDecayFlv_ = evtDecayFlv_[iEvt]; curEvtDecayTime_ = evtDecayTime_[iEvt]; } void LauFlavTag::generateEventInfo(const Flavour trueTagFlv, const Double_t curEvtDecayTime) { curEvtTrueTagFlv_ = trueTagFlv; curEvtDecayFlv_ = Flavour::Unknown; curEvtDecayTime_ = curEvtDecayTime; Double_t randNo{0.0}; Double_t tagEffB0{0.0}; Double_t tagEffB0bar{0.0}; const std::size_t nTaggers { this->getNTaggers() }; for ( std::size_t taggerID{0}; taggerIDgetEtaGen(taggerID); if (this->getUseAveDelta()) { if (tagEff_ave_[taggerID]==nullptr){ const Double_t ave = tagEff_hist_ave_[taggerID]->GetBinContent(tagEff_hist_ave_[taggerID]->FindFixBin(curEvtDecayTime_)); const Double_t delta = tagEff_hist_delta_[taggerID]->GetBinContent(tagEff_hist_delta_[taggerID]->FindFixBin(curEvtDecayTime_)); tagEffB0 = ave + 0.5 * delta; tagEffB0bar = ave - 0.5 * delta; } else { tagEffB0 = tagEff_ave_[taggerID]->unblindValue() + 0.5 * tagEff_delta_[taggerID]->unblindValue(); tagEffB0bar = tagEff_ave_[taggerID]->unblindValue() - 0.5 * tagEff_delta_[taggerID]->unblindValue(); } } else { if (tagEff_B0_[taggerID]==nullptr){ tagEffB0 = tagEff_hist_B0_[taggerID]->GetBinContent(tagEff_hist_B0_[taggerID]->FindFixBin(curEvtDecayTime_)); tagEffB0bar = tagEff_hist_B0bar_[taggerID]->GetBinContent(tagEff_hist_B0bar_[taggerID]->FindFixBin(curEvtDecayTime_)); } else { tagEffB0 = tagEff_B0_[taggerID]->unblindValue(); tagEffB0bar = tagEff_B0bar_[taggerID]->unblindValue(); } } if (curEvtTrueTagFlv_ == Flavour::B) { randNo = LauRandom::randomFun()->Rndm(); // Try to tag in tageff% of cases if (randNo <= tagEffB0) { randNo = LauRandom::randomFun()->Rndm(); // Account for (calibrated) mistag if (randNo > this->getLittleOmega(taggerID, Flavour::B)){ curEvtTagFlv_[taggerID] = Flavour::B; } else { curEvtTagFlv_[taggerID] = Flavour::Bbar; } } else { curEvtTagFlv_[taggerID] = Flavour::Unknown; } } else if (curEvtTrueTagFlv_ == Flavour::Bbar) { randNo = LauRandom::randomFun()->Rndm(); // Try to tag in tageff% of cases if (randNo <= tagEffB0bar) { randNo = LauRandom::randomFun()->Rndm(); // Account for (calibrated) mistag if (randNo > this->getLittleOmega(taggerID, Flavour::Bbar)){ curEvtTagFlv_[taggerID] = Flavour::Bbar; } else { curEvtTagFlv_[taggerID] = Flavour::B; } } else { curEvtTagFlv_[taggerID] = Flavour::Unknown; } } else { std::cerr << "ERROR in LauFlavTag::generateEventInfo : Invalid true tag flavour, should be either B (+1) or Bbar (-1)" << std::endl; gSystem->Exit(EXIT_FAILURE); } } } void LauFlavTag::generateBkgndEventInfo(const std::size_t bkgndID, const Flavour trueTagFlv, const Flavour trueDecayFlv, const Double_t curEvtDecayTime) { if (bkgndID > bkgndIndex_.size()){ std::cerr << "ERROR in LauFlavTag::generateBkgndEventInfo : Invalid backgrond class identifier" << std::endl; gSystem->Exit(EXIT_FAILURE); } curEvtTrueTagFlv_ = trueTagFlv; curEvtDecayFlv_ = trueDecayFlv; curEvtDecayTime_ = curEvtDecayTime; Double_t randNo{0.0}; Double_t tagEffB0{0.0}; Double_t tagEffB0bar{0.0}; const std::size_t nTaggers { this->getNTaggers() }; for ( std::size_t taggerID{0}; taggerID( bkgndDecayFlvDep_[bkgndID] * curEvtDecayFlv_ + 1 ) }; this->getEtaGenBkgnd(taggerID,bkgndID,curEvtDecayFlv_); if (this->getUseAveDelta()) { if (tagEffBkgnd_ave_[bkgndID][taggerID][arrayIndex]==nullptr){ const Double_t ave = tagEffBkgnd_hist_ave_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_ave_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)); const Double_t delta = tagEffBkgnd_hist_delta_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_delta_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)); tagEffB0 = ave + 0.5 * delta; tagEffB0bar = ave - 0.5 * delta; } else { tagEffB0 = tagEffBkgnd_ave_[bkgndID][taggerID][arrayIndex]->unblindValue() + 0.5 * tagEffBkgnd_delta_[bkgndID][taggerID][arrayIndex]->unblindValue(); tagEffB0bar = tagEffBkgnd_ave_[bkgndID][taggerID][arrayIndex]->unblindValue() - 0.5 * tagEffBkgnd_delta_[bkgndID][taggerID][arrayIndex]->unblindValue(); } } else { if (tagEffBkgnd_B0_[bkgndID][taggerID][arrayIndex]==nullptr){ tagEffB0 = tagEffBkgnd_hist_B0_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_B0_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)); tagEffB0bar = tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)); } else { tagEffB0 = tagEffBkgnd_B0_[bkgndID][taggerID][arrayIndex]->unblindValue(); tagEffB0bar = tagEffBkgnd_B0bar_[bkgndID][taggerID][arrayIndex]->unblindValue(); } } if (bkgndTypes_[bkgndID] == BkgndType::Combinatorial){ randNo = LauRandom::randomFun()->Rndm(); if (randNo<=tagEffB0){ curEvtTagFlv_[taggerID] = Flavour::B; } else if(randNo<=(tagEffB0+tagEffB0bar)){ curEvtTagFlv_[taggerID] = Flavour::Bbar; } else { curEvtTagFlv_[taggerID] = Flavour::Unknown; } }else{ if (curEvtTrueTagFlv_ == Flavour::B) { randNo = LauRandom::randomFun()->Rndm(); // Try to tag in tageff% of cases if (randNo <= tagEffB0) { randNo = LauRandom::randomFun()->Rndm(); // Account for (calibrated) mistag if (randNo > this->getLittleOmegaBkgnd(bkgndID, taggerID, Flavour::B, curEvtDecayFlv_)){ curEvtTagFlv_[taggerID] = Flavour::B; } else { curEvtTagFlv_[taggerID] = Flavour::Bbar; } } else { curEvtTagFlv_[taggerID] = Flavour::Unknown; } } else if (curEvtTrueTagFlv_ == Flavour::Bbar) { randNo = LauRandom::randomFun()->Rndm(); // Try to tag in tageff% of cases if (randNo <= tagEffB0bar) { randNo = LauRandom::randomFun()->Rndm(); // Account for (calibrated) mistag if (randNo > this->getLittleOmegaBkgnd(bkgndID, taggerID, Flavour::Bbar, curEvtDecayFlv_)){ curEvtTagFlv_[taggerID] = Flavour::Bbar; } else { curEvtTagFlv_[taggerID] = Flavour::B; } } else { curEvtTagFlv_[taggerID] = Flavour::Unknown; } } else { std::cerr << "ERROR in LauFlavTag::generateBkgndEventInfo : Invalid true tag flavour, should be either B (+1) or Bbar (-1)" << std::endl; gSystem->Exit(EXIT_FAILURE); } } } } Double_t LauFlavTag::getLittleOmega(const std::size_t taggerID, const Flavour flag) const { if ( flag == Flavour::Unknown ){ std::cerr << "ERROR in LauFlavTag::getLittleOmega : Invalid flag, you must request either omega (+1) or omega bar (-1) to be returned" << std::endl; return 0.0; } Double_t calibp0(0.), calibp1(0.), calibp0bar(0.), calibp1bar(0.); //If we are floating average omega and delta omega we need to use those parameters instead if (useAveDelta_){ calibp0 = calib_p0_ave_[taggerID]->unblindValue() + 0.5*calib_p0_delta_[taggerID]->unblindValue(); calibp0bar = calib_p0_ave_[taggerID]->unblindValue() - 0.5*calib_p0_delta_[taggerID]->unblindValue(); calibp1 = calib_p1_ave_[taggerID]->unblindValue() + 0.5*calib_p1_delta_[taggerID]->unblindValue(); calibp1bar = calib_p1_ave_[taggerID]->unblindValue() - 0.5*calib_p1_delta_[taggerID]->unblindValue(); } else { calibp0 = calib_p0_B0_[taggerID]->unblindValue(); calibp0bar = calib_p0_B0bar_[taggerID]->unblindValue(); calibp1 = calib_p1_B0_[taggerID]->unblindValue(); calibp1bar = calib_p1_B0bar_[taggerID]->unblindValue(); } if ( flag == Flavour::B ){ return calibp0 + calibp1 * (curEvtMistag_[taggerID] - perEvtAvgMistag_[taggerID]); } else{ return calibp0bar + calibp1bar * (curEvtMistag_[taggerID] - perEvtAvgMistag_[taggerID]); } return 0.0; } Double_t LauFlavTag::getCapitalOmega(const std::size_t taggerID, const Flavour flag) const { if ( flag == Flavour::Unknown ){ std::cerr << "ERROR in LauFlavTag::getCapitalOmega : Invalid flag, you must request either Omega (+1) or Omega bar (-1) to be returned" << std::endl; return 0.0; } //Delta functions to control which terms contribute const bool delta_p1 { curEvtTagFlv_[taggerID] == Flavour::B }; const bool delta_m1 { curEvtTagFlv_[taggerID] == Flavour::Bbar }; const bool delta_0 { curEvtTagFlv_[taggerID] == Flavour::Unknown }; //Efficiency Double_t eff{0.0}; if (useAveDelta_){ if(flag==Flavour::B){ if (tagEff_ave_[taggerID]==nullptr){ eff = tagEff_hist_ave_[taggerID]->GetBinContent(tagEff_hist_ave_[taggerID]->FindFixBin(curEvtDecayTime_)) + 0.5*tagEff_hist_delta_[taggerID]->GetBinContent(tagEff_hist_delta_[taggerID]->FindFixBin(curEvtDecayTime_)); } else { eff = tagEff_ave_[taggerID]->unblindValue() + 0.5*tagEff_delta_[taggerID]->unblindValue(); } } else { if (tagEff_ave_[taggerID]==nullptr){ eff = tagEff_hist_ave_[taggerID]->GetBinContent(tagEff_hist_ave_[taggerID]->FindFixBin(curEvtDecayTime_)) - 0.5*tagEff_hist_delta_[taggerID]->GetBinContent(tagEff_hist_delta_[taggerID]->FindFixBin(curEvtDecayTime_)); } else { eff = tagEff_ave_[taggerID]->unblindValue() - 0.5*tagEff_delta_[taggerID]->unblindValue(); } } }else{ if(flag==Flavour::B){ if (tagEff_B0_[taggerID]==nullptr){ eff = tagEff_hist_B0_[taggerID]->GetBinContent(tagEff_hist_B0_[taggerID]->FindFixBin(curEvtDecayTime_)); } else { eff = tagEff_B0_[taggerID]->unblindValue(); } }else{ if (tagEff_B0bar_[taggerID]==nullptr){ eff = tagEff_hist_B0bar_[taggerID]->GetBinContent(tagEff_hist_B0bar_[taggerID]->FindFixBin(curEvtDecayTime_)); } else { eff = tagEff_B0bar_[taggerID]->unblindValue(); } } } //Little omega const Double_t omega { this->getLittleOmega(taggerID, flag) }; //Transform to omega prime - TODO isn't this the inverse, getLittleOmega is actually giving us omegaPrime and on the next line we convert back to omega? Double_t omegaPrime { useEtaPrime_ ? (1.0/(1.0+TMath::Exp(-1.0*omega))) : omega }; //little omega must be between 0 and 1. Force this for now, if the fits keep getting stuck can look more closely at it. static Bool_t tooSmallWarningIssued {kFALSE}; static Bool_t tooLargeWarningIssued {kFALSE}; if (omegaPrime < 0.0){ if ( not tooSmallWarningIssued ) { std::cerr << "WARNING in LauFlavTag::getCapitalOmega : The value of little omega is less than 0, shifting to 0\n"; std::cerr << " : Further WARNINGs of this type will be suppressed" << std::endl; tooSmallWarningIssued = kTRUE; } omegaPrime = 0.0; } if (omegaPrime > 1.0){ if ( not tooLargeWarningIssued ) { std::cerr << "WARNING in LauFlavTag::getCapitalOmega : The value of little omega is greater than 1, shifting to 1\n"; std::cerr << " : Further WARNINGs of this type will be suppressed" << std::endl; tooLargeWarningIssued = kTRUE; } omegaPrime = 1.0; } //eta PDF value const std::vector abs { curEvtMistag_[taggerID] }; etaPdfs_[taggerID]->calcLikelihoodInfo(abs); const Double_t h { etaPdfs_[taggerID]->getLikelihood() }; const Double_t u { 2.0 }; // the PDF value for a uniform PDF between 0.0 and 0.5 //If h returns 0 for a tagged event, the event likelihood will be zero, so we should warn that this is where the problem lies if (h==0.0 && not delta_0){ std::cerr << "WARNING in LauFlavTag::getCapitalOmega : The value of the signal eta PDF is zero at eta = " << curEvtMistag_[taggerID] << std::endl; return 0.0; } //Put it together if (flag == Flavour::B){ return (delta_p1*eff*(1.0-omegaPrime) + delta_m1*eff*omegaPrime)*h + delta_0*(1.0-eff)*u; } else { return (delta_m1*eff*(1.0-omegaPrime) + delta_p1*eff*omegaPrime)*h + delta_0*(1.0-eff)*u; } } Double_t LauFlavTag::getLittleOmegaBkgnd(const std::size_t bkgndID, const std::size_t taggerID, const Flavour flag, const Flavour decayFlv) const { if ( flag == Flavour::Unknown ){ std::cerr << "ERROR in LauFlavTag::getLittleOmegaBkgnd : Invalid flag, you must request either omega (+1) or omega bar (-1) to be returned" << std::endl; return 0.0; } Double_t calibp0(0.), calibp1(0.), calibp0bar(0.), calibp1bar(0.); //If we are floating average omega and delta omega we need to use those parameters instead if (useAveDelta_){ calibp0 = calib_p0_ave_[taggerID]->unblindValue() + 0.5*calib_p0_delta_[taggerID]->unblindValue(); calibp0bar = calib_p0_ave_[taggerID]->unblindValue() - 0.5*calib_p0_delta_[taggerID]->unblindValue(); calibp1 = calib_p1_ave_[taggerID]->unblindValue() + 0.5*calib_p1_delta_[taggerID]->unblindValue(); calibp1bar = calib_p1_ave_[taggerID]->unblindValue() - 0.5*calib_p1_delta_[taggerID]->unblindValue(); } else { calibp0 = calib_p0_B0_[taggerID]->unblindValue(); calibp0bar = calib_p0_B0bar_[taggerID]->unblindValue(); calibp1 = calib_p1_B0_[taggerID]->unblindValue(); calibp1bar = calib_p1_B0bar_[taggerID]->unblindValue(); } const std::size_t arrayIndex { static_cast( bkgndDecayFlvDep_[bkgndID] * decayFlv + 1 ) }; if ( flag == Flavour::B ){ return calibp0 + calibp1 * (curEvtMistag_[taggerID] - avgMistagBkgnd_[bkgndID][taggerID][arrayIndex]); } else{ return calibp0bar + calibp1bar * (curEvtMistag_[taggerID] - avgMistagBkgnd_[bkgndID][taggerID][arrayIndex]); } return 0.0; } Double_t LauFlavTag::getCapitalOmegaBkgnd(const std::size_t bkgndID, const std::size_t taggerID, const Flavour flag, const Flavour decayFlv) const { //Delta functions to control which terms contribute const bool delta_p1 { curEvtTagFlv_[taggerID] == Flavour::B }; const bool delta_m1 { curEvtTagFlv_[taggerID] == Flavour::Bbar }; const bool delta_0 { curEvtTagFlv_[taggerID] == Flavour::Unknown }; //Use the right histograms/parameters depending on the decay flavour const std::size_t arrayIndex { static_cast( bkgndDecayFlvDep_[bkgndID] * decayFlv + 1 ) }; //Efficiency Double_t effB0{0.0}, effB0bar{0.0}; if (useAveDelta_){ if (tagEffBkgnd_ave_[bkgndID][taggerID][arrayIndex]==nullptr){ effB0 = tagEffBkgnd_hist_ave_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_ave_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)) + 0.5*tagEffBkgnd_hist_delta_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_delta_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)); effB0bar = tagEffBkgnd_hist_ave_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_ave_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)) - 0.5*tagEffBkgnd_hist_delta_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_delta_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)); } else { effB0 = tagEffBkgnd_ave_[bkgndID][taggerID][arrayIndex]->unblindValue() + 0.5*tagEffBkgnd_delta_[bkgndID][taggerID][arrayIndex]->unblindValue(); effB0bar = tagEffBkgnd_ave_[bkgndID][taggerID][arrayIndex]->unblindValue() - 0.5*tagEffBkgnd_delta_[bkgndID][taggerID][arrayIndex]->unblindValue(); } }else{ if (tagEffBkgnd_B0_[bkgndID][taggerID][arrayIndex]==nullptr){ effB0 = tagEffBkgnd_hist_B0_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_B0_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)); effB0bar = tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][arrayIndex]->GetBinContent(tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][arrayIndex]->FindFixBin(curEvtDecayTime_)); } else { effB0 = tagEffBkgnd_B0_[bkgndID][taggerID][arrayIndex]->unblindValue(); effB0bar = tagEffBkgnd_B0bar_[bkgndID][taggerID][arrayIndex]->unblindValue(); } } //Need to know which background eta PDF to use - bkgndID const std::vector abs { curEvtMistag_[taggerID] }; etaBkgndPdfs_[bkgndID][taggerID][arrayIndex]->calcLikelihoodInfo(abs); const Double_t h { etaBkgndPdfs_[bkgndID][taggerID][arrayIndex]->getLikelihood() }; const Double_t u { 2.0 }; // the PDF value for a uniform PDF between 0.0 and 0.5 switch ( bkgndTypes_[bkgndID] ) { case BkgndType::Combinatorial : // Combinatorial efficiences are defined differently // effB0 - chance to tag any combinatorial event as a B0 // effB0bar - chance to tag any combinatorial event as a B0bar return (delta_p1*effB0 + delta_m1*effB0bar)*h + delta_0*(1.0 - effB0 - effB0bar)*u; case BkgndType::FlavourSpecific : case BkgndType::SelfConjugate : case BkgndType::NonSelfConjugate : { if ( flag == Flavour::Unknown ){ std::cerr << "ERROR in LauFlavTag::getCapitalOmegaBkgnd : Invalid flag, you must request either Omega (+1) or Omega bar (-1) to be returned" << std::endl; return 0.0; } const Double_t omega { this->getLittleOmegaBkgnd(bkgndID, taggerID, flag, decayFlv) }; //Transform to omega prime - TODO isn't this the inverse, getLittleOmega is actually giving us omegaPrime and on the next line we convert back to omega? Double_t omegaPrime { useEtaPrime_ ? (1.0/(1.0+TMath::Exp(-1.0*omega))) : omega }; //little omega must be between 0 and 1. Force this for now, if the fits keep getting stuck can look more closely at it. static Bool_t tooSmallWarningIssued {kFALSE}; static Bool_t tooLargeWarningIssued {kFALSE}; if (omegaPrime < 0.0){ if ( not tooSmallWarningIssued ) { std::cerr << "WARNING in LauFlavTag::getCapitalOmegaBkgnd : The value of little omega is less than 0, shifting to 0\n"; std::cerr << " : Further WARNINGs of this type will be suppressed" << std::endl; tooSmallWarningIssued = kTRUE; } omegaPrime = 0.0; } if (omegaPrime > 1.0){ if ( not tooLargeWarningIssued ) { std::cerr << "WARNING in LauFlavTag::getCapitalOmegaBkgnd : The value of little omega is greater than 1, shifting to 1\n"; std::cerr << " : Further WARNINGs of this type will be suppressed" << std::endl; tooLargeWarningIssued = kTRUE; } omegaPrime = 1.0; } if (flag == Flavour::B){ return (delta_p1*effB0*(1.0-omegaPrime) + delta_m1*effB0*omegaPrime)*h + delta_0*(1.0-effB0)*u; } else { return (delta_m1*effB0bar*(1.0-omegaPrime) + delta_p1*effB0bar*omegaPrime)*h + delta_0*(1.0-effB0bar)*u; } } } //Should never get here, but it keeps GCC happy return 0.0; } void LauFlavTag::setBkgndParams(const TString& bkgndName, const TString& taggerName, LauAbsPdf* etaPdf, const std::pair tagEff) { if ( bkgndIndex_.count(bkgndName) == 0 ) { std::cerr << "ERROR in LauFlavTag::setBkgndParams : Background name \"" << bkgndName << "\" not recognised please check your options" << std::endl; return; } if ( taggerIndex_.count(taggerName) == 0 ){ std::cerr << "ERROR in LauFlavTag::setBkgndParams : Tagger name \"" << taggerName << "\" not recognised please check your options" << std::endl; return; } const std::size_t bkgndID { bkgndIndex_.at(bkgndName) }; const std::size_t taggerID { taggerIndex_.at(taggerName) }; bkgndDecayFlvDep_[bkgndID] = false; etaBkgndPdfs_[bkgndID][taggerID][0] = nullptr; etaBkgndPdfs_[bkgndID][taggerID][1] = etaPdf; etaBkgndPdfs_[bkgndID][taggerID][2] = nullptr; Lau1DHistPdf* etahistpdf = dynamic_cast(etaPdf); if (etahistpdf){ avgMistagBkgnd_[bkgndID][taggerID][0] = 0.0; avgMistagBkgnd_[bkgndID][taggerID][1] = etahistpdf->getMean(); avgMistagBkgnd_[bkgndID][taggerID][2] = 0.0; } else { std::cerr << "WARNING in LauFlavTag::setBkgndParams : Couldn't determine average eta value from PDF. Setting it to 0.4." << std::endl; avgMistagBkgnd_[bkgndID][taggerID][0] = 0.0; avgMistagBkgnd_[bkgndID][taggerID][1] = 0.4; avgMistagBkgnd_[bkgndID][taggerID][2] = 0.0; } if (useAveDelta_){ const TString tagEff_aveName{"tagEff_ave_"+taggerName+"_bkgnd_"+bkgndName}; const TString tagEff_deltaName{"tagEff_delta_"+taggerName+"_bkgnd_"+bkgndName}; tagEffBkgnd_ave_[bkgndID][taggerID][0] = nullptr; tagEffBkgnd_delta_[bkgndID][taggerID][0] = nullptr; tagEffBkgnd_ave_[bkgndID][taggerID][1] = new LauParameter{tagEff_aveName, tagEff.first, 0.0, 1.0, kTRUE}; tagEffBkgnd_delta_[bkgndID][taggerID][1] = new LauParameter{tagEff_deltaName, tagEff.second, -1.0, 1.0, kTRUE}; tagEffBkgnd_ave_[bkgndID][taggerID][2] = nullptr; tagEffBkgnd_delta_[bkgndID][taggerID][2] = nullptr; } else { const TString tagEff_b0Name{"tagEff_b0_"+taggerName+"_bkgnd_"+bkgndName}; const TString tagEff_b0barName{"tagEff_b0bar_"+taggerName+"_bkgnd_"+bkgndName}; tagEffBkgnd_B0_[bkgndID][taggerID][0] = nullptr; tagEffBkgnd_B0bar_[bkgndID][taggerID][0] = nullptr; tagEffBkgnd_B0_[bkgndID][taggerID][1] = new LauParameter{tagEff_b0Name, tagEff.first, 0.0, 1.0, kTRUE}; tagEffBkgnd_B0bar_[bkgndID][taggerID][1] = new LauParameter{tagEff_b0barName, tagEff.second, 0.0, 1.0, kTRUE}; tagEffBkgnd_B0_[bkgndID][taggerID][2] = nullptr; tagEffBkgnd_B0bar_[bkgndID][taggerID][2] = nullptr; } std::cout << "INFO in LauFlavTag::setBkgndParams : Added efficiency parameters and eta PDF for background " << bkgndName << " for tagger " << taggerName << std::endl; } void LauFlavTag::setBkgndParams(const TString& bkgndName, const TString& taggerName, LauAbsPdf* etaPdfB, std::pair tagEffB, LauAbsPdf* etaPdfBbar, std::pair tagEffBbar) { if ( bkgndIndex_.count(bkgndName) == 0 ) { std::cerr << "ERROR in LauFlavTag::setBkgndParams : Background name \"" << bkgndName << "\" not recognised please check your options" << std::endl; return; } if ( taggerIndex_.count(taggerName) == 0 ){ std::cerr << "ERROR in LauFlavTag::setBkgndParams : Tagger name \"" << taggerName << "\" not recognised please check your options" << std::endl; return; } const std::size_t bkgndID { bkgndIndex_.at(bkgndName) }; const std::size_t taggerID { taggerIndex_.at(taggerName) }; bkgndDecayFlvDep_[bkgndID] = true; etaBkgndPdfs_[bkgndID][taggerID][0] = etaPdfBbar; Lau1DHistPdf* etahistpdf = dynamic_cast(etaPdfBbar); if (etahistpdf){ avgMistagBkgnd_[bkgndID][taggerID][0] = etahistpdf->getMean(); } else { std::cerr << "WARNING in LauFlavTag::setBkgndParams : Couldn't determine average eta value from PDF. Setting it to 0.4." << std::endl; avgMistagBkgnd_[bkgndID][taggerID][0] = 0.4; } etaBkgndPdfs_[bkgndID][taggerID][1] = nullptr; avgMistagBkgnd_[bkgndID][taggerID][1] = 0.0; etaBkgndPdfs_[bkgndID][taggerID][2] = etaPdfB; etahistpdf = dynamic_cast(etaPdfB); if (etahistpdf){ avgMistagBkgnd_[bkgndID][taggerID][2] = etahistpdf->getMean(); } else { std::cerr << "WARNING in LauFlavTag::setBkgndParams : Couldn't determine average eta value from PDF. Setting it to 0.4." << std::endl; avgMistagBkgnd_[bkgndID][taggerID][2] = 0.4; } if (useAveDelta_){ TString tagEff_aveName { "tagEff_ave_decayFlvB0bar_"+taggerName+"_bkgnd_"+bkgndName }; TString tagEff_deltaName { "tagEff_delta_decayFlvB0bar_"+taggerName+"_bkgnd_"+bkgndName }; tagEffBkgnd_ave_[bkgndID][taggerID][0] = new LauParameter{tagEff_aveName, tagEffBbar.first, 0.0, 1.0, kTRUE}; tagEffBkgnd_delta_[bkgndID][taggerID][0] = new LauParameter{tagEff_deltaName, tagEffBbar.second, -1.0, 1.0, kTRUE}; tagEffBkgnd_ave_[bkgndID][taggerID][1] = nullptr; tagEffBkgnd_delta_[bkgndID][taggerID][1] = nullptr; tagEff_aveName = "tagEff_ave_decayFlvB0_"+taggerName+"_bkgnd_"+bkgndName; tagEff_deltaName = "tagEff_delta_decayFlvB0_"+taggerName+"_bkgnd_"+bkgndName; tagEffBkgnd_ave_[bkgndID][taggerID][2] = new LauParameter{tagEff_aveName, tagEffB.first, 0.0, 1.0, kTRUE}; tagEffBkgnd_delta_[bkgndID][taggerID][2] = new LauParameter{tagEff_deltaName, tagEffB.second, -1.0, 1.0, kTRUE}; } else { TString tagEff_b0Name { "tagEff_b0_decayFlvB0bar_"+taggerName+"_bkgnd_"+bkgndName }; TString tagEff_b0barName { "tagEff_b0bar_decayFlvB0bar_"+taggerName+"_bkgnd_"+bkgndName }; tagEffBkgnd_B0_[bkgndID][taggerID][0] = new LauParameter{tagEff_b0Name, tagEffBbar.first, 0.0, 1.0, kTRUE}; tagEffBkgnd_B0bar_[bkgndID][taggerID][0] = new LauParameter{tagEff_b0barName, tagEffBbar.second, 0.0, 1.0, kTRUE}; tagEffBkgnd_B0_[bkgndID][taggerID][1] = nullptr; tagEffBkgnd_B0bar_[bkgndID][taggerID][1] = nullptr; tagEff_b0Name = "tagEff_b0_decayFlvB0_"+taggerName+"_bkgnd_"+bkgndName; tagEff_b0barName = "tagEff_b0bar_decayFlvB0_"+taggerName+"_bkgnd_"+bkgndName; tagEffBkgnd_B0_[bkgndID][taggerID][2] = new LauParameter{tagEff_b0Name, tagEffB.first, 0.0, 1.0, kTRUE}; tagEffBkgnd_B0bar_[bkgndID][taggerID][2] = new LauParameter{tagEff_b0barName, tagEffB.second, 0.0, 1.0, kTRUE}; } std::cout << "INFO in LauFlavTag::setBkgndParams : Added efficiency parameters and eta PDF for background " << bkgndName << " for tagger " << taggerName << std::endl; } void LauFlavTag::setBkgndParams(const TString& bkgndName, const TString& taggerName, LauAbsPdf* etaPdf, std::pair tagEff) { if ( bkgndIndex_.count(bkgndName) == 0 ) { std::cerr << "ERROR in LauFlavTag::setBkgndParams : Background name \"" << bkgndName << "\" not recognised please check your options" << std::endl; return; } if (taggerIndex_.count(taggerName)==0){ std::cerr << "ERROR in LauFlavTag::setBkgndParams : Tagger name \"" << taggerName << "\" not recognised please check your options" << std::endl; return; } if (tagEff.first==nullptr || tagEff.second==nullptr){ std::cerr << "ERROR in LauFlavTag::setBkgndParams : Efficiency histogram(s) do not exist please check your options" << std::endl; return; } const std::size_t bkgndID { bkgndIndex_.at(bkgndName) }; const std::size_t taggerID { taggerIndex_.at(taggerName) }; bkgndDecayFlvDep_[bkgndID] = false; etaBkgndPdfs_[bkgndID][taggerID][0] = nullptr; etaBkgndPdfs_[bkgndID][taggerID][1] = etaPdf; etaBkgndPdfs_[bkgndID][taggerID][2] = nullptr; Lau1DHistPdf* etahistpdf = dynamic_cast(etaPdf); if (etahistpdf){ avgMistagBkgnd_[bkgndID][taggerID][0] = 0.0; avgMistagBkgnd_[bkgndID][taggerID][1] = etahistpdf->getMean(); avgMistagBkgnd_[bkgndID][taggerID][2] = 0.0; } else { std::cerr << "WARNING in LauFlavTag::setBkgndParams : Couldn't determine average eta value from PDF. Setting it to 0.4." << std::endl; avgMistagBkgnd_[bkgndID][taggerID][0] = 0.0; avgMistagBkgnd_[bkgndID][taggerID][1] = 0.4; avgMistagBkgnd_[bkgndID][taggerID][2] = 0.0; } if (useAveDelta_){ tagEffBkgnd_hist_ave_[bkgndID][taggerID][0] = nullptr; tagEffBkgnd_hist_delta_[bkgndID][taggerID][0] = nullptr; tagEffBkgnd_hist_ave_[bkgndID][taggerID][1] = tagEff.first; tagEffBkgnd_hist_delta_[bkgndID][taggerID][1] = tagEff.second; tagEffBkgnd_hist_ave_[bkgndID][taggerID][2] = nullptr; tagEffBkgnd_hist_delta_[bkgndID][taggerID][2] = nullptr; } else { tagEffBkgnd_hist_B0_[bkgndID][taggerID][0] = nullptr; tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][0] = nullptr; tagEffBkgnd_hist_B0_[bkgndID][taggerID][1] = tagEff.first; tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][1] = tagEff.second; tagEffBkgnd_hist_B0_[bkgndID][taggerID][2] = nullptr; tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][2] = nullptr; } std::cout << "INFO in LauFlavTag::setBkgndParams : Added efficiency histograms and eta PDF for background " << bkgndName << " for tagger " << taggerName << std::endl; } void LauFlavTag::setBkgndParams(const TString& bkgndName, const TString& taggerName, LauAbsPdf* etaPdfB, std::pair tagEffB, LauAbsPdf* etaPdfBbar, std::pair tagEffBbar) { if ( bkgndIndex_.count(bkgndName) == 0 ) { std::cerr << "ERROR in LauFlavTag::setBkgndParams : Background name \"" << bkgndName << "\" not recognised please check your options" << std::endl; return; } if (taggerIndex_.count(taggerName)==0){ std::cerr << "ERROR in LauFlavTag::setBkgndParams : Tagger name \"" << taggerName << "\" not recognised please check your options" << std::endl; return; } if (tagEffB.first==nullptr || tagEffB.second==nullptr || tagEffBbar.first==nullptr || tagEffBbar.second==nullptr){ std::cerr << "ERROR in LauFlavTag::setBkgndParams : Efficiency histogram(s) do not exist please check your options" << std::endl; return; } const std::size_t bkgndID { bkgndIndex_.at(bkgndName) }; const std::size_t taggerID { taggerIndex_.at(taggerName) }; bkgndDecayFlvDep_[bkgndID] = true; etaBkgndPdfs_[bkgndID][taggerID][0] = etaPdfBbar; Lau1DHistPdf* etahistpdf = dynamic_cast(etaPdfBbar); if (etahistpdf){ avgMistagBkgnd_[bkgndID][taggerID][0] = etahistpdf->getMean(); } else { std::cerr << "WARNING in LauFlavTag::setBkgndParams : Couldn't determine average eta value from PDF. Setting it to 0.4." << std::endl; avgMistagBkgnd_[bkgndID][taggerID][0] = 0.4; } etaBkgndPdfs_[bkgndID][taggerID][1] = nullptr; avgMistagBkgnd_[bkgndID][taggerID][1] = 0.0; etaBkgndPdfs_[bkgndID][taggerID][2] = etaPdfB; etahistpdf = dynamic_cast(etaPdfB); if (etahistpdf){ avgMistagBkgnd_[bkgndID][taggerID][2] = etahistpdf->getMean(); } else { std::cerr << "WARNING in LauFlavTag::setBkgndParams : Couldn't determine average eta value from PDF. Setting it to 0.4." << std::endl; avgMistagBkgnd_[bkgndID][taggerID][2] = 0.4; } if (useAveDelta_){ tagEffBkgnd_hist_ave_[bkgndID][taggerID][0] = tagEffBbar.first; tagEffBkgnd_hist_delta_[bkgndID][taggerID][0] = tagEffBbar.second; tagEffBkgnd_hist_ave_[bkgndID][taggerID][1] = nullptr; tagEffBkgnd_hist_delta_[bkgndID][taggerID][1] = nullptr; tagEffBkgnd_hist_ave_[bkgndID][taggerID][2] = tagEffB.first; tagEffBkgnd_hist_delta_[bkgndID][taggerID][2] = tagEffB.second; } else { tagEffBkgnd_hist_B0_[bkgndID][taggerID][0] = tagEffBbar.first; tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][0] = tagEffBbar.second; tagEffBkgnd_hist_B0_[bkgndID][taggerID][1] = nullptr; tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][1] = nullptr; tagEffBkgnd_hist_B0_[bkgndID][taggerID][2] = tagEffB.first; tagEffBkgnd_hist_B0bar_[bkgndID][taggerID][2] = tagEffB.second; } std::cout << "INFO in LauFlavTag::setBkgndParams : Added efficiency histograms and eta PDF for background " << bkgndName << " for tagger " << taggerName << std::endl; } Double_t LauFlavTag::getEtaGen(const std::size_t taggerID) { const LauFitData data { etaPdfs_[taggerID]->generate(nullptr) }; Double_t etagen { data.at(etaPdfs_[taggerID]->varName()) }; if (etagen > 0.5){etagen = 0.5;} if (etagen < 0.0){etagen = 0.0;} curEvtMistag_[taggerID] = etagen; return etagen; } Double_t LauFlavTag::getEtaGenBkgnd(const std::size_t taggerID, const std::size_t bkgndID, const Flavour decayFlv) { const std::size_t arrayIndex { static_cast( bkgndDecayFlvDep_[bkgndID] * decayFlv + 1 ) }; const LauFitData data { etaBkgndPdfs_[bkgndID][taggerID][arrayIndex]->generate(nullptr) }; Double_t etagen { data.at(etaBkgndPdfs_[bkgndID][taggerID][arrayIndex]->varName()) }; if (etagen > 0.5){etagen = 0.5;} if (etagen < 0.0){etagen = 0.0;} curEvtMistag_[taggerID] = etagen; return etagen; } void LauFlavTag::setTrueTagVarName(TString trueTagVarName) { trueTagVarName_ = std::move(trueTagVarName); } void LauFlavTag::setDecayFlvVarName(TString decayFlvVarName) { decayFlvVarName_ = std::move(decayFlvVarName); } void LauFlavTag::addP0GaussianConstraints(const TString& taggerName, const std::pair constraint1, const std::pair constraint2) { // Does key exist? if (taggerIndex_.count(taggerName)==0){ std::cerr << "ERROR in LauFlavTag::addP0GaussianConstraints : Tagger name not recognised please check your options" << std::endl; std::cerr << "ERROR in LauFlavTag::addP0GaussianConstraints : Constraints have not been applied" << std::endl; return; } // Find index in the vectors from the tagger name const std::size_t pos { taggerIndex_.at(taggerName) }; if (!useAveDelta_){ calib_p0_B0_[pos]->addGaussianConstraint(constraint1.first,constraint1.second); calib_p0_B0bar_[pos]->addGaussianConstraint(constraint2.first,constraint2.second); }else{ calib_p0_ave_[pos]->addGaussianConstraint(constraint1.first,constraint1.second); calib_p0_delta_[pos]->addGaussianConstraint(constraint2.first,constraint2.second); } std::cout << "INFO in LauFlavTag::addP0GaussianConstraints : Added Gaussian constraints for the P0 calibration parameters of tagger " << taggerName << std::endl; } void LauFlavTag::addP1GaussianConstraints(const TString& taggerName, const std::pair constraint1, const std::pair constraint2) { // Does key exist? if (taggerIndex_.count(taggerName)==0){ std::cerr << "ERROR in LauFlavTag::addP1GaussianConstraints : Tagger name not recognised please check your options" << std::endl; std::cerr << "ERROR in LauFlavTag::addP1GaussianConstraints : Constraints have not been applied" << std::endl; return; } // Find index in the vector from the tagger name const std::size_t pos { taggerIndex_.at(taggerName) }; if (!useAveDelta_){ calib_p1_B0_[pos]->addGaussianConstraint(constraint1.first,constraint1.second); calib_p1_B0bar_[pos]->addGaussianConstraint(constraint2.first,constraint2.second); }else{ calib_p1_ave_[pos]->addGaussianConstraint(constraint1.first,constraint1.second); calib_p1_delta_[pos]->addGaussianConstraint(constraint2.first,constraint2.second); } std::cout << "INFO in LauFlavTag::addP1GaussianConstraints : Added Gaussian constraints for the P1 calibration parameters of tagger " << taggerName << std::endl; } void LauFlavTag::addTagEffGaussianConstraints(const TString& taggerName, const std::pair constraint1, const std::pair constraint2) { // Does key exist? if (taggerIndex_.count(taggerName)==0){ std::cerr << "ERROR in LauFlavTag::addTagEffGaussianConstraints : Tagger name not recognised please check your options" << std::endl; std::cerr << "ERROR in LauFlavTag::addTagEffGaussianConstraints : Constraints have not been applied" << std::endl; return; } // Find index in the vector from the tagger name const std::size_t pos { taggerIndex_.at(taggerName) }; if (tagEff_B0_[pos] == nullptr){ std::cerr << "ERROR in LauFlavTag::addTagEffGaussianConstraints : Cannot add Gaussian constraints to a histogram!" << std::endl; return; } if (!useAveDelta_){ tagEff_B0_[pos]->addGaussianConstraint(constraint1.first,constraint1.second); tagEff_B0bar_[pos]->addGaussianConstraint(constraint2.first,constraint2.second); }else{ tagEff_ave_[pos]->addGaussianConstraint(constraint1.first,constraint1.second); tagEff_delta_[pos]->addGaussianConstraint(constraint2.first,constraint2.second); } std::cout << "INFO in LauFlavTag::addTagEffGaussianConstraints : Added Gaussian constraints for the tagging efficiency parameters of tagger " << taggerName << std::endl; } LauParameter* LauFlavTag::findParameter( const TString& taggerName, std::vector& parameters ) { // Check the tagger name is valid auto iter = taggerIndex_.find(taggerName); if ( iter == taggerIndex_.end() ){ return nullptr; } // If so, find the appropriate parameter const std::size_t index { iter->second }; return parameters[index]; } void LauFlavTag::floatCalibParP0B0(const TString& taggerName) { if (useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0B0 : Trying to set calibration parameters for B0/B0bar not average/delta" << std::endl; return; } if (taggerName==""){ // Float parameters for all taggers for (auto& param : calib_p0_B0_){ if (param==nullptr){continue;} param->fixed(kFALSE); } return; } // Float parameter for requested tagger LauParameter* par { this->findParameter( taggerName, calib_p0_B0_ ) }; if ( not par ){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0B0 : Tagger name not recognised please check your options" << std::endl; return; } par->fixed(kFALSE); } void LauFlavTag::floatCalibParP1B0(const TString& taggerName) { if (useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1B0 : Trying to set calibration parameters for B0/B0bar not average/delta" << std::endl; return; } if (taggerName==""){ // Float parameters for all taggers for (auto& param : calib_p1_B0_){ if (param==nullptr){continue;} param->fixed(kFALSE); } return; } // Float parameter for requested tagger LauParameter* par { this->findParameter( taggerName, calib_p1_B0_ ) }; if ( not par ){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1B0 : Tagger name not recognised please check your options" << std::endl; return; } par->fixed(kFALSE); } void LauFlavTag::floatCalibParP0B0bar(const TString& taggerName) { if (useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0B0bar : Trying to set calibration parameters for B0/B0bar not average/delta" << std::endl; return; } if (taggerName==""){ // Float parameters for all taggers for (auto& param : calib_p0_B0bar_){ if (param==nullptr){continue;} param->fixed(kFALSE); } return; } // Float parameter for requested tagger LauParameter* par { this->findParameter( taggerName, calib_p0_B0bar_ ) }; if ( not par ) { std::cerr << "ERROR in LauFlavTag::floatCalibParP0B0bar : Tagger name not recognised please check your options" << std::endl; return; } par->fixed(kFALSE); } void LauFlavTag::floatCalibParP1B0bar(const TString& taggerName) { if (useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1B0bar : Trying to set calibration parameters for B0/B0bar not average/delta" << std::endl; return; } if (taggerName==""){ // Float parameters for all taggers for (auto& param : calib_p1_B0bar_){ if (param==nullptr){continue;} param->fixed(kFALSE); } return; } // Float parameter for requested tagger LauParameter* par { this->findParameter( taggerName, calib_p1_B0bar_ ) }; if ( not par ) { std::cerr << "ERROR in LauFlavTag::floatCalibParP1B0bar : Tagger name not recognised please check your options" << std::endl; return; } par->fixed(kFALSE); } void LauFlavTag::floatCalibParP0Ave(const TString& taggerName) { if (!useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0Ave : Trying to set calibration parameters for average/delta not B0/B0bar" << std::endl; return; } if (taggerName==""){ // Float parameters for all taggers for (auto& param : calib_p0_ave_){ if (param==nullptr){continue;} param->fixed(kFALSE); } return; } // Float parameter for requested tagger LauParameter* par { this->findParameter( taggerName, calib_p0_ave_ ) }; if ( not par ) { std::cerr << "ERROR in LauFlavTag::floatCalibParP0Ave : Tagger name not recognised please check your options" << std::endl; return; } par->fixed(kFALSE); } void LauFlavTag::floatCalibParP0Delta(const TString& taggerName) { if (!useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0Delta : Trying to set calibration parameters for average/delta not B0/B0bar" << std::endl; return; } if (taggerName==""){ // Float parameters for all taggers for (auto& param : calib_p0_delta_){ if (param==nullptr){continue;} param->fixed(kFALSE); } return; } // Float parameter for requested tagger LauParameter* par { this->findParameter( taggerName, calib_p0_delta_ ) }; if ( not par ) { std::cerr << "ERROR in LauFlavTag::floatCalibParP0Delta : Tagger name not recognised please check your options" << std::endl; return; } par->fixed(kFALSE); } void LauFlavTag::floatCalibParP1Ave(const TString& taggerName) { if (!useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1Ave : Trying to set calibration parameters for average/delta not B0/B0bar" << std::endl; return; } if (taggerName==""){ // Float parameters for all taggers for (auto& param : calib_p1_ave_){ if (param==nullptr){continue;} param->fixed(kFALSE); } return; } // Float parameter for requested tagger LauParameter* par { this->findParameter( taggerName, calib_p1_ave_ ) }; if ( not par ) { std::cerr << "ERROR in LauFlavTag::floatCalibParP1Ave : Tagger name not recognised please check your options" << std::endl; return; } par->fixed(kFALSE); } void LauFlavTag::floatCalibParP1Delta(const TString& taggerName) { if (!useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1Delta : Trying to set calibration parameters for average/delta not B0/B0bar" << std::endl; return; } if (taggerName==""){ // Float parameters for all taggers for (auto& param : calib_p1_delta_){ if (param==nullptr){continue;} param->fixed(kFALSE); } return; } // Float parameter for requested tagger LauParameter* par { this->findParameter( taggerName, calib_p1_delta_ ) }; if ( not par ) { std::cerr << "ERROR in LauFlavTag::floatCalibParP1Delta : Tagger name not recognised please check your options" << std::endl; return; } par->fixed(kFALSE); } -void LauFlavTag::floatAllCalibPars(){ +void LauFlavTag::floatAllCalibPars() +{ if (useAveDelta_){ this->floatCalibParP0Ave(); this->floatCalibParP0Delta(); this->floatCalibParP1Ave(); this->floatCalibParP1Delta(); } else { this->floatCalibParP0B0(); this->floatCalibParP1B0(); this->floatCalibParP0B0bar(); this->floatCalibParP1B0bar(); } } + +void LauFlavTag::updatePulls() +{ + const std::size_t nTaggers { this->getNTaggers() }; + if (useAveDelta_){ + for ( std::size_t taggerID{0}; taggerIDfixed() ) { + calib_p0_ave_[taggerID]->updatePull(); + } + if ( not calib_p0_delta_[taggerID]->fixed() ) { + calib_p0_delta_[taggerID]->updatePull(); + } + if ( not calib_p1_ave_[taggerID]->fixed() ) { + calib_p1_ave_[taggerID]->updatePull(); + } + if ( not calib_p1_delta_[taggerID]->fixed() ) { + calib_p1_delta_[taggerID]->updatePull(); + } + } + } else { + for ( std::size_t taggerID{0}; taggerIDfixed() ) { + calib_p0_B0_[taggerID]->updatePull(); + } + if ( not calib_p0_B0bar_[taggerID]->fixed() ) { + calib_p0_B0bar_[taggerID]->updatePull(); + } + if ( not calib_p1_B0_[taggerID]->fixed() ) { + calib_p1_B0_[taggerID]->updatePull(); + } + if ( not calib_p1_B0bar_[taggerID]->fixed() ) { + calib_p1_B0bar_[taggerID]->updatePull(); + } + } + } +} diff --git a/src/LauTimeDepFitModel.cc b/src/LauTimeDepFitModel.cc index 32beb50..d5d48b7 100644 --- a/src/LauTimeDepFitModel.cc +++ b/src/LauTimeDepFitModel.cc @@ -1,3190 +1,3193 @@ /* Copyright 2006 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauTimeDepFitModel.cc \brief File containing implementation of LauTimeDepFitModel class. */ #include #include #include #include #include #include #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 "LauParamFixed.hh" #include "LauPrint.hh" #include "LauRandom.hh" #include "LauScfMap.hh" #include "LauTimeDepFitModel.hh" #include "LauFlavTag.hh" ClassImp(LauTimeDepFitModel) LauTimeDepFitModel::LauTimeDepFitModel(LauIsobarDynamics* modelB0bar, LauIsobarDynamics* modelB0, LauFlavTag* flavTag) : LauAbsFitModel(), sigModelB0bar_(modelB0bar), sigModelB0_(modelB0), kinematicsB0bar_(modelB0bar ? modelB0bar->getKinematics() : 0), kinematicsB0_(modelB0 ? modelB0->getKinematics() : 0), usingBkgnd_(kFALSE), flavTag_(flavTag), curEvtTrueTagFlv_(LauFlavTag::Flavour::Unknown), curEvtDecayFlv_(LauFlavTag::Flavour::Unknown), nSigComp_(0), nSigDPPar_(0), nDecayTimePar_(0), nExtraPdfPar_(0), nNormPar_(0), nCalibPar_(0), nTagEffPar_(0), nEffiPar_(0), nAsymPar_(0), coeffsB0bar_(0), coeffsB0_(0), coeffPars_(0), fitFracB0bar_(0), fitFracB0_(0), fitFracAsymm_(0), acp_(0), meanEffB0bar_("meanEffB0bar",0.0,0.0,1.0), meanEffB0_("meanEffB0",0.0,0.0,1.0), DPRateB0bar_("DPRateB0bar",0.0,0.0,100.0), DPRateB0_("DPRateB0",0.0,0.0,100.0), signalEvents_(0), signalAsym_(0), cpevVarName_(""), cpEigenValue_(CPEven), evtCPEigenVals_(0), deltaM_("deltaM",0.0), deltaGamma_("deltaGamma",0.0), tau_("tau",LauConstants::tauB0), phiMix_("phiMix", 2.0*LauConstants::beta, -LauConstants::threePi, LauConstants::threePi, kFALSE), sinPhiMix_("sinPhiMix", TMath::Sin(2.0*LauConstants::beta), -1.0, 1.0, kFALSE), cosPhiMix_("cosPhiMix", TMath::Cos(2.0*LauConstants::beta), -1.0, 1.0, kFALSE), useSinCos_(kFALSE), phiMixComplex_(TMath::Cos(-2.0*LauConstants::beta),TMath::Sin(-2.0*LauConstants::beta)), signalDecayTimePdf_(), BkgndTypes_(flavTag_->getBkgndTypes()), BkgndDecayTimePdfs_(), curEvtDecayTime_(0.0), curEvtDecayTimeErr_(0.0), sigExtraPdf_(), AProd_("AProd",0.0,-1.0,1.0,kTRUE), iterationsMax_(100000000), nGenLoop_(0), ASq_(0.0), aSqMaxVar_(0.0), aSqMaxSet_(1.25), storeGenAmpInfo_(kFALSE), signalTree_(), reuseSignal_(kFALSE), sigDPLike_(0.0), sigExtraLike_(0.0), sigTotalLike_(0.0) { // Set up ftag here? this->setBkgndClassNames(flavTag_->getBkgndNames()); const std::size_t nBkgnds { this->nBkgndClasses() }; if ( nBkgnds > 0 ){ usingBkgnd_ = kTRUE; for ( std::size_t iBkgnd{0}; iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass { this->bkgndClassName( iBkgnd ) }; AProdBkgnd_[iBkgnd] = new LauParameter("AProd_"+bkgndClass,0.0,-1.0,1.0,kTRUE); } } // Make sure that the integration scheme will be symmetrised sigModelB0bar_->forceSymmetriseIntegration(kTRUE); sigModelB0_->forceSymmetriseIntegration(kTRUE); } LauTimeDepFitModel::~LauTimeDepFitModel() { for ( LauAbsPdf* pdf : sigExtraPdf_ ) { delete pdf; } for(auto& data : bkgndTree_){ delete data; } } void LauTimeDepFitModel::setupBkgndVectors() { UInt_t nBkgnds { this->nBkgndClasses() }; AProdBkgnd_.resize( nBkgnds ); BkgndDPModelsB_.resize( nBkgnds ); BkgndDPModelsBbar_.resize( nBkgnds ); BkgndDecayTimePdfs_.resize( nBkgnds ); BkgndPdfs_.resize( nBkgnds ); bkgndEvents_.resize( nBkgnds ); bkgndAsym_.resize( nBkgnds ); bkgndTree_.resize( nBkgnds ); reuseBkgnd_.resize( nBkgnds ); bkgndDPLike_.resize( nBkgnds ); bkgndExtraLike_.resize( nBkgnds ); bkgndTotalLike_.resize( nBkgnds ); } void LauTimeDepFitModel::setNSigEvents(LauParameter* nSigEvents) { if ( nSigEvents == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( signalEvents_ != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl; return; } if ( signalAsym_ != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl; return; } signalEvents_ = nSigEvents; signalEvents_->name("signalEvents"); Double_t value = nSigEvents->value(); signalEvents_->range(-2.0*(TMath::Abs(value)+1.0),2.0*(TMath::Abs(value)+1.0)); signalAsym_ = new LauParameter("signalAsym",0.0,-1.0,1.0,kTRUE); } void LauTimeDepFitModel::setNSigEvents(LauParameter* nSigEvents, LauParameter* sigAsym) { if ( nSigEvents == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The event LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( sigAsym == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The asym LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( signalEvents_ != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl; return; } if ( signalAsym_ != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl; return; } signalEvents_ = nSigEvents; signalEvents_->name("signalEvents"); Double_t value = nSigEvents->value(); signalEvents_->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0)); signalAsym_ = sigAsym; signalAsym_->name("signalAsym"); signalAsym_->range(-1.0,1.0); } void LauTimeDepFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents) { if ( nBkgndEvents == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBgkndEvents : The background yield LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; gSystem->Exit(EXIT_FAILURE); } UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() ); if ( bkgndEvents_[bkgndID] != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl; return; } if ( bkgndAsym_[bkgndID] != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl; return; } nBkgndEvents->name( nBkgndEvents->name()+"Events" ); if ( nBkgndEvents->isLValue() ) { Double_t value = nBkgndEvents->value(); LauParameter* yield = dynamic_cast( nBkgndEvents ); yield->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0)); } bkgndEvents_[bkgndID] = nBkgndEvents; bkgndAsym_[bkgndID] = new LauParameter(nBkgndEvents->name()+"Asym",0.0,-1.0,1.0,kTRUE); } void LauTimeDepFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym) { if ( nBkgndEvents == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : The background yield LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( bkgndAsym == 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : The background asym LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; gSystem->Exit(EXIT_FAILURE); } UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() ); if ( bkgndEvents_[bkgndID] != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl; return; } if ( bkgndAsym_[bkgndID] != 0 ) { std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl; return; } nBkgndEvents->name( nBkgndEvents->name()+"Events" ); if ( nBkgndEvents->isLValue() ) { Double_t value = nBkgndEvents->value(); LauParameter* yield = dynamic_cast( nBkgndEvents ); yield->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0)); } bkgndEvents_[bkgndID] = nBkgndEvents; bkgndAsym->name( nBkgndEvents->name()+"Asym" ); if ( bkgndAsym->isLValue() ) { LauParameter* asym = dynamic_cast( bkgndAsym ); asym->range(-1.0, 1.0); } bkgndAsym_[bkgndID] = bkgndAsym; } void LauTimeDepFitModel::setSignalDtPdf(LauDecayTimePdf* pdf) { if (pdf==0) { std::cerr<<"ERROR in LauTimeDepFitModel::setSignalDtPdf : The PDF pointer is null, not adding it."<validBkgndClass( bkgndClass) ) { std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDtPdf : Invalid background class \"" << bkgndClass << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; return; } UInt_t bkgndID = this->bkgndClassID( bkgndClass ); BkgndDecayTimePdfs_[bkgndID] = pdf; usingBkgnd_ = kTRUE; } void LauTimeDepFitModel::setBkgndDPModels(const TString& bkgndClass, LauAbsBkgndDPModel* BModel, LauAbsBkgndDPModel* BbarModel) { if (BModel==nullptr) { std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDPModels : the model pointer is null for the particle model." << std::endl; return; } // check that this background name is valid if ( ! this->validBkgndClass( bkgndClass) ) { std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDPModels : Invalid background class \"" << bkgndClass << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; return; } UInt_t bkgndID = this->bkgndClassID( bkgndClass ); BkgndDPModelsB_[bkgndID] = BModel; if (BbarModel==nullptr) { std::cout << "INFO in LauTimeDepFitModel::setBkgndDPModels : the model pointer is null for the anti-particle model. Using only the particle model." << std::endl; BkgndDPModelsBbar_[bkgndID] = nullptr; } else { BkgndDPModelsBbar_[bkgndID] = BbarModel; } usingBkgnd_ = kTRUE; } void LauTimeDepFitModel::setSignalPdfs(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 (pdf==0) { std::cerr<<"ERROR in LauTimeDepFitModel::setSignalPdfs : The PDF pointer is null."<validBkgndClass( bkgndClass ) ) { std::cerr << "ERROR in LauTimeDepFitModel::setBkgndPdf : Invalid background class \"" << bkgndClass << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; return; } UInt_t bkgndID = this->bkgndClassID( bkgndClass ); BkgndPdfs_[bkgndID].push_back(pdf); usingBkgnd_ = kTRUE; } void LauTimeDepFitModel::setPhiMix(const Double_t phiMix, const Bool_t fixPhiMix, const Bool_t useSinCos) { phiMix_.value(phiMix); phiMix_.initValue(phiMix); phiMix_.genValue(phiMix); phiMix_.fixed(fixPhiMix); const Double_t sinPhiMix = TMath::Sin(phiMix); sinPhiMix_.value(sinPhiMix); sinPhiMix_.initValue(sinPhiMix); sinPhiMix_.genValue(sinPhiMix); sinPhiMix_.fixed(fixPhiMix); const Double_t cosPhiMix = TMath::Cos(phiMix); cosPhiMix_.value(cosPhiMix); cosPhiMix_.initValue(cosPhiMix); cosPhiMix_.genValue(cosPhiMix); cosPhiMix_.fixed(fixPhiMix); useSinCos_ = useSinCos; phiMixComplex_.setRealPart(cosPhiMix); phiMixComplex_.setImagPart(-1.0*sinPhiMix); } void LauTimeDepFitModel::initialise() { // From the initial parameter values calculate the coefficients // so they can be passed to the signal model this->updateCoeffs(); // Initialisation if (this->useDP() == kTRUE) { this->initialiseDPModels(); } // Flavour tagging //flavTag_->initialise(); // Decay-time PDFs signalDecayTimePdf_->initialise(); //Initialise for backgrounds if necessary for (auto& pdf : BkgndDecayTimePdfs_){ pdf->initialise(); } if (!this->useDP() && sigExtraPdf_.empty()) { std::cerr<<"ERROR in LauTimeDepFitModel::initialise : Signal model doesn't exist for any variable."<Exit(EXIT_FAILURE); } if (this->useDP() == kTRUE) { // Check that we have all the Dalitz-plot models if ((sigModelB0bar_ == 0) || (sigModelB0_ == 0)) { std::cerr<<"ERROR in LauTimeDepFitModel::initialise : the pointer to one (particle or anti-particle) of the signal DP models is null."<Exit(EXIT_FAILURE); } } // Next check that, if a given component is being used we've got the // right number of PDFs for all the variables involved // TODO - should probably check variable names and so on as well //UInt_t nsigpdfvars(0); //for ( LauPdfPList::const_iterator pdf_iter = sigExtraPdf_.begin(); pdf_iter != sigExtraPdf_.end(); ++pdf_iter ) { // std::vector varNames = (*pdf_iter)->varNames(); // for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { // if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { // ++nsigpdfvars; // } // } //} //if (usingBkgnd_) { // for (LauBkgndPdfsList::const_iterator bgclass_iter = BkgndPdfsB0_.begin(); bgclass_iter != BkgndPdfsB0_.end(); ++bgclass_iter) { // UInt_t nbkgndpdfvars(0); // const LauPdfPList& pdfList = (*bgclass_iter); // for ( LauPdfPList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter ) { // std::vector varNames = (*pdf_iter)->varNames(); // for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { // if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { // ++nbkgndpdfvars; // } // } // } // if (nbkgndpdfvars != nsigpdfvars) { // std::cerr << "ERROR in LauTimeDepFitModel::initialise : There are " << nsigpdfvars << " signal PDF variables but " << nbkgndpdfvars << " bkgnd PDF variables." << std::endl; // gSystem->Exit(EXIT_FAILURE); // } // } //} // Clear the vectors of parameter information so we can start from scratch this->clearFitParVectors(); // Set the fit parameters for signal and background models this->setSignalDPParameters(); // Set the fit parameters for the decay time models this->setDecayTimeParameters(); // Set the fit parameters for the extra PDFs this->setExtraPdfParameters(); // Set the initial bg and signal events this->setFitNEvents(); // Handle flavour-tagging calibration parameters this->setCalibParams(); // Add tagging efficiency parameters this->setTagEffParams(); //Asymmetry terms AProd and in setAsymmetries()? this->setAsymParams(); // Check that we have the expected number of fit variables const LauParameterPList& fitVars = this->fitPars(); if (fitVars.size() != (nSigDPPar_ + nDecayTimePar_ + nExtraPdfPar_ + nNormPar_ + nCalibPar_ + nTagEffPar_ + nEffiPar_ + nAsymPar_)) { std::cerr<<"ERROR in LauTimeDepFitModel::initialise : Number of fit parameters not of expected size."<Exit(EXIT_FAILURE); } if (sigModelB0_ == 0) { std::cerr<<"ERROR in LauTimeDepFitModel::initialiseDPModels : B0 signal DP model doesn't exist"<Exit(EXIT_FAILURE); } // Need to check that the number of components we have and that the dynamics has matches up const UInt_t nAmpB0bar = sigModelB0bar_->getnTotAmp(); const UInt_t nAmpB0 = sigModelB0_->getnTotAmp(); if ( nAmpB0bar != nAmpB0 ) { std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Unequal number of signal DP components in the particle and anti-particle models: " << nAmpB0bar << " != " << nAmpB0 << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( nAmpB0bar != nSigComp_ ) { std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Number of signal DP components in the model (" << nAmpB0bar << ") not equal to number of coefficients supplied (" << nSigComp_ << ")." << std::endl; gSystem->Exit(EXIT_FAILURE); } std::cout<<"INFO in LauTimeDepFitModel::initialiseDPModels : Initialising signal DP model"<initialise(coeffsB0bar_); sigModelB0_->initialise(coeffsB0_); fifjEffSum_.clear(); fifjEffSum_.resize(nSigComp_); for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) { fifjEffSum_[iAmp].resize(nSigComp_); } // calculate the integrals of the A*Abar terms this->calcInterferenceTermIntegrals(); this->calcInterferenceTermNorm(); // Add backgrounds if (usingBkgnd_ == kTRUE) { for (auto& model : BkgndDPModelsB_){ model->initialise(); } for (auto& model : BkgndDPModelsBbar_){ if (model != nullptr) { model->initialise(); } } } } void LauTimeDepFitModel::calcInterferenceTermIntegrals() { const std::vector& integralInfoListB0bar = sigModelB0bar_->getIntegralInfos(); const std::vector& integralInfoListB0 = sigModelB0_->getIntegralInfos(); // TODO should check (first time) that they match in terms of number of entries in the vectors and that each entry has the same number of points, ranges, weights etc. LauComplex A, Abar, fifjEffSumTerm; for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) { for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) { fifjEffSum_[iAmp][jAmp].zero(); } } const UInt_t nIntegralRegions = integralInfoListB0bar.size(); for ( UInt_t iRegion(0); iRegion < nIntegralRegions; ++iRegion ) { const LauDPPartialIntegralInfo* integralInfoB0bar = integralInfoListB0bar[iRegion]; const LauDPPartialIntegralInfo* integralInfoB0 = integralInfoListB0[iRegion]; const UInt_t nm13Points = integralInfoB0bar->getnm13Points(); const UInt_t nm23Points = integralInfoB0bar->getnm23Points(); for (UInt_t m13 = 0; m13 < nm13Points; ++m13) { for (UInt_t m23 = 0; m23 < nm23Points; ++m23) { const Double_t weight = integralInfoB0bar->getWeight(m13,m23); const Double_t eff = integralInfoB0bar->getEfficiency(m13,m23); const Double_t effWeight = eff*weight; for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) { A = integralInfoB0->getAmplitude(m13, m23, iAmp); for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) { Abar = integralInfoB0bar->getAmplitude(m13, m23, jAmp); fifjEffSumTerm = Abar*A.conj(); fifjEffSumTerm.rescale(effWeight); fifjEffSum_[iAmp][jAmp] += fifjEffSumTerm; } } } } } } void LauTimeDepFitModel::calcInterferenceTermNorm() { const std::vector& fNormB0bar = sigModelB0bar_->getFNorm(); const std::vector& fNormB0 = sigModelB0_->getFNorm(); LauComplex norm; for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) { for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) { LauComplex coeffTerm = coeffsB0bar_[jAmp]*coeffsB0_[iAmp].conj(); coeffTerm *= fifjEffSum_[iAmp][jAmp]; coeffTerm.rescale(fNormB0bar[jAmp] * fNormB0[iAmp]); norm += coeffTerm; } } norm *= phiMixComplex_; interTermReNorm_ = 2.0*norm.re(); interTermImNorm_ = 2.0*norm.im(); } void LauTimeDepFitModel::setAmpCoeffSet(LauAbsCoeffSet* coeffSet) { // Is there a component called compName in the signal models? TString compName = coeffSet->name(); TString conjName = sigModelB0bar_->getConjResName(compName); const LauDaughters* daughtersB0bar = sigModelB0bar_->getDaughters(); const LauDaughters* daughtersB0 = sigModelB0_->getDaughters(); const Bool_t conjugate = daughtersB0bar->isConjugate( daughtersB0 ); if ( ! sigModelB0bar_->hasResonance(compName) ) { if ( ! sigModelB0bar_->hasResonance(conjName) ) { std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0bar signal DP model doesn't contain component \""<name( compName ); } if ( conjugate ) { if ( ! sigModelB0_->hasResonance(conjName) ) { std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<hasResonance(compName) ) { std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<name() == compName) { std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : Have already set coefficients for \""<index(nSigComp_); coeffPars_.push_back(coeffSet); TString parName = coeffSet->baseName(); parName += "FitFracAsym"; fitFracAsymm_.push_back(LauParameter(parName, 0.0, -1.0, 1.0)); acp_.push_back(coeffSet->acp()); ++nSigComp_; std::cout<<"INFO in LauTimeDepFitModel::setAmpCoeffSet : Added coefficients for component \""<acp(); LauAsymmCalc asymmCalc(fitFracB0bar_[i][i].value(), fitFracB0_[i][i].value()); Double_t asym = asymmCalc.getAsymmetry(); fitFracAsymm_[i].value(asym); if (initValues) { fitFracAsymm_[i].genValue(asym); fitFracAsymm_[i].initValue(asym); } } } void LauTimeDepFitModel::setSignalDPParameters() { // Set the fit parameters for the signal model. nSigDPPar_ = 0; if ( ! this->useDP() ) { return; } std::cout << "INFO in LauTimeDepFitModel::setSignalDPParameters : Setting the initial fit parameters for the signal DP model." << std::endl; // Place isobar coefficient parameters in vector of fit variables for (UInt_t i = 0; i < nSigComp_; ++i) { LauParameterPList pars = coeffPars_[i]->getParameters(); nSigDPPar_ += this->addFitParameters( pars, kTRUE ); } // 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 LauParameterPList& sigDPParsB0bar = sigModelB0bar_->getFloatingParameters(); LauParameterPList& sigDPParsB0 = sigModelB0_->getFloatingParameters(); nSigDPPar_ += this->addResonanceParameters( sigDPParsB0bar ); nSigDPPar_ += this->addResonanceParameters( sigDPParsB0 ); } UInt_t LauTimeDepFitModel::addFitParameters(LauDecayTimePdf* decayTimePdf) { return this->addFitParameters( decayTimePdf->getParameters(), kTRUE ); } UInt_t LauTimeDepFitModel::addFitParameters(std::vector& decayTimePdfList) { UInt_t nParsAdded{0}; for ( auto decayTimePdf : decayTimePdfList ) { nParsAdded += this->addFitParameters( decayTimePdf ); } return nParsAdded; } void LauTimeDepFitModel::setDecayTimeParameters() { nDecayTimePar_ = 0; std::cout << "INFO in LauTimeDepFitModel::setDecayTimeParameters : Setting the initial fit parameters of the DecayTime Pdfs." << std::endl; // Loop over the Dt PDFs nDecayTimePar_ += this->addFitParameters( signalDecayTimePdf_ ); if (usingBkgnd_){ nDecayTimePar_ += this->addFitParameters(BkgndDecayTimePdfs_); } if (useSinCos_) { nDecayTimePar_ += this->addFitParameters( &sinPhiMix_ ); nDecayTimePar_ += this->addFitParameters( &cosPhiMix_ ); } else { nDecayTimePar_ += this->addFitParameters( &phiMix_ ); } } 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; std::cout << "INFO in LauTimeDepFitModel::setExtraPdfParameters : Setting the initial fit parameters of the extra Pdfs." << std::endl; nExtraPdfPar_ += this->addFitParameters(sigExtraPdf_); if (usingBkgnd_ == kTRUE) { for (auto& pdf : BkgndPdfs_){ nExtraPdfPar_ += this->addFitParameters(pdf); } } } void LauTimeDepFitModel::setFitNEvents() { nNormPar_ = 0; std::cout << "INFO in LauTimeDepFitModel::setFitNEvents : Setting the initial fit parameters of the signal and background yields." << std::endl; // Initialise the total number of events to be the sum of all the hypotheses Double_t nTotEvts = signalEvents_->value(); this->eventsPerExpt(TMath::FloorNint(nTotEvts)); // 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..."<addFitParameters( signalEvents_ ); } else { std::cout<<"INFO in LauTimeDepFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..."<useDP() == kFALSE) { nNormPar_ += this->addFitParameters( signalAsym_ ); } // TODO arguably should delegate this //LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac(); // tagging-category fractions for signal events //for (LauTagCatParamMap::iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) { // if (iter == signalTagCatFrac.begin()) { // continue; // } // LauParameter* par = &((*iter).second); // fitVars.push_back(par); // ++nNormPar_; //} // Backgrounds if (usingBkgnd_ == kTRUE) { nNormPar_ += this->addFitParameters( bkgndEvents_ ); nNormPar_ += this->addFitParameters( bkgndAsym_ ); } } void LauTimeDepFitModel::setAsymParams() { nAsymPar_ = 0; //Signal nAsymPar_ += this->addFitParameters( &AProd_ ); //Background(s) nAsymPar_ += this->addFitParameters( AProdBkgnd_ ); } void LauTimeDepFitModel::setTagEffParams() { nTagEffPar_ = 0; Bool_t useAltPars = flavTag_->getUseAveDelta(); std::cout << "INFO in LauTimeDepFitModel::setTagEffParams : Setting the initial fit parameters for flavour tagging efficiencies." << std::endl; if (useAltPars){ std::vector tageff_ave = flavTag_->getTagEffAve(); std::vector tageff_delta = flavTag_->getTagEffDelta(); nTagEffPar_ += this->addFitParameters( tageff_ave ); nTagEffPar_ += this->addFitParameters( tageff_delta ); } else { std::vector tageff_b0 = flavTag_->getTagEffB0(); std::vector tageff_b0bar = flavTag_->getTagEffB0bar(); nTagEffPar_ += this->addFitParameters( tageff_b0 ); nTagEffPar_ += this->addFitParameters( tageff_b0bar ); } if (usingBkgnd_){ if (useAltPars){ auto tageff_ave = flavTag_->getTagEffBkgndAve(); auto tageff_delta = flavTag_->getTagEffBkgndDelta(); for(auto& innerVec : tageff_ave){ nTagEffPar_ += this->addFitParameters( innerVec ); } for(auto& innerVec : tageff_delta){ nTagEffPar_ += this->addFitParameters( innerVec ); } } else { auto tageff_b0 = flavTag_->getTagEffBkgndB0(); auto tageff_b0bar = flavTag_->getTagEffBkgndB0bar(); for(auto& innerVec : tageff_b0){ nTagEffPar_ += this->addFitParameters( innerVec ); } for(auto& innerVec : tageff_b0bar){ nTagEffPar_ += this->addFitParameters( innerVec ); } } } } void LauTimeDepFitModel::setCalibParams() { nCalibPar_ = 0; Bool_t useAltPars = flavTag_->getUseAveDelta(); std::cout << "INFO in LauTimeDepFitModel::setCalibParams : Setting the initial fit parameters of the flavour tagging calibration parameters." << std::endl; if (useAltPars){ std::vector p0pars_ave = flavTag_->getCalibP0Ave(); std::vector p0pars_delta = flavTag_->getCalibP0Delta(); std::vector p1pars_ave = flavTag_->getCalibP1Ave(); std::vector p1pars_delta = flavTag_->getCalibP1Delta(); nCalibPar_ += this->addFitParameters( p0pars_ave ); nCalibPar_ += this->addFitParameters( p0pars_delta ); nCalibPar_ += this->addFitParameters( p1pars_ave ); nCalibPar_ += this->addFitParameters( p1pars_delta ); } else { std::vector p0pars_b0 = flavTag_->getCalibP0B0(); std::vector p0pars_b0bar = flavTag_->getCalibP0B0bar(); std::vector p1pars_b0 = flavTag_->getCalibP1B0(); std::vector p1pars_b0bar = flavTag_->getCalibP1B0bar(); nCalibPar_ += this->addFitParameters( p0pars_b0 ); nCalibPar_ += this->addFitParameters( p0pars_b0bar ); nCalibPar_ += this->addFitParameters( p1pars_b0 ); nCalibPar_ += this->addFitParameters( p1pars_b0bar ); } } void LauTimeDepFitModel::setExtraNtupleVars() { // Set-up other parameters derived from the fit results, e.g. fit fractions. if (this->useDP() != kTRUE) { return; } // First clear the vectors so we start from scratch this->clearExtraVarVectors(); LauParameterList& extraVars = this->extraPars(); // Add the B0 and B0bar fit fractions for each signal component fitFracB0bar_ = sigModelB0bar_->getFitFractions(); if (fitFracB0bar_.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); igetFitFractions(); if (fitFracB0_.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); icalcAsymmetries(kTRUE); // Add the Fit Fraction asymmetry for each signal component for (UInt_t i = 0; i < nSigComp_; i++) { extraVars.push_back(fitFracAsymm_[i]); } // Add the calculated CP asymmetry for each signal component for (UInt_t i = 0; i < nSigComp_; i++) { extraVars.push_back(acp_[i]); } // Now add in the DP efficiency values Double_t initMeanEffB0bar = sigModelB0bar_->getMeanEff().initValue(); meanEffB0bar_.value(initMeanEffB0bar); meanEffB0bar_.initValue(initMeanEffB0bar); meanEffB0bar_.genValue(initMeanEffB0bar); extraVars.push_back(meanEffB0bar_); Double_t initMeanEffB0 = sigModelB0_->getMeanEff().initValue(); meanEffB0_.value(initMeanEffB0); meanEffB0_.initValue(initMeanEffB0); meanEffB0_.genValue(initMeanEffB0); extraVars.push_back(meanEffB0_); // Also add in the DP rates Double_t initDPRateB0bar = sigModelB0bar_->getDPRate().initValue(); DPRateB0bar_.value(initDPRateB0bar); DPRateB0bar_.initValue(initDPRateB0bar); DPRateB0bar_.genValue(initDPRateB0bar); extraVars.push_back(DPRateB0bar_); Double_t initDPRateB0 = sigModelB0_->getDPRate().initValue(); DPRateB0_.value(initDPRateB0); DPRateB0_.initValue(initDPRateB0); DPRateB0_.genValue(initDPRateB0); extraVars.push_back(DPRateB0_); } void LauTimeDepFitModel::setAsymmetries(const Double_t AProd, const Bool_t AProdFix){ AProd_.value(AProd); AProd_.fixed(AProdFix); } void LauTimeDepFitModel::setBkgndAsymmetries(const TString& bkgndClass, const Double_t AProd, const Bool_t AProdFix){ // check that this background name is valid if ( ! this->validBkgndClass( bkgndClass) ) { std::cerr << "ERROR in LauTimeDepFitModel::setBkgndAsymmetries : Invalid background class \"" << bkgndClass << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; return; } UInt_t bkgndID = this->bkgndClassID( bkgndClass ); AProdBkgnd_[bkgndID]->value( AProd ); AProdBkgnd_[bkgndID]->genValue( AProd ); AProdBkgnd_[bkgndID]->initValue( AProd ); AProdBkgnd_[bkgndID]->fixed( AProdFix ); } 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 signalDecayTimePdf_->updatePulls(); // and for backgrounds if required if (usingBkgnd_){ for (auto& pdf : BkgndDecayTimePdfs_){ pdf->updatePulls(); } } + // Finalise the pulls on the flavour tagging parameters + flavTag_->updatePulls(); + if (useSinCos_) { if ( not sinPhiMix_.fixed() ) { sinPhiMix_.updatePull(); cosPhiMix_.updatePull(); } } else { this->checkMixingPhase(); } if (usingBkgnd_ == kTRUE) { for (auto& params : bkgndEvents_){ std::vector parameters = params->getPars(); for ( LauParameter* parameter : parameters ) { parameter->updatePull(); } } for (auto& params : bkgndAsym_){ std::vector parameters = params->getPars(); for ( LauParameter* parameter : parameters ) { parameter->updatePull(); } } } // Update the pulls on all the extra PDFs' parameters this->updateFitParameters(sigExtraPdf_); if (usingBkgnd_ == kTRUE) { for (auto& pdf : BkgndPdfs_){ this->updateFitParameters(pdf); } } // Fill the fit results to the ntuple // update the coefficients and then calculate the fit fractions and ACP's if (this->useDP() == kTRUE) { this->updateCoeffs(); sigModelB0bar_->updateCoeffs(coeffsB0bar_); sigModelB0bar_->calcExtraInfo(); sigModelB0_->updateCoeffs(coeffsB0_); sigModelB0_->calcExtraInfo(); LauParArray fitFracB0bar = sigModelB0bar_->getFitFractions(); if (fitFracB0bar.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } LauParArray fitFracB0 = sigModelB0_->getFitFractions(); if (fitFracB0.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); igetMeanEff().value()); meanEffB0_.value(sigModelB0_->getMeanEff().value()); DPRateB0bar_.value(sigModelB0bar_->getDPRate().value()); DPRateB0_.value(sigModelB0_->getDPRate().value()); this->calcAsymmetries(); // Then store the final fit parameters, and any extra parameters for // the signal model (e.g. fit fractions, FF asymmetries, ACPs, mean efficiency and DP rate) this->clearExtraVarVectors(); LauParameterList& extraVars = this->extraPars(); for (UInt_t i(0); iprintFitFractions(std::cout); this->printAsymmetries(std::cout); } const LauParameterPList& fitVars = this->fitPars(); const LauParameterList& extraVars = this->extraPars(); LauFitNtuple* ntuple = this->fitNtuple(); ntuple->storeParsAndErrors(fitVars, extraVars); // find out the correlation matrix for the parameters ntuple->storeCorrMatrix(this->iExpt(), this->fitStatus(), this->covarianceMatrix()); // Fill the data into ntuple ntuple->updateFitNtuple(); // Print out the partial fit fractions, phases and the // averaged efficiency, reweighted by the dynamics (and anything else) if (this->writeLatexTable()) { TString sigOutFileName(tablePrefixName); sigOutFileName += "_"; sigOutFileName += this->iExpt(); sigOutFileName += "Expt.tex"; this->writeOutTable(sigOutFileName); } } void LauTimeDepFitModel::printFitFractions(std::ostream& output) { // Print out Fit Fractions, total DP rate and mean efficiency // First for the B0bar events for (UInt_t i = 0; i < nSigComp_; i++) { const TString compName(coeffPars_[i]->name()); output<<"B0bar FitFraction for component "<useDP() == kTRUE) { // print the fit coefficients in one table coeffPars_.front()->printTableHeading(fout); for (UInt_t i = 0; i < nSigComp_; i++) { coeffPars_[i]->printTableRow(fout); } fout<<"\\hline"<name(); resName = resName.ReplaceAll("_", "\\_"); fout< =$ & $"; print.printFormat(fout, meanEffB0bar_.value()); fout << "$ & $"; print.printFormat(fout, meanEffB0_.value()); fout << "$ & & \\\\" << std::endl; if (useSinCos_) { fout << "$\\sinPhiMix =$ & $"; print.printFormat(fout, sinPhiMix_.value()); fout << " \\pm "; print.printFormat(fout, sinPhiMix_.error()); fout << "$ & & & & & & & \\\\" << std::endl; fout << "$\\cosPhiMix =$ & $"; print.printFormat(fout, cosPhiMix_.value()); fout << " \\pm "; print.printFormat(fout, cosPhiMix_.error()); fout << "$ & & & & & & & \\\\" << std::endl; } else { fout << "$\\phiMix =$ & $"; print.printFormat(fout, phiMix_.value()); fout << " \\pm "; print.printFormat(fout, phiMix_.error()); fout << "$ & & & & & & & \\\\" << std::endl; } fout << "\\hline \n\\end{tabular}" << std::endl; } if (!sigExtraPdf_.empty()) { fout<<"\\begin{tabular}{|l|c|}"<printFitParameters(sigExtraPdf_, fout); if (usingBkgnd_ == kTRUE && !BkgndPdfs_.empty()) { fout << "\\hline" << std::endl; fout << "\\Extra Background PDFs' Parameters: & \\\\" << std::endl; for (auto& pdf : BkgndPdfs_){ this->printFitParameters(pdf, fout); } } fout<<"\\hline \n\\end{tabular}"<updateSigEvents(); // Check whether we want to have randomised initial fit parameters for the signal model if (this->useRandomInitFitPars() == kTRUE) { this->randomiseInitFitPars(); } } void LauTimeDepFitModel::randomiseInitFitPars() { // Only randomise those parameters that are not fixed! std::cout<<"INFO in LauTimeDepFitModel::randomiseInitFitPars : Randomising the initial values of the coefficients of the DP components (and phiMix)..."<randomiseInitValues(); } phiMix_.randomiseValue(-LauConstants::pi, LauConstants::pi); if (useSinCos_) { sinPhiMix_.initValue(TMath::Sin(phiMix_.initValue())); cosPhiMix_.initValue(TMath::Cos(phiMix_.initValue())); } } LauTimeDepFitModel::LauGenInfo LauTimeDepFitModel::eventsToGenerate() { // Determine the number of events to generate for each hypothesis // If we're smearing then smear each one individually // NB this individual smearing has to be done individually per tagging category as well LauGenInfo nEvtsGen; // 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 ); } //TOD sigAysm doesn't do anything here? Double_t sigAsym(0.0); if (this->useDP() == kFALSE) { sigAsym = signalAsym_->genValue(); //TODO fill in here if we care } 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. if (this->doPoissonSmearing()) { nEvts = LauRandom::randomFun()->Poisson(signalEvents_->genValue()); } nEvtsGen["signal"] = std::make_pair( nEvts, evtWeight ); } std::cout<<"INFO in LauTimeDepFitModel::eventsToGenerate : Generating toy MC with:"<bkgndClassName(bkgndID)<<" background events = "<genValue()<eventsToGenerate(); Bool_t genOK(kTRUE); Int_t evtNum(0); const UInt_t nBkgnds = this->nBkgndClasses(); std::vector bkgndClassNames(nBkgnds); std::vector bkgndClassNamesGen(nBkgnds); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name( this->bkgndClassName(iBkgnd) ); bkgndClassNames[iBkgnd] = name; bkgndClassNamesGen[iBkgnd] = "gen"+name; } // Loop over the hypotheses and generate the appropriate number of // events for each one for (auto& hypo : nEvts){ // find the category of events (e.g. signal) const TString& evtCategory(hypo.first); // Type const TString& type(hypo.first); // Number of events Int_t nEvtsGen( hypo.second.first ); // get the event weight for this category const Double_t evtWeight( hypo.second.second ); for (Int_t iEvt(0); iEvtsetGenNtupleDoubleBranchValue( "evtWeight", evtWeight ); if (evtCategory == "signal") { this->setGenNtupleIntegerBranchValue("genSig",1); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], 0 ); } // All the generate*Event() methods have to fill in curEvtDecayTime_ and curEvtDecayTimeErr_ // In addition, generateSignalEvent has to decide on the tag and fill in curEvtTagFlv_ genOK = this->generateSignalEvent(); } else { this->setGenNtupleIntegerBranchValue("genSig",0); UInt_t bkgndID(0); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { Int_t gen(0); if ( bkgndClassNames[iBkgnd] == type ) { gen = 1; bkgndID = iBkgnd; } this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], gen ); } genOK = this->generateBkgndEvent(bkgndID); } if (!genOK) { // If there was a problem with the generation then break out and return. // The problem model will have adjusted itself so that all should be OK next time. break; } if (this->useDP() == kTRUE) { this->setDPDtBranchValues(); // store DP, decay time and tagging variables in the ntuple } // Store the event's tag and tagging category this->setGenNtupleIntegerBranchValue("cpEigenvalue", cpEigenValue_); const TString& trueTagVarName { flavTag_->getTrueTagVarName() }; if ( trueTagVarName != "" ) { this->setGenNtupleIntegerBranchValue(trueTagVarName, curEvtTrueTagFlv_); } if ( cpEigenValue_ == QFS ) { const TString& decayFlvVarName { flavTag_->getDecayFlvVarName() }; if ( decayFlvVarName == "" ) { std::cerr<<"ERROR in LauTimeDepFitModel::genExpt : Decay flavour variable not set for QFS decay, see LauFlavTag::setDecayFlvVarName()."<Exit(EXIT_FAILURE); } else { this->setGenNtupleIntegerBranchValue(decayFlvVarName, curEvtDecayFlv_); } } const std::vector& tagVarNames { flavTag_->getTagVarNames() }; const std::vector& mistagVarNames { flavTag_->getMistagVarNames() }; // Loop over the taggers - values set via generateSignalEvent const std::size_t nTaggers {flavTag_->getNTaggers()}; for (std::size_t i=0; isetGenNtupleIntegerBranchValue(tagVarNames[i], curEvtTagFlv_[i]); this->setGenNtupleDoubleBranchValue(mistagVarNames[i], curEvtMistag_[i]); } // Store the event number (within this experiment) // and then increment it this->setGenNtupleIntegerBranchValue("iEvtWithinExpt",evtNum); ++evtNum; // Write the values into the tree this->fillGenNtupleBranches(); // Print an occasional progress message if (iEvt%1000 == 0) {std::cout<<"INFO in LauTimeDepFitModel::genExpt : Generated event number "<useDP() && genOK) { sigModelB0bar_->checkToyMC(kTRUE); sigModelB0_->checkToyMC(kTRUE); std::cout<<"aSqMaxSet = "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } LauParArray fitFracB0 = sigModelB0_->getFitFractions(); if (fitFracB0.size() != nSigComp_) { std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); igetMeanEff().value()); meanEffB0_.value(sigModelB0_->getMeanEff().value()); DPRateB0bar_.value(sigModelB0bar_->getDPRate().value()); DPRateB0_.value(sigModelB0_->getDPRate().value()); } } // If we're reusing embedded events or if the generation is being // reset then clear the lists of used events if (reuseSignal_ || !genOK) { if (signalTree_) { signalTree_->clearUsedList(); } } for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { LauEmbeddedData* data = bkgndTree_[bkgndID]; if (reuseBkgnd_[bkgndID] || !genOK) { if (data) { data->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); if (this->useDP()) { if (signalTree_) { signalTree_->getEmbeddedEvent(kinematicsB0bar_); //curEvtTagFlv_ = TMath::Nint(signalTree_->getValue("tagFlv")); curEvtDecayTimeErr_ = signalTree_->getValue(signalDecayTimePdf_->varErrName()); curEvtDecayTime_ = signalTree_->getValue(signalDecayTimePdf_->varName()); if (signalTree_->haveBranch("mcMatch")) { Int_t match = TMath::Nint(signalTree_->getValue("mcMatch")); if (match) { this->setGenNtupleIntegerBranchValue("genTMSig",1); this->setGenNtupleIntegerBranchValue("genSCFSig",0); } else { this->setGenNtupleIntegerBranchValue("genTMSig",0); this->setGenNtupleIntegerBranchValue("genSCFSig",1); } } } else { nGenLoop_ = 0; // Now generate from the combined DP / decay-time PDF while (generatedEvent == kFALSE && nGenLoop_ < iterationsMax_) { curEvtTrueTagFlv_ = LauFlavTag::Flavour::Unknown; curEvtDecayFlv_ = LauFlavTag::Flavour::Unknown; // First choose the true tag, accounting for the production asymmetry // CONVENTION WARNING regarding meaning of sign of AProd Double_t random = LauRandom::randomFun()->Rndm(); if (random <= 0.5 * ( 1.0 - AProd_.unblindValue() ) ) { curEvtTrueTagFlv_ = LauFlavTag::Flavour::B; } else { curEvtTrueTagFlv_ = LauFlavTag::Flavour::Bbar; } // 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); // Generate decay time const Double_t tMin = signalDecayTimePdf_->minAbscissa(); const Double_t tMax = signalDecayTimePdf_->maxAbscissa(); curEvtDecayTime_ = LauRandom::randomFun()->Uniform(tMin,tMax); // Generate the decay time error (NB the kTRUE forces the generation of a new value) curEvtDecayTimeErr_ = signalDecayTimePdf_->generateError(kTRUE); // Calculate all the decay time info signalDecayTimePdf_->calcLikelihoodInfo(curEvtDecayTime_,curEvtDecayTimeErr_); // Retrieve the amplitudes and efficiency from the dynamics const LauComplex& Abar { sigModelB0bar_->getEvtDPAmp() }; const LauComplex& A { sigModelB0_->getEvtDPAmp() }; const Double_t ASq { A.abs2() }; const Double_t AbarSq { Abar.abs2() }; const Double_t dpEff { sigModelB0bar_->getEvtEff() }; // Also retrieve all the decay time terms const Double_t dtCos { signalDecayTimePdf_->getCosTerm() }; const Double_t dtSin { signalDecayTimePdf_->getSinTerm() }; const Double_t dtCosh { signalDecayTimePdf_->getCoshTerm() }; const Double_t dtSinh { signalDecayTimePdf_->getSinhTerm() }; // and the decay time acceptance const Double_t dtEff { signalDecayTimePdf_->getEffiTerm() }; if ( cpEigenValue_ == QFS) { // Calculate the total intensities for each flavour-specific final state const Double_t ATotSq { ( ASq * dtCosh + curEvtTrueTagFlv_ * ASq * dtCos ) * dpEff * dtEff }; const Double_t AbarTotSq { ( AbarSq * dtCosh - curEvtTrueTagFlv_ * AbarSq * dtCos ) * dpEff * dtEff }; const Double_t ASumSq { ATotSq + AbarTotSq }; // Finally we throw the dice to see whether this event should be generated (and, if so, which final state) const Double_t randNum = LauRandom::randomFun()->Rndm(); if (randNum <= ASumSq / aSqMaxSet_ ) { generatedEvent = kTRUE; nGenLoop_ = 0; if (ASumSq > aSqMaxVar_) {aSqMaxVar_ = ASumSq;} if ( randNum <= ATotSq / aSqMaxSet_ ) { curEvtDecayFlv_ = LauFlavTag::Flavour::B; } else { curEvtDecayFlv_ = LauFlavTag::Flavour::Bbar; } // Generate the flavour tagging information from the true tag // (we do this after accepting the event to save time) flavTag_->generateEventInfo( curEvtTrueTagFlv_, curEvtDecayTime_ ); curEvtTagFlv_ = flavTag_->getCurEvtTagFlv(); curEvtMistag_ = flavTag_->getCurEvtMistag(); } else { nGenLoop_++; } } else { // Calculate the DP terms const Double_t aSqSum { ASq + AbarSq }; const Double_t aSqDif { ASq - AbarSq }; const LauComplex inter { Abar * A.conj() * phiMixComplex_ }; const Double_t interTermIm { ( cpEigenValue_ == CPEven ) ? 2.0 * inter.im() : -2.0 * inter.im() }; const Double_t interTermRe { ( cpEigenValue_ == CPEven ) ? 2.0 * inter.re() : -2.0 * inter.re() }; // Combine DP and decay-time info for all terms const Double_t coshTerm { aSqSum * dtCosh }; const Double_t sinhTerm { interTermRe * dtSinh }; const Double_t cosTerm { aSqDif * dtCos }; const Double_t sinTerm { interTermIm * dtSin }; // Sum to obtain the total and multiply by the efficiency // Multiplying the cos and sin terms by the true flavour at production const Double_t ATotSq { ( coshTerm + sinhTerm + curEvtTrueTagFlv_ * ( cosTerm - sinTerm ) ) * dpEff * dtEff }; //Finally we throw the dice to see whether this event should be generated const Double_t randNum = LauRandom::randomFun()->Rndm(); if (randNum <= ATotSq/aSqMaxSet_ ) { generatedEvent = kTRUE; nGenLoop_ = 0; if (ATotSq > aSqMaxVar_) {aSqMaxVar_ = ATotSq;} // Generate the flavour tagging information from the true tag // (we do this after accepting the event to save time) flavTag_->generateEventInfo( curEvtTrueTagFlv_, curEvtDecayTime_ ); curEvtTagFlv_ = flavTag_->getCurEvtTagFlv(); curEvtMistag_ = flavTag_->getCurEvtMistag(); } else { nGenLoop_++; } } } // end of while !generatedEvent loop } // end of if (signalTree_) else control } else { if ( signalTree_ ) { signalTree_->getEmbeddedEvent(0); //curEvtTagFlv_ = TMath::Nint(signalTree_->getValue("tagFlv")); curEvtDecayTimeErr_ = signalTree_->getValue(signalDecayTimePdf_->varErrName()); curEvtDecayTime_ = signalTree_->getValue(signalDecayTimePdf_->varName()); } } // Check whether we have generated the toy MC OK. if (nGenLoop_ >= iterationsMax_) { aSqMaxSet_ = 1.01 * aSqMaxVar_; genOK = kFALSE; std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Hit max iterations: setting aSqMaxSet_ to "< aSqMaxSet_) { aSqMaxSet_ = 1.01 * aSqMaxVar_; genOK = kFALSE; std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Found a larger ASq value: setting aSqMaxSet_ to "<updateKinematics(kinematicsB0bar_->getm13Sq(), kinematicsB0bar_->getm23Sq() ); this->generateExtraPdfValues(sigExtraPdf_, signalTree_); } // Check for problems with the embedding if (signalTree_ && (signalTree_->nEvents() == signalTree_->nUsedEvents())) { std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events."<clearUsedList(); } return genOK; } Bool_t LauTimeDepFitModel::generateBkgndEvent(UInt_t bkgndID) { // Generate Bkgnd event Bool_t genOK{kTRUE}; //Check necessary ingredients are in place //TODO these checks should be part of a general sanity check during the initialisation phase if (BkgndDPModelsB_[bkgndID] == nullptr){ std::cerr << "ERROR in LauTimeDepFitModel::generateBkgndEvent : Dalitz plot model is missing" << std::endl; gSystem->Exit(EXIT_FAILURE); } if (BkgndDecayTimePdfs_[bkgndID] == nullptr){ std::cerr << "ERROR in LauTimeDepFitModel::generateBkgndEvent : Decay time model is missing" << std::endl; gSystem->Exit(EXIT_FAILURE); } //TODO restore the ability to embed events from an external source //LauAbsBkgndDPModel* model(0); //LauEmbeddedData* embeddedData(0); //LauPdfPList* extraPdfs(0); //LauKinematics* kinematics(0); //model = BkgndDPModels_[bkgndID]; //if (this->enableEmbedding()) { // // find the right embedded data for the current tagging category // LauTagCatEmbDataMap::const_iterator emb_iter = bkgndTree_[bkgndID].find(curEvtTagCat_); // embeddedData = (emb_iter != bkgndTree_[bkgndID].end()) ? emb_iter->second : 0; //} //extraPdfs = &BkgndPdfs_[bkgndID]; //kinematics = kinematicsB0bar_; //if (this->useDP()) { // if (embeddedData) { // embeddedData->getEmbeddedEvent(kinematics); // } else { // if (model == 0) { // const TString& bkgndClass = this->bkgndClassName(bkgndID); // std::cerr << "ERROR in LauCPFitModel::generateBkgndEvent : Can't find the DP model for background class \"" << bkgndClass << "\"." << std::endl; // gSystem->Exit(EXIT_FAILURE); // } // genOK = model->generate(); // } //} else { // if (embeddedData) { // embeddedData->getEmbeddedEvent(0); // } //} //if (genOK) { // this->generateExtraPdfValues(extraPdfs, embeddedData); //} //// Check for problems with the embedding //if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) { // const TString& bkgndClass = this->bkgndClassName(bkgndID); // std::cerr << "WARNING in LauCPFitModel::generateBkgndEvent : Source of embedded " << bkgndClass << " events used up, clearing the list of used events." << std::endl; // embeddedData->clearUsedList(); //} // LauKinematics* kinematics{nullptr}; switch ( BkgndTypes_[bkgndID] ) { case LauFlavTag::BkgndType::Combinatorial: { // First choose the true tag, accounting for the production asymmetry // CONVENTION WARNING regarding meaning of sign of AProd // NB the true tag doesn't really mean anything for combinatorial background Double_t random = LauRandom::randomFun()->Rndm(); if ( random <= 0.5 * ( 1.0 - AProdBkgnd_[bkgndID]->unblindValue() ) ) { curEvtTrueTagFlv_ = LauFlavTag::Flavour::B; } else { curEvtTrueTagFlv_ = LauFlavTag::Flavour::Bbar; } kinematics = kinematicsB0_; if ( cpEigenValue_ == CPEigenvalue::QFS ) { if ( BkgndDPModelsBbar_[bkgndID] != nullptr ) { // generate the true decay flavour and the corresponding DP position // (the supply of two DP models indicates a possible asymmetry) const Double_t rateB { BkgndDPModelsB_[bkgndID]->getPdfNorm() }; const Double_t rateBbar { BkgndDPModelsBbar_[bkgndID]->getPdfNorm() }; const Double_t ADet { ( rateBbar - rateB ) / ( rateBbar + rateB ) }; random = LauRandom::randomFun()->Rndm(); if ( random <= 0.5 * ( 1.0 - ADet ) ) { curEvtDecayFlv_ = LauFlavTag::Flavour::B; BkgndDPModelsB_[bkgndID]->generate(); kinematics = kinematicsB0_; } else { curEvtDecayFlv_ = LauFlavTag::Flavour::Bbar; BkgndDPModelsBbar_[bkgndID]->generate(); kinematics = kinematicsB0bar_; } } else { // generate the true decay flavour // (the supply of only a single model indicates no asymmetry) random = LauRandom::randomFun()->Rndm(); if ( random <= 0.5 ) { curEvtDecayFlv_ = LauFlavTag::Flavour::B; } else { curEvtDecayFlv_ = LauFlavTag::Flavour::Bbar; } // generate the DP position BkgndDPModelsB_[bkgndID]->generate(); } } else { // mark that the decay flavour is unknown curEvtDecayFlv_ = LauFlavTag::Flavour::Unknown; // generate the DP position BkgndDPModelsB_[bkgndID]->generate(); } // generate decay time and its error curEvtDecayTimeErr_ = BkgndDecayTimePdfs_[bkgndID]->generateError(kTRUE); curEvtDecayTime_ = BkgndDecayTimePdfs_[bkgndID]->generate( kinematics ); // generate the flavour tagging information from the true tag and decay flavour // (we do this after accepting the event to save time) flavTag_->generateBkgndEventInfo( bkgndID, curEvtTrueTagFlv_, curEvtDecayFlv_, curEvtDecayTime_ ); curEvtTagFlv_ = flavTag_->getCurEvtTagFlv(); curEvtMistag_ = flavTag_->getCurEvtMistag(); break; } case LauFlavTag::BkgndType::FlavourSpecific: { const LauDecayTime::FuncType dtType { BkgndDecayTimePdfs_[bkgndID]->getFuncType() }; if ( dtType == LauDecayTime::FuncType::ExpTrig or dtType == LauDecayTime::FuncType::ExpHypTrig ) { nGenLoop_ = 0; Bool_t generatedEvent{kFALSE}; kinematics = kinematicsB0bar_; do { curEvtTrueTagFlv_ = LauFlavTag::Flavour::Unknown; curEvtDecayFlv_ = LauFlavTag::Flavour::Unknown; // First choose the true tag, accounting for the production asymmetry // CONVENTION WARNING regarding meaning of sign of AProd Double_t random = LauRandom::randomFun()->Rndm(); if (random <= 0.5 * ( 1.0 - AProdBkgnd_[bkgndID]->unblindValue() ) ) { curEvtTrueTagFlv_ = LauFlavTag::Flavour::B; } else { curEvtTrueTagFlv_ = LauFlavTag::Flavour::Bbar; } // Generate the DP position Double_t m13Sq{0.0}, m23Sq{0.0}; kinematics->genFlatPhaseSpace(m13Sq, m23Sq); // Next, calculate the total A^2 and Abar^2 for the given DP position BkgndDPModelsB_[bkgndID]->calcLikelihoodInfo(m13Sq, m23Sq); BkgndDPModelsBbar_[bkgndID]->calcLikelihoodInfo(m13Sq, m23Sq); // Generate decay time const Double_t tMin = BkgndDecayTimePdfs_[bkgndID]->minAbscissa(); const Double_t tMax = BkgndDecayTimePdfs_[bkgndID]->maxAbscissa(); curEvtDecayTime_ = LauRandom::randomFun()->Uniform(tMin,tMax); // Generate the decay time error (NB the kTRUE forces the generation of a new value) curEvtDecayTimeErr_ = BkgndDecayTimePdfs_[bkgndID]->generateError(kTRUE); // Calculate all the decay time info BkgndDecayTimePdfs_[bkgndID]->calcLikelihoodInfo(curEvtDecayTime_,curEvtDecayTimeErr_); // Retrieve the DP intensities const Double_t ASq { BkgndDPModelsB_[bkgndID]->getUnNormValue() }; const Double_t AbarSq { BkgndDPModelsBbar_[bkgndID]->getUnNormValue() }; // Also retrieve all the decay time terms const Double_t dtCos { BkgndDecayTimePdfs_[bkgndID]->getCosTerm() }; const Double_t dtCosh { BkgndDecayTimePdfs_[bkgndID]->getCoshTerm() }; // and the decay time acceptance const Double_t dtEff { BkgndDecayTimePdfs_[bkgndID]->getEffiTerm() }; // Calculate the total intensities for each flavour-specific final state const Double_t ATotSq { ( ASq * dtCosh + curEvtTrueTagFlv_ * ASq * dtCos ) * dtEff }; const Double_t AbarTotSq { ( AbarSq * dtCosh - curEvtTrueTagFlv_ * AbarSq * dtCos ) * dtEff }; const Double_t ASumSq { ATotSq + AbarTotSq }; // TODO - check if this really is the max possible const Double_t ASumSqMax { 2.0 * ( BkgndDPModelsB_[bkgndID]->getMaxHeight() + BkgndDPModelsBbar_[bkgndID]->getMaxHeight() ) }; // Finally we throw the dice to see whether this event should be generated (and, if so, which final state) const Double_t randNum = LauRandom::randomFun()->Rndm(); if (randNum <= ASumSq / ASumSqMax ) { generatedEvent = kTRUE; nGenLoop_ = 0; if (ASumSq > ASumSqMax) { std::cerr << "WARNING in LauTimeDepFitModel::generateBkgndEvent : ASumSq > ASumSqMax" << std::endl; } if ( randNum <= ATotSq / ASumSqMax ) { curEvtDecayFlv_ = LauFlavTag::Flavour::B; } else { curEvtDecayFlv_ = LauFlavTag::Flavour::Bbar; } // Generate the flavour tagging information from the true tag and decay flavour // (we do this after accepting the event to save time) flavTag_->generateBkgndEventInfo( bkgndID, curEvtTrueTagFlv_, curEvtDecayFlv_, curEvtDecayTime_ ); curEvtTagFlv_ = flavTag_->getCurEvtTagFlv(); curEvtMistag_ = flavTag_->getCurEvtMistag(); } else { nGenLoop_++; } } while (generatedEvent == kFALSE && nGenLoop_ < iterationsMax_); } else { // Hist, Delta, Exp, DeltaExp decay-time types // Since there are no oscillations for these decay-time types, // the true decay flavour must be equal to the true tag flavour // First choose the true tag and decay flavour, accounting for both the production and detection asymmetries // CONVENTION WARNING regarding meaning of sign of AProd and ADet const Double_t AProd { AProdBkgnd_[bkgndID]->unblindValue() }; const Double_t rateB { BkgndDPModelsB_[bkgndID]->getPdfNorm() }; const Double_t rateBbar { BkgndDPModelsBbar_[bkgndID]->getPdfNorm() }; const Double_t ADet { ( rateBbar - rateB ) / ( rateBbar + rateB ) }; const Double_t random = LauRandom::randomFun()->Rndm(); // TODO - is this the correct way to combine the production and detection asymmetries? if ( random <= 0.5 * ( 1.0 - AProd ) * ( 1.0 - ADet ) ) { curEvtDecayFlv_ = curEvtTrueTagFlv_ = LauFlavTag::Flavour::B; } else { curEvtDecayFlv_ = curEvtTrueTagFlv_ = LauFlavTag::Flavour::Bbar; } // generate the DP position if ( curEvtDecayFlv_ == LauFlavTag::Flavour::B ) { BkgndDPModelsB_[bkgndID]->generate(); kinematics = kinematicsB0_; } else { BkgndDPModelsBbar_[bkgndID]->generate(); kinematics = kinematicsB0bar_; } // generate decay time and its error curEvtDecayTimeErr_ = BkgndDecayTimePdfs_[bkgndID]->generateError(kTRUE); curEvtDecayTime_ = BkgndDecayTimePdfs_[bkgndID]->generate( kinematics ); // generate the flavour tagging information from the true tag and decay flavour // (we do this after accepting the event to save time) flavTag_->generateBkgndEventInfo( bkgndID, curEvtTrueTagFlv_, curEvtDecayFlv_, curEvtDecayTime_ ); curEvtTagFlv_ = flavTag_->getCurEvtTagFlv(); curEvtMistag_ = flavTag_->getCurEvtMistag(); } break; } case LauFlavTag::BkgndType::SelfConjugate: // TODO break; case LauFlavTag::BkgndType::NonSelfConjugate: // TODO break; } if ( genOK ) { // Make sure both kinematics objects are up-to-date kinematicsB0_->updateKinematics(kinematics->getm13Sq(), kinematics->getm23Sq() ); kinematicsB0bar_->updateKinematics(kinematics->getm13Sq(), kinematics->getm23Sq() ); this->generateExtraPdfValues(BkgndPdfs_[bkgndID], bkgndTree_[bkgndID]); } return genOK; } void LauTimeDepFitModel::setupGenNtupleBranches() { // Setup the required ntuple branches this->addGenNtupleDoubleBranch("evtWeight"); this->addGenNtupleIntegerBranch("genSig"); this->addGenNtupleIntegerBranch("cpEigenvalue"); const TString& trueTagVarName { flavTag_->getTrueTagVarName() }; if ( trueTagVarName != "" ) { this->addGenNtupleIntegerBranch(trueTagVarName); } if ( cpEigenValue_ == QFS ) { const TString& decayFlvVarName { flavTag_->getDecayFlvVarName() }; if ( decayFlvVarName == "" ) { std::cerr<<"ERROR in LauTimeDepFitModel::setupGenNtupleBranches : Decay flavour variable not set for QFS decay, see LauFlavTag::setDecayFlvVarName()."<Exit(EXIT_FAILURE); } else { this->addGenNtupleIntegerBranch(decayFlvVarName); } } const std::vector& tagVarNames { flavTag_->getTagVarNames() }; const std::vector& mistagVarNames { flavTag_->getMistagVarNames() }; const std::size_t nTaggers {flavTag_->getNTaggers()}; for (std::size_t taggerID{0}; taggerIDaddGenNtupleIntegerBranch(tagVarNames[taggerID]); this->addGenNtupleDoubleBranch(mistagVarNames[taggerID]); } if (this->useDP() == kTRUE) { // Let's add the decay time variables. this->addGenNtupleDoubleBranch(signalDecayTimePdf_->varName()); if ( signalDecayTimePdf_->varErrName() != "" ) { this->addGenNtupleDoubleBranch(signalDecayTimePdf_->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 for ( const LauAbsPdf* pdf : sigExtraPdf_ ) { const std::vector varNames{ pdf->varNames() }; for ( const TString& varName : varNames ) { if ( varName != "m13Sq" && varName != "m23Sq" ) { this->addGenNtupleDoubleBranch( varName ); } } } } void LauTimeDepFitModel::setDPDtBranchValues() { // Store the decay time variables. this->setGenNtupleDoubleBranchValue(signalDecayTimePdf_->varName(),curEvtDecayTime_); if ( signalDecayTimePdf_->varErrName() != "" ) { this->setGenNtupleDoubleBranchValue(signalDecayTimePdf_->varErrName(),curEvtDecayTimeErr_); } // CONVENTION WARNING // TODO check - for now use B0 for any tags //LauKinematics* kinematics(0); //if (curEvtTagFlv_[position]<0) { LauKinematics* 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(LauPdfPList& extraPdfs, LauEmbeddedData* embeddedData) { // CONVENTION WARNING LauKinematics* kinematics = kinematicsB0_; //LauKinematics* kinematics(0); //if (curEvtTagFlv_<0) { // kinematics = kinematicsB0_; //} else { // kinematics = kinematicsB0bar_; //} // Generate from the extra PDFs for (auto& pdf : extraPdfs){ LauFitData genValues; if (embeddedData) { genValues = embeddedData->getValues( pdf->varNames() ); } else { genValues = pdf->generate(kinematics); } for (auto& var : genValues){ TString varName = var.first; if ( varName != "m13Sq" && varName != "m23Sq" ) { Double_t value = var.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->calcInterferenceTermNorm(); } // Update the decay time normalisation if ( signalDecayTimePdf_ ) { signalDecayTimePdf_->propagateParUpdates(); } // TODO // - maybe also need to add an update of the background decay time PDFs here // Update the signal events from the background numbers if not doing an extended fit // And update the tagging category fractions this->updateSigEvents(); } void LauTimeDepFitModel::updateSigEvents() { // The background parameters will have been set from Minuit. // We need to update the signal events using these. if (!this->doEMLFit()) { Double_t nTotEvts = this->eventsPerExpt(); Double_t signalEvents = nTotEvts; signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts); for (auto& nBkgndEvents : bkgndEvents_){ if ( nBkgndEvents->isLValue() ) { LauParameter* yield = dynamic_cast( nBkgndEvents ); yield->range(-2.0*nTotEvts,2.0*nTotEvts); } } // Subtract background events (if any) from signal. if (usingBkgnd_ == kTRUE) { for (auto& nBkgndEvents : bkgndEvents_){ signalEvents -= nBkgndEvents->value(); } } if ( ! signalEvents_->fixed() ) { signalEvents_->value(signalEvents); } } } void LauTimeDepFitModel::cacheInputFitVars() { // Fill the internal data trees of the signal and background models. // Note that we store the events of both charges in both the // negative and the positive models. It's only later, at the stage // when the likelihood is being calculated, that we separate them. LauFitDataTree* inputFitData = this->fitData(); evtCPEigenVals_.clear(); const Bool_t hasCPEV = ( (cpevVarName_ != "") && inputFitData->haveBranch( cpevVarName_ ) ); UInt_t nEvents = inputFitData->nEvents(); evtCPEigenVals_.reserve( nEvents ); LauFitData::const_iterator fitdata_iter; for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) { const LauFitData& dataValues = inputFitData->getData(iEvt); // if the CP-eigenvalue is in the data use those, otherwise use the default if ( hasCPEV ) { fitdata_iter = dataValues.find( cpevVarName_ ); const Int_t cpEV = static_cast( fitdata_iter->second ); if ( cpEV == 1 ) { cpEigenValue_ = CPEven; } else if ( cpEV == -1 ) { cpEigenValue_ = CPOdd; } else if ( cpEV == 0 ) { cpEigenValue_ = QFS; } else { std::cerr<<"WARNING in LauTimeDepFitModel::cacheInputFitVars : Unknown value: "<useDP() == kTRUE) { // DecayTime and SigmaDecayTime signalDecayTimePdf_->cacheInfo(*inputFitData); // cache all the backgrounds too for(auto& bg : BkgndDecayTimePdfs_) {bg->cacheInfo(*inputFitData);} } // Flavour tagging information flavTag_->cacheInputFitVars(inputFitData,signalDecayTimePdf_->varName()); // ...and then the extra PDFs if (not sigExtraPdf_.empty()){ this->cacheInfo(sigExtraPdf_, *inputFitData); } if(usingBkgnd_ == kTRUE){ for (auto& pdf : BkgndPdfs_){ this->cacheInfo(pdf, *inputFitData); } } if (this->useDP() == kTRUE) { sigModelB0bar_->fillDataTree(*inputFitData); sigModelB0_->fillDataTree(*inputFitData); if (usingBkgnd_ == kTRUE) { for (auto& model : BkgndDPModelsB_){ model->fillDataTree(*inputFitData); } for (auto& model : BkgndDPModelsBbar_){ if (model != nullptr) { model->fillDataTree(*inputFitData); } } } } } Double_t LauTimeDepFitModel::getTotEvtLikelihood(const UInt_t 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 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_ * sigExtraLike_; Double_t signalEvents = signalEvents_->unblindValue(); // TODO - consider what to do here - do we even want the option not to use the DP in this model? //if ( not this->useDP() ) { //signalEvents *= 0.5 * (1.0 + curEvtTagFlv_ * signalAsym_->unblindValue()); //} // Construct the total event likelihood Double_t likelihood { sigLike * signalEvents }; if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { const Double_t bkgndEvents { bkgndEvents_[bkgndID]->unblindValue() }; likelihood += bkgndEvents*bkgndDPLike_[bkgndID]*bkgndExtraLike_[bkgndID]; } } return likelihood; } Double_t LauTimeDepFitModel::getEventSum() const { Double_t eventSum(0.0); eventSum += signalEvents_->unblindValue(); if (usingBkgnd_) { for ( const auto& yieldPar : bkgndEvents_ ) { eventSum += yieldPar->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. if ( ! this->useDP() ) { // There's always going to be a term in the likelihood for the // signal, so we'd better not zero it. sigDPLike_ = 1.0; const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { if (usingBkgnd_ == kTRUE) { bkgndDPLike_[bkgndID] = 1.0; } else { bkgndDPLike_[bkgndID] = 0.0; } } return; } // Calculate event quantities // Get the DP dynamics, decay time, and flavour tagging to calculate // everything required for the likelihood calculation sigModelB0bar_->calcLikelihoodInfo(iEvt); sigModelB0_->calcLikelihoodInfo(iEvt); signalDecayTimePdf_->calcLikelihoodInfo(static_cast(iEvt)); flavTag_->updateEventInfo(iEvt); // Retrieve the amplitudes and efficiency from the dynamics LauComplex Abar { sigModelB0bar_->getEvtDPAmp() }; LauComplex A { sigModelB0_->getEvtDPAmp() }; const Double_t dpEff { sigModelB0bar_->getEvtEff() }; // If this is a QFS decay, one of the DP amplitudes needs to be zeroed curEvtDecayFlv_ = LauFlavTag::Flavour::Unknown; if (cpEigenValue_ == QFS){ curEvtDecayFlv_ = flavTag_->getCurEvtDecayFlv(); if ( curEvtDecayFlv_ == LauFlavTag::Flavour::B ) { Abar.zero(); } else if ( curEvtDecayFlv_ == LauFlavTag::Flavour::Bbar ) { A.zero(); } else { std::cerr<<"ERROR in LauTimeDepFitModel::getEvtDPDtLikelihood : Decay flavour must be known for QFS decays."<Exit(EXIT_FAILURE); } } // Next calculate the DP terms const Double_t aSqSum { A.abs2() + Abar.abs2() }; const Double_t aSqDif { A.abs2() - Abar.abs2() }; Double_t interTermRe { 0.0 }; Double_t interTermIm { 0.0 }; if ( cpEigenValue_ != QFS ) { const LauComplex inter { Abar * A.conj() * phiMixComplex_ }; if ( cpEigenValue_ == CPEven ) { interTermIm = 2.0 * inter.im(); interTermRe = 2.0 * inter.re(); } else { interTermIm = -2.0 * inter.im(); interTermRe = -2.0 * inter.re(); } } // First get all the decay time terms // TODO Backgrounds // Get the decay time acceptance const Double_t dtEff { signalDecayTimePdf_->getEffiTerm() }; // Get all the decay time terms const Double_t dtCos { signalDecayTimePdf_->getCosTerm() }; const Double_t dtSin { signalDecayTimePdf_->getSinTerm() }; const Double_t dtCosh { signalDecayTimePdf_->getCoshTerm() }; const Double_t dtSinh { signalDecayTimePdf_->getSinhTerm() }; // Get the decay time error term const Double_t dtErrLike { signalDecayTimePdf_->getErrTerm() }; // Get flavour tagging terms Double_t omega{1.0}; Double_t omegabar{1.0}; const std::size_t nTaggers { flavTag_->getNTaggers() }; for (std::size_t taggerID{0}; taggerIDgetCapitalOmega(taggerID, LauFlavTag::Flavour::B); omegabar *= flavTag_->getCapitalOmega(taggerID, LauFlavTag::Flavour::Bbar); } const Double_t prodAsym { AProd_.unblindValue() }; const Double_t ftOmegaHyp { ((1.0 - prodAsym)*omega + (1.0 + prodAsym)*omegabar) }; const Double_t ftOmegaTrig { ((1.0 - prodAsym)*omega - (1.0 + prodAsym)*omegabar) }; const Double_t coshTerm { ftOmegaHyp * dtCosh * aSqSum }; const Double_t sinhTerm { ftOmegaHyp * dtSinh * interTermRe }; const Double_t cosTerm { ftOmegaTrig * dtCos * aSqDif }; const Double_t sinTerm { ftOmegaTrig * dtSin * interTermIm }; // Combine all terms to get the total amplitude squared const Double_t ASq { coshTerm + sinhTerm + cosTerm - sinTerm }; // Calculate the DP and time normalisation const Double_t normASqSum { sigModelB0_->getDPNorm() + sigModelB0bar_->getDPNorm() }; const Double_t normASqDiff { sigModelB0_->getDPNorm() - sigModelB0bar_->getDPNorm() }; Double_t normInterTermRe { 0.0 }; Double_t normInterTermIm { 0.0 }; if ( cpEigenValue_ != QFS ) { // TODO - double check this sign flipping here (it's presumably right but...) normInterTermRe = ( cpEigenValue_ == CPOdd ) ? -1.0 * interTermReNorm_ : interTermReNorm_; normInterTermIm = ( cpEigenValue_ == CPOdd ) ? -1.0 * interTermImNorm_ : interTermImNorm_; } const Double_t normCoshTerm { signalDecayTimePdf_->getNormTermCosh() }; const Double_t normSinhTerm { signalDecayTimePdf_->getNormTermSinh() }; const Double_t normCosTerm { signalDecayTimePdf_->getNormTermCos() }; const Double_t normSinTerm { signalDecayTimePdf_->getNormTermSin() }; const Double_t normHyp { normASqSum * normCoshTerm + normInterTermRe * normSinhTerm }; const Double_t normTrig { - prodAsym * ( normASqDiff * normCosTerm + normInterTermIm * normSinTerm ) }; // Combine all terms to get the total normalisation const Double_t norm { 2.0 * ( normHyp + normTrig ) }; // Multiply the squared-amplitude by the efficiency (DP and decay time) and decay-time error likelihood // and normalise to obtain the signal likelihood sigDPLike_ = ( ASq * dpEff * dtEff * dtErrLike ) / norm; // Background part // TODO move to new function as getEvtBkgndLikelihoods? const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { if ( not usingBkgnd_ ) { bkgndDPLike_[bkgndID] = 0.0; continue; } Double_t omegaBkgnd{1.0}; Double_t omegaBarBkgnd{1.0}; BkgndDecayTimePdfs_[bkgndID]->calcLikelihoodInfo(static_cast(iEvt)); // Consider background type switch ( BkgndTypes_[bkgndID] ) { case LauFlavTag::BkgndType::Combinatorial : { // For combinatorial background the DP and decay-time models factorise completely, just mulitply them // Start with the DP likelihood... if ( (cpEigenValue_ == QFS) and BkgndDPModelsBbar_[bkgndID] != nullptr ) { //Flavour specific (with possible detection asymmetry) if ( curEvtDecayFlv_ == LauFlavTag::Flavour::B ) { bkgndDPLike_[bkgndID] = BkgndDPModelsB_[bkgndID]->getUnNormValue(iEvt); } else { bkgndDPLike_[bkgndID] = BkgndDPModelsBbar_[bkgndID]->getUnNormValue(iEvt); } bkgndDPLike_[bkgndID] /= ( BkgndDPModelsB_[bkgndID]->getPdfNorm() + BkgndDPModelsBbar_[bkgndID]->getPdfNorm() ); } else { bkgndDPLike_[bkgndID] = 0.5 * BkgndDPModelsB_[bkgndID]->getLikelihood(iEvt); } // ...include the decay time... switch( BkgndDecayTimePdfs_[bkgndID]->getFuncType() ) { case LauDecayTime::FuncType::Hist : bkgndDPLike_[bkgndID] *= BkgndDecayTimePdfs_[bkgndID]->getHistTerm(); break; case LauDecayTime::FuncType::Exp : bkgndDPLike_[bkgndID] *= ( BkgndDecayTimePdfs_[bkgndID]->getExpTerm() / BkgndDecayTimePdfs_[bkgndID]->getNormTermExp() ); break; // TODO - any other decay time function types that make sense for combinatorial? // - should also have a set of checks in initialise that we have everything we need for the backgrounds and that the various settings make sense default : // TODO as per comment just above, once the above mentioned checks are implemented this error message can be removed std::cerr << "WARNING in LauTimeDepFitModel::getEvtDPDtLikelihood : bkgnd types other than Hist and Exp don't make sense for combinatorial!" << std::endl; break; } // ...include flavour tagging for (std::size_t taggerID{0}; taggerIDgetCapitalOmegaBkgnd(bkgndID, taggerID, LauFlavTag::Flavour::B, curEvtDecayFlv_); } bkgndDPLike_[bkgndID] *= omegaBkgnd; break; } case LauFlavTag::BkgndType::FlavourSpecific : { // DP terms needed by all decay-time cases Double_t Asq { BkgndDPModelsB_[bkgndID]->getUnNormValue(iEvt) }; Double_t Asqbar { BkgndDPModelsBbar_[bkgndID]->getUnNormValue(iEvt) }; if ( cpEigenValue_ == QFS ){ // If the signal is flavour-specific we can know which DP to use, so zero the other one if ( curEvtDecayFlv_ == LauFlavTag::Flavour::B ) { Asqbar = 0.0; } else if ( curEvtDecayFlv_ == LauFlavTag::Flavour::Bbar ) { Asq = 0.0; } } const Double_t AsqSum { Asq + Asqbar }; // DP norm terms needed by all decay-time cases const Double_t AsqNorm { BkgndDPModelsB_[bkgndID]->getPdfNorm() }; const Double_t AsqbarNorm { BkgndDPModelsBbar_[bkgndID]->getPdfNorm() }; const Double_t AsqNormSum { AsqNorm + AsqbarNorm }; // FT terms needed by all decay-time cases omegaBkgnd = omegaBarBkgnd = 1.0; for (std::size_t taggerID{0}; taggerIDgetCapitalOmegaBkgnd(bkgndID, taggerID, LauFlavTag::Flavour::B, curEvtDecayFlv_); omegaBarBkgnd *= flavTag_->getCapitalOmegaBkgnd(bkgndID, taggerID, LauFlavTag::Flavour::Bbar, curEvtDecayFlv_); } const Double_t AProd { AProdBkgnd_[bkgndID]->unblindValue() }; const Double_t ftOmegaHypBkgnd { (1.0 - AProd)*omegaBkgnd + (1.0 + AProd)*omegaBarBkgnd }; switch( BkgndDecayTimePdfs_[bkgndID]->getFuncType() ) { case LauDecayTime::FuncType::Hist: // DP and decay-time still factorise { // Start with the DP terms... bkgndDPLike_[bkgndID] = AsqSum / AsqNormSum; // ...include the decay time... bkgndDPLike_[bkgndID] *= BkgndDecayTimePdfs_[bkgndID]->getHistTerm(); // ...include flavour tagging bkgndDPLike_[bkgndID] *= ( 0.5 * ftOmegaHypBkgnd ); break; } case LauDecayTime::FuncType::Exp : // DP and decay-time still factorise { // Start with the DP terms... bkgndDPLike_[bkgndID] = AsqSum / AsqNormSum; // ...include the decay time... bkgndDPLike_[bkgndID] *= ( BkgndDecayTimePdfs_[bkgndID]->getExpTerm() / BkgndDecayTimePdfs_[bkgndID]->getNormTermExp() ); // ...include flavour tagging bkgndDPLike_[bkgndID] *= ( 0.5 * ftOmegaHypBkgnd ); break; } case LauDecayTime::FuncType::ExpTrig: // DP and decay-time don't factorise case LauDecayTime::FuncType::ExpHypTrig: { // DP and FT terms specific to this case const Double_t AsqDiff { Asq - Asqbar }; const Double_t AsqNormDiff { AsqNorm - AsqbarNorm }; //TODO check this shouldn't be `fabs`ed const Double_t ftOmegaTrigBkgnd { (1.0 - AProd)*omegaBkgnd - (1.0 + AProd)*omegaBarBkgnd }; // decay time terms // Sin and Sinh terms are ignored: the FS modes can't exhibit TD CPV const Double_t dtCoshBkgnd { BkgndDecayTimePdfs_[bkgndID]->getCoshTerm() }; const Double_t dtCosBkgnd { BkgndDecayTimePdfs_[bkgndID]->getCosTerm() }; const Double_t normCoshTermBkgnd { BkgndDecayTimePdfs_[bkgndID]->getNormTermCosh() }; const Double_t normCosTermBkgnd { BkgndDecayTimePdfs_[bkgndID]->getNormTermCos() }; // Combine the DP, FT, and decay time terms const Double_t coshTermBkgnd { ftOmegaHypBkgnd * dtCoshBkgnd * AsqSum }; const Double_t cosTermBkgnd { ftOmegaTrigBkgnd * dtCosBkgnd * AsqDiff }; // Similarly for the normalisation, see Laura note eq. 41 const Double_t normBkgnd { 2.0 * ( (normCoshTermBkgnd * AsqNormSum) - AProd*(normCosTermBkgnd * AsqNormDiff) ) }; bkgndDPLike_[bkgndID] = (coshTermBkgnd + cosTermBkgnd) / normBkgnd; break; } case LauDecayTime::FuncType::Delta: case LauDecayTime::FuncType::DeltaExp: // TODO as per comment above, once the checks in initialise are implemented this error message can be removed std::cerr << "WARNING in LauTimeDepFitModel::getEvtDPDtLikelihood : bkgnd types Delta and DeltaExp don't make sense!" << std::endl; break; } break; } case LauFlavTag::BkgndType::SelfConjugate : //Copy this from the CPeigenstate signal case std::cerr << "WARNING in LauTimeDepFitModel::getEvtDPDtLikelihood : SelfConjugate states aren't implemented yet!" << std::endl; bkgndDPLike_[bkgndID] = 0.0; break; case LauFlavTag::BkgndType::NonSelfConjugate : // TODO this has been ignored for now since it's not used in the B->Dpipi case std::cerr << "WARNING in LauTimeDepFitModel::getEvtDPDtLikelihood : NonSelfConjugate states aren't implemented yet!" << std::endl; bkgndDPLike_[bkgndID] = 0.0; break; } // Get the decay time acceptance const Double_t dtEffBkgnd { BkgndDecayTimePdfs_[bkgndID]->getEffiTerm() }; // Get the decay time error term const Double_t dtErrLikeBkgnd { BkgndDecayTimePdfs_[bkgndID]->getErrTerm() }; // Include these terms in the background likelihood bkgndDPLike_[bkgndID] *= ( dtEffBkgnd * dtErrLikeBkgnd ); } } 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 if ( not sigExtraPdf_.empty() ) { sigExtraLike_ = this->prodPdfValue( sigExtraPdf_, iEvt ); } // Background const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { if (usingBkgnd_) { bkgndExtraLike_[bkgndID] = this->prodPdfValue( BkgndPdfs_[bkgndID], iEvt ); } else { bkgndExtraLike_[bkgndID] = 0.0; } } } void LauTimeDepFitModel::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(const TString& fileName, const TString& treeName, const Bool_t reuseEventsWithinEnsemble, const Bool_t reuseEventsWithinExperiment, const Bool_t useReweighting) { if (signalTree_) { std::cerr<<"ERROR in LauTimeDepFitModel::embedSignal : Already embedding signal from file."<findBranches(); if (!dataOK) { delete signalTree_; signalTree_ = 0; std::cerr<<"ERROR in LauTimeDepFitModel::embedSignal : Problem creating data tree for embedding."<enableEmbedding(kTRUE); } void LauTimeDepFitModel::embedBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName, const Bool_t reuseEventsWithinEnsemble, const Bool_t reuseEventsWithinExperiment, const Bool_t useReweighting) { if ( ! this->validBkgndClass( bkgndClass ) ) { std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; return; } UInt_t bkgndID = this->bkgndClassID( bkgndClass ); LauEmbeddedData* bkgTree = bkgndTree_[bkgndID]; if (bkgTree) { std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Already embedding background from a file." << std::endl; return; } bkgTree = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment); Bool_t dataOK = bkgTree->findBranches(); if (!dataOK) { delete bkgTree; bkgTree = 0; std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Problem creating data tree for embedding." << std::endl; return; } reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble; useReweighting_ = useReweighting; this->enableEmbedding(kTRUE); } void LauTimeDepFitModel::setupSPlotNtupleBranches() { // add branches for storing the experiment number and the number of // the event within the current experiment this->addSPlotNtupleIntegerBranch("iExpt"); this->addSPlotNtupleIntegerBranch("iEvtWithinExpt"); // Store the efficiency of the event (for inclusive BF calculations). if (this->storeDPEff()) { this->addSPlotNtupleDoubleBranch("efficiency"); } // Store the total event likelihood for each species. this->addSPlotNtupleDoubleBranch("sigTotalLike"); if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name( this->bkgndClassName(iBkgnd) ); name += "TotalLike"; this->addSPlotNtupleDoubleBranch(name); } } // Store the DP likelihoods if (this->useDP()) { this->addSPlotNtupleDoubleBranch("sigDPLike"); if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name( this->bkgndClassName(iBkgnd) ); name += "DPLike"; this->addSPlotNtupleDoubleBranch(name); } } } // Store the likelihoods for each extra PDF this->addSPlotNtupleBranches(sigExtraPdf_, "sig"); if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); this->addSPlotNtupleBranches(BkgndPdfs_[iBkgnd], bkgndClass); } } } void LauTimeDepFitModel::addSPlotNtupleBranches(const LauPdfPList& extraPdfs, const TString& prefix) { // Loop through each of the PDFs for ( const LauAbsPdf* pdf : extraPdfs ) { // 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}; const std::vector varNames { pdf->varNames() }; for ( const TString& varName : varNames ) { if ( varName != "m13Sq" && varName != "m23Sq" ) { ++nVars; } } if ( nVars == 1 ) { // If the PDF only has one variable then // simply add one branch for that variable TString name{prefix}; name += pdf->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 ( const TString& varName : varNames ) { if ( varName != "m13Sq" && varName != "m23Sq" ) { allVars += varName; TString name{prefix}; name += varName; 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."<calcLikelihoodInfo(iEvt); extraLike = pdf->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}; const std::vector varNames { pdf->varNames() }; for ( const TString& varName : varNames ) { if ( varName != "m13Sq" && varName != "m23Sq" ) { ++nVars; } } if ( nVars == 1 ) { // If the PDF only has one variable then // simply store the value for that variable TString name{prefix}; name += pdf->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 ( const TString& varName : varNames ) { if ( varName != "m13Sq" && varName != "m23Sq" ) { allVars += varName; TString name{prefix}; name += varName; name += "Like"; const Double_t indivLike = pdf->getLikelihood( varName ); this->setSPlotNtupleDoubleBranchValue(name, indivLike); } } TString name{prefix}; name += allVars; name += "Like"; this->setSPlotNtupleDoubleBranchValue(name, extraLike); } else { std::cerr<<"WARNING in LauAllFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs."<useDP()) { nameSet.insert("DP"); } for ( const LauAbsPdf* pdf : sigExtraPdf_ ) { // Loop over the variables involved in each PDF const std::vector varNames { pdf->varNames() }; for ( const TString& varName : varNames ) { // If they are not DP coordinates then add them if ( varName != "m13Sq" && varName != "m23Sq" ) { nameSet.insert( varName ); } } } return nameSet; } LauSPlot::NumbMap LauTimeDepFitModel::freeSpeciesNames() const { LauSPlot::NumbMap numbMap; if (!signalEvents_->fixed() && this->doEMLFit()) { numbMap["sig"] = signalEvents_->genValue(); } if ( usingBkgnd_ ) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); const LauAbsRValue* par = bkgndEvents_[iBkgnd]; if (!par->fixed()) { numbMap[bkgndClass] = par->genValue(); if ( ! par->isLValue() ) { std::cerr << "WARNING in LauTimeDepFitModel::freeSpeciesNames : \"" << par->name() << "\" is a LauFormulaPar, which implies it is perhaps not entirely free to float in the fit, so the sWeight calculation may not be reliable" << std::endl; } } } } return numbMap; } LauSPlot::NumbMap LauTimeDepFitModel::fixdSpeciesNames() const { LauSPlot::NumbMap numbMap; if (signalEvents_->fixed() && this->doEMLFit()) { numbMap["sig"] = signalEvents_->genValue(); } if ( usingBkgnd_ ) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); const LauAbsRValue* par = bkgndEvents_[iBkgnd]; if (par->fixed()) { numbMap[bkgndClass] = par->genValue(); } } } return numbMap; } LauSPlot::TwoDMap LauTimeDepFitModel::twodimPDFs() const { LauSPlot::TwoDMap twodimMap; for ( const LauAbsPdf* pdf : sigExtraPdf_ ) { // Count the number of input variables that are not DP variables UInt_t nVars{0}; const std::vector varNames { pdf->varNames() }; for ( const TString& varName : varNames ) { if ( varName != "m13Sq" && varName != "m23Sq" ) { ++nVars; } } if ( nVars == 2 ) { twodimMap.insert( std::make_pair( "sig", std::make_pair( varNames[0], varNames[1] ) ) ); } } if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); for ( const LauAbsPdf* pdf : BkgndPdfs_[iBkgnd] ) { // Count the number of input variables that are not DP variables UInt_t nVars{0}; const std::vector varNames { pdf->varNames() }; for ( const TString& varName : varNames ) { if ( varName != "m13Sq" && varName != "m23Sq" ) { ++nVars; } } if ( nVars == 2 ) { twodimMap.insert( std::make_pair( bkgndClass, std::make_pair( varNames[0], varNames[1] ) ) ); } } } } return twodimMap; } void LauTimeDepFitModel::storePerEvtLlhds() { std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Storing per-event likelihood values..."<fitData(); // if we've not been using the DP model then we need to cache all // the info here so that we can get the efficiency from it if (!this->useDP() && this->storeDPEff()) { sigModelB0bar_->initialise(coeffsB0bar_); sigModelB0_->initialise(coeffsB0_); sigModelB0bar_->fillDataTree(*inputFitData); sigModelB0_->fillDataTree(*inputFitData); } UInt_t evtsPerExpt(this->eventsPerExpt()); LauIsobarDynamics* sigModel(sigModelB0bar_); for (UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt) { // Find out whether we have B0bar or B0 flavTag_->updateEventInfo(iEvt); curEvtTagFlv_ = flavTag_->getCurEvtTagFlv(); curEvtMistag_ = flavTag_->getCurEvtMistag(); // the DP information this->getEvtDPDtLikelihood(iEvt); if (this->storeDPEff()) { if (!this->useDP()) { sigModel->calcLikelihoodInfo(iEvt); } this->setSPlotNtupleDoubleBranchValue("efficiency",sigModel->getEvtEff()); } if (this->useDP()) { sigTotalLike_ = sigDPLike_; this->setSPlotNtupleDoubleBranchValue("sigDPLike",sigDPLike_); if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name = this->bkgndClassName(iBkgnd); name += "DPLike"; this->setSPlotNtupleDoubleBranchValue(name,bkgndDPLike_[iBkgnd]); } } } else { sigTotalLike_ = 1.0; } // the signal PDF values sigTotalLike_ *= this->setSPlotNtupleBranchValues(sigExtraPdf_, "sig", iEvt); // the background PDF values if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); LauPdfPList& pdfs = BkgndPdfs_[iBkgnd]; bkgndTotalLike_[iBkgnd] *= this->setSPlotNtupleBranchValues(pdfs, bkgndClass, iEvt); } } // the total likelihoods this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_); if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name = this->bkgndClassName(iBkgnd); name += "TotalLike"; this->setSPlotNtupleDoubleBranchValue(name,bkgndTotalLike_[iBkgnd]); } } // fill the tree this->fillSPlotNtupleBranches(); } std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Finished storing per-event likelihood values."<