diff --git a/inc/LauFlavTag.hh b/inc/LauFlavTag.hh index 4b343a2..2445ad1 100644 --- a/inc/LauFlavTag.hh +++ b/inc/LauFlavTag.hh @@ -1,542 +1,563 @@ /* 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 "LauParameter.hh" #include "TString.h" -#include "LauParameter.hh" +#include + +#include +#include class LauAbsPdf; class LauFitDataTree; +class TH1; 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 + //! Add a tagger /*! \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 + //! Add a tagger + /*! + \param [in] config the JSON object containing the tagging variable name, the mistag variable name, the tagging efficiency, calib_p0, and calib_p1 information + \param [in] etapdf the mistag distribution for the tagger + */ + void addTagger(const nlohmann::json& config, LauAbsPdf* etapdf); + + //! Add a tagger /*! \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); + //! Add a tagger + /*! + \param [in] config the JSON object containing the tagging variable name, the mistag variable name, the tagging efficiency, calib_p0, and calib_p1 information + \param [in] etapdf the mistag distribution for the tagger + \param [in] tagEff tagging efficiency histograms - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag + */ + void addTagger(const nlohmann::json& config, LauAbsPdf* etapdf, const std::pair tagEff); + //! 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 \param [in] decayFlv the decay flavour of the B */ 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(); + //! Write tagger information to JSON + void writeTaggersToJson(const TString& jsonFileName) const; + 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 }; //! \cond DOXYGEN_IGNORE // map LauKMatrixPropagator::KMatrixChannels values to JSON as strings LAURA_NLOHMANN_JSON_SERIALIZE_ENUM( LauFlavTag::BkgndType, { {LauFlavTag::BkgndType::Combinatorial, "Combinatorial"}, {LauFlavTag::BkgndType::FlavourSpecific, "FlavourSpecific"}, {LauFlavTag::BkgndType::SelfConjugate, "SelfConjugate"}, {LauFlavTag::BkgndType::NonSelfConjugate, "NonSelfConjugate"}, }) //! \endcond #endif diff --git a/src/LauFlavTag.cc b/src/LauFlavTag.cc index 4f86f60..c45df27 100644 --- a/src/LauFlavTag.cc +++ b/src/LauFlavTag.cc @@ -1,1445 +1,1589 @@ /* 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 "LauFlavTag.hh" + +#include "Lau1DHistPdf.hh" +#include "LauAbsPdf.hh" +#include "LauJsonTools.hh" +#include "LauRandom.hh" #include "TMath.h" #include "TString.h" #include "TSystem.h" -#include "Lau1DHistPdf.hh" -#include "LauAbsPdf.hh" -#include "LauFlavTag.hh" -#include "LauRandom.hh" +#include + +#include +#include +#include ClassImp(LauFlavTag) 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 nlohmann::json& config, LauAbsPdf* etapdf) +{ + using LauJsonTools::JsonType; + using LauJsonTools::ElementNameType; + using LauJsonTools::checkObjectElements; + using LauJsonTools::getValue; + + if ( ! config.is_object() ) { + std::cerr << "ERROR in LauFlavTag::addTagger : invalid JSON object supplied" << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + + const std::string suffix1 { useAveDelta_ ? "_ave" : "_b0" }; + const std::string suffix2 { useAveDelta_ ? "_delta" : "_b0bar" }; + + const std::vector mandatoryElements { + std::make_pair("name", JsonType::String), + std::make_pair("dec_name", JsonType::String), + std::make_pair("eta_name", JsonType::String), + std::make_pair("eff"+suffix1, JsonType::Number), + std::make_pair("eff"+suffix2, JsonType::Number), + std::make_pair("calib_p0"+suffix1, JsonType::Number), + std::make_pair("calib_p0"+suffix2, JsonType::Number), + std::make_pair("calib_p1"+suffix1, JsonType::Number), + std::make_pair("calib_p1"+suffix2, JsonType::Number) + }; + if ( ! checkObjectElements( config, mandatoryElements ) ) { + std::cerr << "ERROR in LauFlavTag::addTagger : invalid JSON object supplied" << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + + // Extract the necessary information from the JSON object + const TString name { getValue( config, "name" ) }; + const TString tagVarName { getValue( config, "dec_name" ) }; + const TString mistagVarName { getValue( config, "eta_name" ) }; + + auto tagEff = std::make_pair( + getValue( config, "eff"+suffix1 ), + getValue( config, "eff"+suffix2 ) ); + + auto calib_p0 = std::make_pair( + getValue( config, "calib_p0"+suffix1 ), + getValue( config, "calib_p0"+suffix2 ) ); + + auto calib_p1 = std::make_pair( + getValue( config, "calib_p1"+suffix1 ), + getValue( config, "calib_p1"+suffix2 ) ); + + // Self-message to actually add the tagger + this->addTagger( name, tagVarName, mistagVarName, etapdf, tagEff, calib_p0, calib_p1 ); +} + 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::addTagger(const nlohmann::json& config, LauAbsPdf* etapdf, const std::pair tagEff) +{ + using LauJsonTools::JsonType; + using LauJsonTools::ElementNameType; + using LauJsonTools::checkObjectElements; + using LauJsonTools::getValue; + + if ( ! config.is_object() ) { + std::cerr << "ERROR in LauFlavTag::addTagger : invalid JSON object supplied" << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + + const std::string suffix1 { useAveDelta_ ? "_ave" : "_b0" }; + const std::string suffix2 { useAveDelta_ ? "_delta" : "_b0bar" }; + + const std::vector mandatoryElements { + std::make_pair("name", JsonType::String), + std::make_pair("dec_name", JsonType::String), + std::make_pair("eta_name", JsonType::String), + std::make_pair("calib_p0"+suffix1, JsonType::Number), + std::make_pair("calib_p0"+suffix2, JsonType::Number), + std::make_pair("calib_p1"+suffix1, JsonType::Number), + std::make_pair("calib_p1"+suffix2, JsonType::Number) + }; + if ( ! checkObjectElements( config, mandatoryElements ) ) { + std::cerr << "ERROR in LauFlavTag::addTagger : invalid JSON object supplied" << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + + // Extract the necessary information from the JSON object + const TString name { getValue( config, "name" ) }; + const TString tagVarName { getValue( config, "dec_name" ) }; + const TString mistagVarName { getValue( config, "eta_name" ) }; + + auto calib_p0 = std::make_pair( + getValue( config, "calib_p0"+suffix1 ), + getValue( config, "calib_p0"+suffix2 ) ); + + auto calib_p1 = std::make_pair( + getValue( config, "calib_p1"+suffix1 ), + getValue( config, "calib_p1"+suffix2 ) ); + + // Self-message to actually add the tagger + this->addTagger( name, tagVarName, mistagVarName, etapdf, tagEff, calib_p0, calib_p1 ); +} + 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(0.5); 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,0.0,1.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,2.0,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,0.0,1.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,2.0,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,0.0,1.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,2.0,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,-1.0,1.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,-1.0,1.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() { 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(); } } } } + +void LauFlavTag::writeTaggersToJson(const TString& jsonFileName) const +{ + using nlohmann::json; + + json config = json::array(); + + for ( auto& [ name, index ] : taggerIndex_ ) { + json taggerConf = json::object(); + + taggerConf["name"] = name; + taggerConf["dec_name"] = tagVarNames_[index]; + taggerConf["eta_name"] = mistagVarNames_[index]; + + if ( useAveDelta_ ) { + if ( tagEff_ave_[index] != nullptr && tagEff_delta_[index] != nullptr ) { + taggerConf["eff_ave"] = tagEff_ave_[index]->value(); + taggerConf["eff_delta"] = tagEff_delta_[index]->value(); + } + taggerConf["calib_p0_ave"] = calib_p0_ave_[index]->value(); + taggerConf["calib_p0_delta"] = calib_p0_delta_[index]->value(); + taggerConf["calib_p1_ave"] = calib_p1_ave_[index]->value(); + taggerConf["calib_p1_delta"] = calib_p1_delta_[index]->value(); + } else { + if ( tagEff_B0_[index] != nullptr && tagEff_B0bar_[index] != nullptr ) { + taggerConf["eff_b0"] = tagEff_B0_[index]->value(); + taggerConf["eff_b0bar"] = tagEff_B0bar_[index]->value(); + } + taggerConf["calib_p0_b0"] = calib_p0_B0_[index]->value(); + taggerConf["calib_p0_b0bar"] = calib_p0_B0bar_[index]->value(); + taggerConf["calib_p1_b0"] = calib_p1_B0_[index]->value(); + taggerConf["calib_p1_b0bar"] = calib_p1_B0bar_[index]->value(); + } + + config.push_back( taggerConf ); + } + + const bool writeOK { LauJsonTools::writeJsonFile( config, jsonFileName.Data() ) }; + if ( ! writeOK ) { + std::cerr << "ERROR in LauFlavTag::writeTaggersToJson : could not successfully write to file \"" << jsonFileName << std::endl; + } +}