diff --git a/inc/LauTimeDepFitModel.hh b/inc/LauTimeDepFitModel.hh index 532341f..c4bddb2 100644 --- a/inc/LauTimeDepFitModel.hh +++ b/inc/LauTimeDepFitModel.hh @@ -1,721 +1,724 @@ /* Copyright 2006 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauTimeDepFitModel.hh \brief File containing declaration of LauTimeDepFitModel class. */ /*! \class LauTimeDepFitModel \brief Class for defining a time-dependent fit model. LauTimeDepFitModel is a class that allows the user to define a three-body Dalitz plot according to the isobar model, i.e. defining a set of resonances that have complex amplitudes that can interfere with each other. It extends the LauSimpleFitModel and LauCPFitModel models in that it allows the fitting of time-dependent particle/antiparticle decays to flavour-conjugate Dalitz plots, including their interference through mixing. */ #ifndef LAU_TIMEDEP_FIT_MODEL #define LAU_TIMEDEP_FIT_MODEL #include #include #include #include #include "TString.h" #include "TStopwatch.h" #include "TSystem.h" #include "LauAbsFitModel.hh" #include "LauConstants.hh" #include "LauEmbeddedData.hh" #include "LauParameter.hh" #include "LauFlavTag.hh" #include "LauCategoryFlavTag.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); //! Set the signal event yield and asymmetry /*! \param [in] nSigEvents contains the signal yield and option to fix it \param [in] sigAsym contains the signal asymmetry and option to fix it */ virtual void setNSigEvents(LauParameter* nSigEvents, LauParameter* sigAsym); //! Set the number of background events /*! The name of the parameter must be that of the corresponding background category (so that it can be correctly assigned) \param [in] nBkgndEvents contains the name, yield and option to fix the yield of the background */ virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents); //! Set the background event yield and asymmetry /*! \param [in] nBkgEvents contains the background yield and option to fix it \param [in] BkgAsym contains the background asymmetry and option to fix it */ virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym); //! Set the background DP models (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 evVarName the variable name in the data tree that specifies the CP eigenvalue */ inline void setCPEigenvalue( const CPEigenvalue defaultCPEV, const TString& cpevVarName = "" ) { cpEigenValue_ = defaultCPEV; cpevVarName_ = cpevVarName; } //! Set the DP amplitude coefficients /*! \param [in] coeffSet the set of coefficients */ void setAmpCoeffSet(LauAbsCoeffSet* coeffSet); //! Set the decay time PDFs /*! \param [in] position the tagger position in the vectors \param [in] pdf the signal decay time PDF */ void setSignalDtPdf(LauDecayTimePdf* pdf); //! Set the decay time PDFs /*! \param [in] tagCat the tagging category for which the PDF should be used \param [in] pdf the background decay time PDF */ void setBkgndDtPdf(const TString& bkgndClass, LauDecayTimePdf* pdf); //! Set the signal PDF for a given variable /*! \param [in] tagCat the tagging category for which the PDF should be used \param [in] pdf the PDF to be added to the signal model */ void setSignalPdfs(LauAbsPdf* pdf); //! Set the background PDF /*! \param [in] bkgndClass the name of the background class \param [in] pdf the PDF to be added to the background model */ void setBkgndPdf(const TString& bkgndClass, LauAbsPdf* pdf); void setSignalFlavTagPdfs( const Int_t tagCat, LauAbsPdf* pdf); void setBkgdFlavTagPdfs( const TString name, LauAbsPdf* pdf); //! Embed full simulation events for the signal, rather than generating toy from the PDFs /*! \param [in] tagCat the tagging category for which the file should be used \param [in] fileName the name of the file containing the events \param [in] treeName the name of the tree \param [in] reuseEventsWithinEnsemble \param [in] reuseEventsWithinExperiment \param [in] useReweighting */ void embedSignal(const TString& fileName, const TString& treeName, const Bool_t reuseEventsWithinEnsemble, const Bool_t reuseEventsWithinExperiment = kFALSE); void embedBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName, Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment = kFALSE); //! Set the value of the mixing phase /*! \param [in] phiMix the value of the mixing phase \param [in] fixPhiMix whether the value should be fixed or floated \param [in] useSinCos whether to use the sine and cosine as separate parameters or to just use the mixing phase itself */ void setPhiMix(const Double_t phiMix, const Bool_t fixPhiMix, const Bool_t useSinCos = kFALSE); //! Initialise the fit virtual void initialise(); //! Initialise the signal DP model virtual void initialiseDPModels(); //! Recalculate Normalization the signal DP models virtual void recalculateNormalisation(); //! Update the coefficients virtual void updateCoeffs(); // Toy MC generation and fitting overloaded functions virtual Bool_t genExpt(); //! Set the maximum value of A squared to be used in the accept/reject /*! \param [in] value the new value */ inline void setASqMaxValue(const Double_t value) {aSqMaxSet_ = value;} //! Weight events based on the DP model /*! \param [in] dataFileName the name of the data file \param [in] dataTreeName the name of the data tree */ virtual void weightEvents( const TString& dataFileName, const TString& dataTreeName ); //! Calculate things that depend on the fit parameters after they have been updated by Minuit virtual void propagateParUpdates(); //! Read in the input fit data variables, e.g. m13Sq and m23Sq virtual void cacheInputFitVars(); //! Check the initial fit parameters virtual void checkInitFitParams(); //! Get the fit results and store them /*! \param [in] tablePrefixName prefix for the name of the output file */ virtual void finaliseFitResults(const TString& tablePrefixName); //! Save the pdf Plots for all the resonances of experiment number fitExp /*! TODO - not working in this model!! \param [in] label prefix for the file name to be saved */ virtual void savePDFPlots(const TString& label); //! Save the pdf Plots for the sum of ressonances correspondint to "sin" of experiment number fitExp /*! TODO - not working in this model!! \param [in] label prefix for the file name to be saved \param [in] spin spin of the wave to be saved */ virtual void savePDFPlotsWave(const TString& label, const Int_t& spin); //! Print the fit fractions, total DP rate and mean efficiency /*! \param [out] output the stream to which to print */ virtual void printFitFractions(std::ostream& output); //! Print the asymmetries /*! \param [out] output the stream to which to print */ virtual void printAsymmetries(std::ostream& output); //! Write the fit results in latex table format /*! \param [in] outputFile the name of the output file */ virtual void writeOutTable(const TString& outputFile); //! Store the per event likelihood values virtual void storePerEvtLlhds(); // Methods to do with calculating the likelihood functions // and manipulating the fitting parameters. //! Get the total likelihood for each event /*! \param [in] iEvt the event number */ virtual Double_t getTotEvtLikelihood(const UInt_t iEvt); //! Calculate the signal and background likelihoods for the DP for a given event /*! \param [in] iEvt the event number */ virtual void getEvtDPDtLikelihood(const UInt_t iEvt); //! Determine the signal and background likelihood for the extra variables for a given event /*! \param [in] iEvt the event number */ virtual void getEvtExtraLikelihoods(const UInt_t iEvt); virtual void getEvtFlavTagLikelihood(const UInt_t iEvt); //! Get the total number of events /*! \return the total number of events */ virtual Double_t getEventSum() const; //! Set the fit parameters for the DP model void setSignalDPParameters(); //! Set the fit parameters for the decay time PDFs void setDecayTimeParameters(); //! Set the fit parameters for the extra PDFs void setExtraPdfParameters(); //! Set the initial yields void setFitNEvents(); //! Set the calibration parameters void setCalibParams(); //! Set the tagging efficiency parameters void setTagEffParams(); //! Set the efficiency parameters void setEffiParams(); //! Set the asymmetry parameters void setAsymParams(); //! Set the tagging asymmetry parameters void setFlavTagAsymParams(); //! Set-up other parameters that are derived from the fit results, e.g. fit fractions void setExtraNtupleVars(); //! Set production and detections asymmetries + /*! + \param [in] AProd is the production asymmetry definied as (N_BqBar-N_Bq)/(N_BqBar+N_Bq) + */ void setAsymmetries(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(); //! Return the map of signal decay time PDFs typedef std::map< Int_t, LauDecayTimePdf*> LauTagCatDtPdfMap; LauDecayTimePdf* getSignalDecayTimePdf(){return signalDecayTimePdf_;} //! Return the map of background decay time PDFs std::vector getBackgroundDecayTimePdfs(){return BkgndDecayTimePdfs_;} protected: typedef std::map< Int_t, LauParameter> LauTagCatParamMap; typedef std::map< Int_t, LauPdfList > LauTagCatPdfListMap; typedef std::map< Int_t, LauAbsPdf* > LauTagCatPdfMap; typedef std::map< TString, LauAbsPdf* > LauBkgdPdfMap; typedef std::map< Int_t, Int_t > LauTagCatEventsMap; typedef std::map< Int_t, LauEmbeddedData* > LauTagCatEmbDataMap; //typedef std::map< Int_t, std::pair > LauTaggerGenInfo; //typedef std::map< std::pair, LauTaggerGenInfo > LauGenInfo; typedef std::map< TString, std::pair > LauGenInfo; typedef std::vector LauTagCatEmbDataMapList; typedef std::vector LauBkgndDPModelList; typedef std::vector LauBkgndPdfsList; typedef std::vector LauBkgndYieldList; typedef std::vector LauBkgndReuseEventsList; //! Determine the number of events to generate for each hypothesis LauGenInfo eventsToGenerate(); //! Generate signal event Bool_t generateSignalEvent(); //! Generate background event /*! \param [in] bgID ID number of the background class */ Bool_t generateBkgndEvent(UInt_t bgID); //! Setup the required ntuple branches void setupGenNtupleBranches(); //! Store all of the DP and decay time information void setDPDtBranchValues(); //! Generate from the extra PDFs /*! \param [in] extraPdfs the list of extra PDFs \param [in] embeddedData the embedded data sample */ void generateExtraPdfValues(LauPdfList& extraPdfs, LauEmbeddedData* embeddedData); //! Add sPlot branches for the extra PDFs /*! \param [in] extraPdfs the list of extra PDFs \param [in] prefix the list of prefixes for the branch names */ void addSPlotNtupleBranches(const LauPdfList& extraPdfs, const TString& prefix); //! Set the branches for the sPlot ntuple with extra PDFs /*! \param [in] extraPdfs the list of extra PDFs \param [in] prefix the list of prefixes for the branch names \param [in] iEvt the event number */ Double_t setSPlotNtupleBranchValues(LauPdfList& extraPdfs, const TString& prefix, 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; } //! Add the parameters from each PDF into the fit /*! \param [in] theMap the container of PDFs */ UInt_t addParametersToFitList(LauPdfList* theList); //! Add the parameters from each decay time PDF into the fit /*! \param [in] theMap the container of PDFs */ UInt_t addParametersToFitList(std::vector theVector); //! Calculate the component integrals of the interference term void calcInterferenceTermIntegrals(); //! Calculate the total integral of the interference term void calcInterTermNorm(); private: //! Dalitz plot PDF for the antiparticle LauIsobarDynamics* sigModelB0bar_; //! Dalitz plot PDF for the particle LauIsobarDynamics* sigModelB0_; //! Kinematics object for antiparticle LauKinematics* kinematicsB0bar_; //! Kinematics object for particle LauKinematics* kinematicsB0_; //! The 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 LauPdfList sigExtraPdf_; //! eta PDFs for each TagCat LauPdfList sigFlavTagPdf_; //! eta PDFs for each background LauBkgdPdfMap bkgdFlavTagPdf_; //! Production asymmetry between B0 and B0bar LauParameter AProd_; // 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_; //! Vector of booleans to reuse background events LauBkgndReuseEventsList reuseBkgnd_; // Likelihood values //! Signal DP likelihood value Double_t sigDPLike_; //! Signal likelihood from extra PDFs Double_t sigExtraLike_; Double_t sigFlavTagLike_; Double_t bkgdFlavTagLike_; //! Total signal likelihood Double_t sigTotalLike_; //! Background DP likelihood value(s) 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/Lau1DHistPdf.cc b/src/Lau1DHistPdf.cc index 23a6cc6..cc54568 100644 --- a/src/Lau1DHistPdf.cc +++ b/src/Lau1DHistPdf.cc @@ -1,267 +1,268 @@ /* 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 Lau1DHistPdf.cc \brief File containing implementation of Lau1DHistPdf class. */ #include #include #include #include "TAxis.h" #include "TH1.h" #include "TRandom.h" #include "TSystem.h" #include "Lau1DHistPdf.hh" #include "LauRandom.hh" class LauParameter; ClassImp(Lau1DHistPdf) Lau1DHistPdf::Lau1DHistPdf(const TString& theVarName, const TH1* hist, Double_t minAbscissa, Double_t maxAbscissa, Bool_t useInterpolation, Bool_t fluctuateBins) : LauAbsPdf(theVarName, std::vector(), minAbscissa, maxAbscissa), hist_(hist ? dynamic_cast(hist->Clone()) : 0), useInterpolation_(useInterpolation), fluctuateBins_(fluctuateBins), nBins_(0), axisMin_(0.0), axisMax_(0.0), axisRange_(0.0) { // Constructor // Set the directory for the histogram hist_->SetDirectory(0); // Save various attributes of the histogram nBins_ = hist_->GetNbinsX(); TAxis* xAxis = hist_->GetXaxis(); axisMin_ = xAxis->GetXmin(); axisMax_ = xAxis->GetXmax(); axisRange_ = axisMax_ - axisMin_; // Check that axis range corresponds to range of abscissa if (TMath::Abs(this->getMinAbscissa() - axisMin_)>1e-6) { std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis minimum: " << axisMin_ << " does not correspond to abscissa minimum: " << this->getMinAbscissa() << "." << std::endl; gSystem->Exit(EXIT_FAILURE); } if (TMath::Abs(this->getMaxAbscissa() - axisMax_)>1e-6) { std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis maximum: " << axisMax_ << " does not correspond to abscissa maximum: " << this->getMaxAbscissa() << "." << std::endl; gSystem->Exit(EXIT_FAILURE); } // If the bins are to be fluctuated then do so now before // calculating anything that depends on the bin content. if (fluctuateBins) { this->doBinFluctuation(); } // Calculate the PDF normalisation. this->calcNorm(); // And check it is OK. this->checkNormalisation(); } Lau1DHistPdf::~Lau1DHistPdf() { // Destructor delete hist_; hist_ = 0; } void Lau1DHistPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ ) { if (this->heightUpToDate()) { return; } // Get the maximum height of the histogram Int_t maxBin = hist_->GetMaximumBin(); Double_t height = hist_->GetBinContent(maxBin); this->setMaxHeight(height); } void Lau1DHistPdf::calcNorm() { // Calculate the histogram normalisation. // Loop over the range to get the total area. // Just sum the contributions up using 1e-3 increments of the range. // Multiply the end result by dx. Double_t dx(1e-3*axisRange_); Double_t area(0.0); Double_t x(axisMin_ + dx/2.0); while (x > axisMin_ && x < axisMax_) { area += this->interpolate(x); x += dx; } Double_t norm = area*dx; this->setNorm(norm); } void Lau1DHistPdf::checkNormalisation() { Double_t dx(1e-3*axisRange_); Double_t area(0.0); Double_t areaNoNorm(0.0); const Bool_t useInterpolationOrig { useInterpolation_ }; Double_t x(axisMin_ + dx/2.0); while (x > axisMin_ && x < axisMax_) { area += this->interpolateNorm(x); useInterpolation_ = kFALSE; areaNoNorm += this->interpolate(x); useInterpolation_ = useInterpolationOrig; x += dx; } Double_t norm = area*dx; std::cout << "INFO in Lau1DHistPdf::checkNormalisation : Area = " << area << ", dx = " << dx << std::endl; std::cout << " : Area with no norm = " << areaNoNorm << "*dx = " << areaNoNorm*dx << std::endl; std::cout << " : The total area of the normalised histogram PDF is " << norm << std::endl; } Double_t Lau1DHistPdf::getBinHistValue(Int_t bin) const { // Check that bin is in range [1 , nBins_] if ((bin < 1) || (bin > nBins_)) { return 0.0; } Double_t value = static_cast(hist_->GetBinContent(bin)); // protect against negative values if ( value < 0.0 ) { - std::cerr << "WARNING in Lau1DHistPdf::getBinHistValue : Negative bin content set to zero!" << std::endl; + std::cerr << "WARNING in Lau1DHistPdf::getBinHistValue : Negative bin content set to zero!" + << (hist_->GetName() ? Form(" Culprit histogram: %s",hist_->GetName()) : "") << std::endl; value = 0.0; } return value; } Double_t Lau1DHistPdf::interpolateNorm(Double_t x) const { // Get the normalised interpolated value. Double_t value = this->interpolate(x); Double_t norm = this->getNorm(); return value/norm; } void Lau1DHistPdf::calcLikelihoodInfo(const LauAbscissas& abscissas) { // Check that the given abscissa is within the allowed range if (!this->checkRange(abscissas)) { gSystem->Exit(EXIT_FAILURE); } // Get our abscissa Double_t abscissa = abscissas[0]; // Calculate the interpolated value Double_t value = this->interpolate(abscissa); this->setUnNormPDFVal(value); } Double_t Lau1DHistPdf::interpolate(Double_t x) const { // This function returns the interpolated value of the histogram function // for the given value of x by finding the adjacent bins and extrapolating // using weights based on the inverse distance of the point from the adajcent // bin centres. // Find the histogram bin Int_t bin = hist_->FindFixBin(x); // Ask whether we want to do the interpolation, since this method is // not reliable for low statistics histograms. if (useInterpolation_ == kFALSE) { return this->getBinHistValue(bin); } // Find the bin centres (actual co-ordinate positions, not histogram indices) Double_t cbinx = hist_->GetBinCenter(bin); // Find the adjacent bins Double_t deltax = x - cbinx; Int_t bin_adj(0); if (deltax > 0.0) { bin_adj = bin + 1; } else { bin_adj = bin - 1; } Bool_t isBoundary(kFALSE); if ( bin_adj > nBins_ || bin_adj < 1 ) { isBoundary = kTRUE; } // At the edges, do no interpolation, use entry in bin. if (isBoundary == kTRUE) { return this->getBinHistValue(bin); } // Linear interpolation using inverse distance as weights. // Find the adjacent bin centre Double_t cbinx_adj = hist_->GetBinCenter(bin_adj); Double_t deltax_adj = cbinx_adj - x; Double_t dx0 = TMath::Abs(deltax); Double_t dx1 = TMath::Abs(deltax_adj); Double_t denom = dx0 + dx1; Double_t value0 = this->getBinHistValue(bin); Double_t value1 = this->getBinHistValue(bin_adj); Double_t value = value0*dx1 + value1*dx0; value /= denom; return value; } void Lau1DHistPdf::doBinFluctuation() { TRandom* random = LauRandom::randomFun(); for (Int_t bin(0); binGetBinContent(bin+1); Double_t currentError = hist_->GetBinError(bin+1); Double_t newContent = random->Gaus(currentContent,currentError); if (newContent<0.0) { hist_->SetBinContent(bin+1,0.0); } else { hist_->SetBinContent(bin+1,newContent); } } } diff --git a/src/LauFlavTag.cc b/src/LauFlavTag.cc index de1dac4..91f6e99 100644 --- a/src/LauFlavTag.cc +++ b/src/LauFlavTag.cc @@ -1,1106 +1,1107 @@ /* 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) LauFlavTag::LauFlavTag(const Bool_t useAveDelta, const Bool_t useEtaPrime, const std::map bkgndInfo) : useAveDelta_(useAveDelta), useEtaPrime_(useEtaPrime) { // Put map values into vectors for (const auto &bkgnd : bkgndInfo){ bkgndNames_.push_back(bkgnd.first); bkgndTypes_.push_back(bkgnd.second); std::cout<< "INFO in LauFlavTag::LauFlavTag : adding background " << bkgnd.first << " of type " << bkgnd.second < 0){ if (!useAveDelta_){ tagEffBkgnd_B0_.clear(); tagEffBkgnd_B0_.resize(nBkgnds); tagEffBkgnd_B0bar_.clear(); tagEffBkgnd_B0bar_.resize(nBkgnds); tagEffBkgnd_hist_B0_.clear(); tagEffBkgnd_hist_B0_.resize(nBkgnds); tagEffBkgnd_hist_B0bar_.clear(); tagEffBkgnd_hist_B0bar_.resize(nBkgnds); } else { tagEffBkgnd_ave_.clear(); tagEffBkgnd_ave_.resize(nBkgnds); tagEffBkgnd_delta_.clear(); tagEffBkgnd_delta_.resize(nBkgnds); tagEffBkgnd_hist_ave_.clear(); tagEffBkgnd_hist_ave_.resize(nBkgnds); tagEffBkgnd_hist_delta_.clear(); tagEffBkgnd_hist_delta_.resize(nBkgnds); } etaBkgndPdfs_.clear(); etaBkgndPdfs_.resize(nBkgnds); } } 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 ( taggerPosition_.find(name) != taggerPosition_.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 ULong_t position { tagVarNames_.size() }; // Update map to relate tagger name and position in the vectors taggerPosition_[name] = position; // 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,position,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_[position]->initValue(tagEff.first); tagEff_B0_[position]->genValue(tagEff.first); tagEff_B0_[position]->fixed(kTRUE); if (tagEff.second==-1.0){ tagEff_B0bar_.push_back(tagEff_B0_[position]->createClone(tagEff_b0barName)); } else { LauParameter* tageffb0bar = new LauParameter(tagEff_b0barName,tagEff.second,0.0,1.0,kTRUE); tagEff_B0bar_.push_back(tageffb0bar); tagEff_B0bar_[position]->initValue(tagEff.second); tagEff_B0bar_[position]->genValue(tagEff.second); tagEff_B0bar_[position]->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_[position]->initValue(tagEff.first); tagEff_ave_[position]->genValue(tagEff.first); tagEff_ave_[position]->fixed(kTRUE); LauParameter* tageffdelta = new LauParameter(tagEff_deltaName,tagEff.second,-1.0,1.0,kTRUE); tagEff_delta_.push_back(tageffdelta); tagEff_delta_[position]->initValue(tagEff.second); tagEff_delta_[position]->genValue(tagEff.second); tagEff_delta_[position]->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 ( taggerPosition_.find(name) != taggerPosition_.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 ULong_t position { tagVarNames_.size() }; // Update map to relate tagger name and position in the vectors taggerPosition_[name] = position; // 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,position,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 (bkgndNames_.size()>0){ if (!useAveDelta_){ //Loop over the outer vector and extend the inner vectors for (auto& innerVec : tagEffBkgnd_B0_){ innerVec.push_back(nullptr); } for (auto& innerVec : tagEffBkgnd_hist_B0_){ innerVec.push_back(nullptr); } for (auto& innerVec : tagEffBkgnd_B0bar_){ innerVec.push_back(nullptr); } for (auto& innerVec : tagEffBkgnd_hist_B0bar_){ innerVec.push_back(nullptr); } } else { //Loop over the outer vector and extend the inner vectors for (auto& innerVec : tagEffBkgnd_ave_){ innerVec.push_back(nullptr); } for (auto& innerVec : tagEffBkgnd_hist_ave_){ innerVec.push_back(nullptr); } for (auto& innerVec : tagEffBkgnd_delta_){ innerVec.push_back(nullptr); } for (auto& innerVec : tagEffBkgnd_hist_delta_){ innerVec.push_back(nullptr); } } for (auto& innerVec : etaBkgndPdfs_){ innerVec.push_back(nullptr); } } } void LauFlavTag::setupCalibParams(const TString& name, const ULong_t position, const std::pair calib_p0, const std::pair calib_p1) { if (!useAveDelta_){ TString calib_p0_b0Name("calib_p0_b0_"+name); TString calib_p0_b0barName("calib_p0_b0bar_"+name); TString calib_p1_b0Name("calib_p1_b0_"+name); TString calib_p1_b0barName("calib_p1_b0bar_"+name); LauParameter* calibp0b0 = new LauParameter(calib_p0_b0Name,calib_p0.first,-10.0,10.0,kTRUE); calib_p0_B0_.push_back(calibp0b0); calib_p0_B0_[position]->initValue(calib_p0.first); calib_p0_B0_[position]->genValue(calib_p0.first); calib_p0_B0_[position]->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_[position]->initValue(calib_p1.first); calib_p1_B0_[position]->genValue(calib_p1.first); calib_p1_B0_[position]->fixed(kTRUE); if (calib_p0.second==-1.0 && calib_p1.second==-1.0){ calib_p0_B0bar_.push_back(calib_p0_B0_[position]->createClone(calib_p0_b0barName)); calib_p1_B0bar_.push_back(calib_p1_B0_[position]->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_[position]->initValue(calib_p0.second); calib_p0_B0bar_[position]->genValue(calib_p0.second); calib_p0_B0bar_[position]->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_[position]->initValue(calib_p1.second); calib_p1_B0bar_[position]->genValue(calib_p1.second); calib_p1_B0bar_[position]->fixed(kTRUE); } } else { //Use average and delta variables TString calib_p0_aveName("calib_p0_ave_"+name); TString calib_p0_deltaName("calib_p0_delta_"+name); TString calib_p1_aveName("calib_p1_ave_"+name); TString calib_p1_deltaName("calib_p1_delta_"+name); LauParameter* calibp0ave = new LauParameter(calib_p0_aveName,calib_p0.first,-10.0,10.0,kTRUE); calib_p0_ave_.push_back(calibp0ave); calib_p0_ave_[position]->initValue(calib_p0.first); calib_p0_ave_[position]->genValue(calib_p0.first); calib_p0_ave_[position]->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_[position]->initValue(calib_p1.first); calib_p1_ave_[position]->genValue(calib_p1.first); calib_p1_ave_[position]->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_[position]->initValue(calib_p0.second); calib_p0_delta_[position]->genValue(calib_p0.second); calib_p0_delta_[position]->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_[position]->initValue(calib_p1.second); calib_p1_delta_[position]->genValue(calib_p1.second); calib_p1_delta_[position]->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 (ULong_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 ULong_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 (ULong_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 (ULong_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 ULong_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 ULong_t nTaggers { this->getNTaggers() }; for ( ULong_t position{0}; positiongetEtaGen(position); if (this->getUseAveDelta()) { if (tagEff_ave_[position]==nullptr){ const Double_t ave = tagEff_hist_ave_[position]->GetBinContent(tagEff_hist_ave_[position]->FindFixBin(curEvtDecayTime_)); const Double_t delta = tagEff_hist_delta_[position]->GetBinContent(tagEff_hist_delta_[position]->FindFixBin(curEvtDecayTime_)); tagEffB0 = ave + 0.5 * delta; tagEffB0bar = ave - 0.5 * delta; } else { tagEffB0 = tagEff_ave_[position]->unblindValue() + 0.5 * tagEff_delta_[position]->unblindValue(); tagEffB0bar = tagEff_ave_[position]->unblindValue() - 0.5 * tagEff_delta_[position]->unblindValue(); } } else { if (tagEff_B0_[position]==nullptr){ tagEffB0 = tagEff_hist_B0_[position]->GetBinContent(tagEff_hist_B0_[position]->FindFixBin(curEvtDecayTime_)); tagEffB0bar = tagEff_hist_B0bar_[position]->GetBinContent(tagEff_hist_B0bar_[position]->FindFixBin(curEvtDecayTime_)); } else { tagEffB0 = tagEff_B0_[position]->unblindValue(); tagEffB0bar = tagEff_B0bar_[position]->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(position, Flavour::B)){ curEvtTagFlv_[position] = Flavour::B; } else { curEvtTagFlv_[position] = Flavour::Bbar; } } else { curEvtTagFlv_[position] = 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(position, Flavour::Bbar)){ curEvtTagFlv_[position] = Flavour::Bbar; } else { curEvtTagFlv_[position] = Flavour::B; } } else { curEvtTagFlv_[position] = 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 Flavour trueTagFlv, const Double_t curEvtDecayTime, const ULong_t bkgndID) { if (bkgndID < 0 || bkgndID > bkgndNames_.size()){ std::cerr << "ERROR in LauFlavTag::generateBkgndEventInfo : Invalid backgrond class identifier" << std::endl; gSystem->Exit(EXIT_FAILURE); } curEvtTrueTagFlv_ = trueTagFlv; curEvtDecayFlv_ = Flavour::Unknown; curEvtDecayTime_ = curEvtDecayTime; Double_t randNo{0.0}; Double_t tagEffB0{0.0}; Double_t tagEffB0bar{0.0}; const ULong_t nTaggers { this->getNTaggers() }; for ( ULong_t position{0}; positiongetEtaGenBkgnd(position,bkgndID); //TODO If bkgnd is signal like should these parameters be clones of signal TagEff etc? //TODO Or call generateEventInfo() instead? if (this->getUseAveDelta()) { if (tagEffBkgnd_ave_[bkgndID][position]==nullptr){ const Double_t ave = tagEffBkgnd_hist_ave_[bkgndID][position]->GetBinContent(tagEffBkgnd_hist_ave_[bkgndID][position]->FindFixBin(curEvtDecayTime_)); const Double_t delta = tagEffBkgnd_hist_delta_[bkgndID][position]->GetBinContent(tagEffBkgnd_hist_delta_[bkgndID][position]->FindFixBin(curEvtDecayTime_)); tagEffB0 = ave + 0.5 * delta; tagEffB0bar = ave - 0.5 * delta; } else { tagEffB0 = tagEffBkgnd_ave_[bkgndID][position]->unblindValue() + 0.5 * tagEffBkgnd_delta_[bkgndID][position]->unblindValue(); tagEffB0bar = tagEffBkgnd_ave_[bkgndID][position]->unblindValue() - 0.5 * tagEffBkgnd_delta_[bkgndID][position]->unblindValue(); } } else { if (tagEffBkgnd_B0_[bkgndID][position]==nullptr){ tagEffB0 = tagEffBkgnd_hist_B0_[bkgndID][position]->GetBinContent(tagEffBkgnd_hist_B0_[bkgndID][position]->FindFixBin(curEvtDecayTime_)); tagEffB0bar = tagEffBkgnd_hist_B0bar_[bkgndID][position]->GetBinContent(tagEffBkgnd_hist_B0bar_[bkgndID][position]->FindFixBin(curEvtDecayTime_)); } else { tagEffB0 = tagEffBkgnd_B0_[bkgndID][position]->unblindValue(); tagEffB0bar = tagEffBkgnd_B0bar_[bkgndID][position]->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 mistag - use eta not littleOmega for now (littleOmega only for SignalLike bkgnd?) //if (randNo > this->getLittleOmega(position, Flavour::B)){ if (randNo > curEvtMistag_[position]){ curEvtTagFlv_[position] = Flavour::B; } else { curEvtTagFlv_[position] = Flavour::Bbar; } } else { curEvtTagFlv_[position] = 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(position, Flavour::Bbar)){ if (randNo > curEvtMistag_[position]){ curEvtTagFlv_[position] = Flavour::Bbar; } else { curEvtTagFlv_[position] = Flavour::B; } } else { curEvtTagFlv_[position] = 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 ULong_t position, 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_[position]->unblindValue() + 0.5*calib_p0_delta_[position]->unblindValue(); calibp0bar = calib_p0_ave_[position]->unblindValue() - 0.5*calib_p0_delta_[position]->unblindValue(); calibp1 = calib_p1_ave_[position]->unblindValue() + 0.5*calib_p1_delta_[position]->unblindValue(); calibp1bar = calib_p1_ave_[position]->unblindValue() - 0.5*calib_p1_delta_[position]->unblindValue(); } else { calibp0 = calib_p0_B0_[position]->unblindValue(); calibp0bar = calib_p0_B0bar_[position]->unblindValue(); calibp1 = calib_p1_B0_[position]->unblindValue(); calibp1bar = calib_p1_B0bar_[position]->unblindValue(); } if ( flag == Flavour::B ){ return calibp0 + calibp1 * (curEvtMistag_[position] - perEvtAvgMistag_[position]); } else{ return calibp0bar + calibp1bar * (curEvtMistag_[position] - perEvtAvgMistag_[position]); } return 0.0; } Double_t LauFlavTag::getCapitalOmega(const ULong_t position, 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 Int_t deltap1(0), deltam1(0), delta0(0); if (curEvtTagFlv_[position] == Flavour::Bbar){ deltam1 = 1; } else if(curEvtTagFlv_[position] == Flavour::B){ deltap1 = 1; } else{ delta0 = 1; } //Efficiency Double_t eff(0.0); if (useAveDelta_){ if(flag==Flavour::B){ if (tagEff_ave_[position]==nullptr){ eff = tagEff_hist_ave_[position]->GetBinContent(tagEff_hist_ave_[position]->FindFixBin(curEvtDecayTime_)) + 0.5*tagEff_hist_delta_[position]->GetBinContent(tagEff_hist_delta_[position]->FindFixBin(curEvtDecayTime_)); } else { eff = tagEff_ave_[position]->unblindValue() + 0.5*tagEff_delta_[position]->unblindValue(); } } else { if (tagEff_ave_[position]==nullptr){ eff = tagEff_hist_ave_[position]->GetBinContent(tagEff_hist_ave_[position]->FindFixBin(curEvtDecayTime_)) - 0.5*tagEff_hist_delta_[position]->GetBinContent(tagEff_hist_delta_[position]->FindFixBin(curEvtDecayTime_)); } else { eff = tagEff_ave_[position]->unblindValue() - 0.5*tagEff_delta_[position]->unblindValue(); } } }else{ if(flag==Flavour::B){ if (tagEff_B0_[position]==nullptr){ eff = tagEff_hist_B0_[position]->GetBinContent(tagEff_hist_B0_[position]->FindFixBin(curEvtDecayTime_)); } else { eff = tagEff_B0_[position]->unblindValue(); } }else{ if (tagEff_B0bar_[position]==nullptr){ eff = tagEff_hist_B0bar_[position]->GetBinContent(tagEff_hist_B0bar_[position]->FindFixBin(curEvtDecayTime_)); } else { eff = tagEff_B0bar_[position]->unblindValue(); } } } //Little omega Double_t omega = this->getLittleOmega(position, flag); Double_t omegaPrime(0.); //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? if (useEtaPrime_){ omegaPrime = (1/(1+TMath::Exp(-1.0*omega))); }else{ omegaPrime = 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. if (omegaPrime < 0.0){ std::cerr << "WARNING in LauFlavTag::getCapitalOmega : The value of little omega is less than 0, shifting to 0" << std::endl; omegaPrime = 0.0; } if (omegaPrime > 1.0){ std::cerr << "WARNING in LauFlavTag::getCapitalOmega : The value of little omega is greater than 1, shifting to 1" << std::endl; omegaPrime = 1.0; } //eta PDF value std::vector abs; abs.push_back(curEvtMistag_[position]); etaPdfs_[position]->calcLikelihoodInfo(abs); Double_t h { etaPdfs_[position]->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 if (h==0 && delta0==0){ std::cerr << "WARNING in LauFlavTag::getCapitalOmega : The value of the eta PDF is zero at eta = " << curEvtMistag_[position] << ", shifting to 0.1" << std::endl; h=0.1; } //Put it together if (flag == Flavour::B){ return (deltap1*eff*(1-omegaPrime) + deltam1*eff*omegaPrime)*h + delta0*(1-eff)*u; } else { return (deltam1*eff*(1-omegaPrime) + deltap1*eff*omegaPrime)*h + delta0*(1-eff)*u; } } Double_t LauFlavTag::getCapitalOmegaBkgnd(const ULong_t position, const Flavour flag, const UInt_t classID) const { //Fill in with the various options of flag = +-1, type = signal-like, combinatorial etc 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; } //Delta functions to control which terms contribute Int_t deltap1(0), deltam1(0), delta0(0); if (curEvtTagFlv_[position] == Flavour::Bbar){ deltam1 = 1; } else if(curEvtTagFlv_[position] == Flavour::B){ deltap1 = 1; } else{ delta0 = 1; } //Efficiency Double_t effB0(0.0), effB0bar(0.0); if (useAveDelta_){ if(flag==Flavour::B){ if (tagEffBkgnd_ave_[classID][position]==nullptr){ effB0 = tagEffBkgnd_hist_ave_[classID][position]->GetBinContent(tagEffBkgnd_hist_ave_[classID][position]->FindFixBin(curEvtDecayTime_)) + 0.5*tagEffBkgnd_hist_delta_[classID][position]->GetBinContent(tagEffBkgnd_hist_delta_[classID][position]->FindFixBin(curEvtDecayTime_)); } else { effB0 = tagEffBkgnd_ave_[classID][position]->unblindValue() + 0.5*tagEffBkgnd_delta_[classID][position]->unblindValue(); } } else { if (tagEffBkgnd_ave_[classID][position]==nullptr){ effB0bar = tagEffBkgnd_hist_ave_[classID][position]->GetBinContent(tagEffBkgnd_hist_ave_[classID][position]->FindFixBin(curEvtDecayTime_)) - 0.5*tagEffBkgnd_hist_delta_[classID][position]->GetBinContent(tagEffBkgnd_hist_delta_[classID][position]->FindFixBin(curEvtDecayTime_)); } else { effB0bar = tagEffBkgnd_ave_[classID][position]->unblindValue() - 0.5*tagEffBkgnd_delta_[classID][position]->unblindValue(); } } }else{ if(flag==Flavour::B){ if (tagEffBkgnd_B0_[classID][position]==nullptr){ effB0 = tagEffBkgnd_hist_B0_[classID][position]->GetBinContent(tagEffBkgnd_hist_B0_[classID][position]->FindFixBin(curEvtDecayTime_)); } else { effB0 = tagEffBkgnd_B0_[classID][position]->unblindValue(); } }else{ if (tagEffBkgnd_B0bar_[classID][position]==nullptr){ effB0bar = tagEffBkgnd_hist_B0bar_[classID][position]->GetBinContent(tagEffBkgnd_hist_B0bar_[classID][position]->FindFixBin(curEvtDecayTime_)); } else { effB0bar = tagEffBkgnd_B0bar_[classID][position]->unblindValue(); } } } //Need to know which background eta PDF to use - classID std::vector abs; abs.push_back(curEvtMistag_[position]); etaBkgndPdfs_[classID][position]->calcLikelihoodInfo(abs); Double_t h { etaBkgndPdfs_[classID][position]->getLikelihood() }; const Double_t u { 2.0 }; // the PDF value for a uniform PDF between 0.0 and 0.5 if (bkgndTypes_[classID] == BkgndType::Combinatorial){ return (deltap1*effB0 + deltam1*effB0bar)*h + delta0*(1.0 - 0.5*(effB0 + effB0bar))*u; } else { return 1.0; } } void LauFlavTag::setBkgndParams(const TString& bkgndName, const TString& taggerName, LauAbsPdf* etaPdf, std::pair tagEff) { if (taggerPosition_.count(taggerName)==0){ std::cerr << "ERROR in LauFlavTag::setBkgndParams : Tagger name not recognised please check your options" << std::endl; return; } Int_t bkgndID(-1); for (ULong_t i=0; iname("tagEff_ave_"+taggerName+"_bkgnd_"+bkgndName); tagEffBkgnd_ave_[bkgndID][position] = tagEffB0; tagEffB0bar->name("tagEff_delta_"+taggerName+"_bkgnd_"+bkgndName); tagEffB0bar->range(-1.0,1.0); tagEffBkgnd_delta_[bkgndID][position] = tagEffB0bar; } else { tagEffBkgnd_B0_[bkgndID][position] = tagEffB0; tagEffBkgnd_B0bar_[bkgndID][position] = tagEffB0bar; } 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 (taggerPosition_.count(taggerName)==0){ std::cerr << "ERROR in LauFlavTag::setBkgndParams : Tagger name 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; } Int_t bkgndID(-1); for (ULong_t i=0; igenerate(nullptr) }; Double_t etagen { data.at(etaPdfs_[position]->varName()) }; if (etagen > 0.5){etagen = 0.5;} if (etagen < 0.0){etagen = 0.0;} curEvtMistag_[position] = etagen; return etagen; } Double_t LauFlavTag::getEtaGenBkgnd(const ULong_t position, const ULong_t bkgndID) { LauFitData data { etaBkgndPdfs_[bkgndID][position]->generate(nullptr) }; Double_t etagen { data.at(etaBkgndPdfs_[bkgndID][position]->varName()) }; if (etagen > 0.5){etagen = 0.5;} if (etagen < 0.0){etagen = 0.0;} curEvtMistag_[position] = etagen; return etagen; } void LauFlavTag::setTrueTagVarName(TString trueTagVarName){ trueTagVarName_ = std::move(trueTagVarName); } void LauFlavTag::setDecayFlvVarName(TString decayFlvVarName){ decayFlvVarName_ = std::move(decayFlvVarName); } void LauFlavTag::addP0GaussianConstraints(TString name, std::pair constraint1, std::pair constraint2){ //Does key exist? if (taggerPosition_.count(name)==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 position in the vector from the tagger name Double_t pos = taggerPosition_.at(name); 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 " << name << std::endl; } void LauFlavTag::addP1GaussianConstraints(TString name, std::pair constraint1, std::pair constraint2){ //Does key exist? if (taggerPosition_.count(name)==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 position in the vector from the tagger name Double_t pos = taggerPosition_.at(name); 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 " << name << std::endl; } void LauFlavTag::addTagEffGaussianConstraints(TString name, std::pair constraint1, std::pair constraint2){ //Does key exist? if (taggerPosition_.count(name)==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 position in the vector from the tagger name Double_t pos = taggerPosition_.at(name); 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 " << name << std::endl; } void LauFlavTag::floatCalibParP0B0(const TString name){ if (useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0B0 : Trying to set calibration parameters for B0/B0bar not average/delta" << std::endl; return; } if (name==""){ for (auto& param : calib_p0_B0_){ if (param==nullptr){continue;} param->fixed(kFALSE); } } else { //Does key exist? if (taggerPosition_.count(name)==0){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0B0 : Tagger name not recognised please check your options" << std::endl; return; } //Find position in the vector from the tagger name Double_t position = taggerPosition_.at(name); calib_p0_B0_[position]->fixed(kFALSE); } } void LauFlavTag::floatCalibParP1B0(const TString name){ if (useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1B0 : Trying to set calibration parameters for B0/B0bar not average/delta" << std::endl; return; } if (name==""){ for (auto& param : calib_p1_B0_){ if (param==nullptr){continue;} param->fixed(kFALSE); } } else { //Does key exist? if (taggerPosition_.count(name)==0){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1B0 : Tagger name not recognised please check your options" << std::endl; return; } //Find position in the vector from the tagger name Double_t position = taggerPosition_.at(name); calib_p1_B0_[position]->fixed(kFALSE); } } void LauFlavTag::floatCalibParP0B0bar(const TString name){ if (useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0B0bar : Trying to set calibration parameters for B0/B0bar not average/delta" << std::endl; return; } if (name==""){ for (auto& param : calib_p0_B0bar_){ if (param==nullptr){continue;} param->fixed(kFALSE); } } else { //Does key exist? if (taggerPosition_.count(name)==0){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0B0bar : Tagger name not recognised please check your options" << std::endl; return; } //Find position in the vector from the tagger name Double_t position = taggerPosition_.at(name); calib_p0_B0bar_[position]->fixed(kFALSE); } } void LauFlavTag::floatCalibParP1B0bar(const TString name){ if (useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1B0bar : Trying to set calibration parameters for B0/B0bar not average/delta" << std::endl; return; } if (name==""){ for (auto& param : calib_p1_B0bar_){ if (param==nullptr){continue;} param->fixed(kFALSE); } } else { //Does key exist? if (taggerPosition_.count(name)==0){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1B0bar : Tagger name not recognised please check your options" << std::endl; return; } //Find position in the vector from the tagger name Double_t position = taggerPosition_.at(name); calib_p1_B0bar_[position]->fixed(kFALSE); } } void LauFlavTag::floatCalibParP0Ave(const TString name){ if (!useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0Ave : Trying to set calibration parameters for average/delta not B0/B0bar" << std::endl; return; } if (name==""){ for (auto& param : calib_p0_ave_){ if (param==nullptr){continue;} param->fixed(kFALSE); } } else { //Does key exist? if (taggerPosition_.count(name)==0){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0Ave : Tagger name not recognised please check your options" << std::endl; return; } //Find position in the vector from the tagger name Double_t position = taggerPosition_.at(name); calib_p0_ave_[position]->fixed(kFALSE); } } void LauFlavTag::floatCalibParP0Delta(const TString name){ if (!useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0Delta : Trying to set calibration parameters for average/delta not B0/B0bar" << std::endl; return; } if (name==""){ for (auto& param : calib_p0_delta_){ if (param==nullptr){continue;} param->fixed(kFALSE); } } else { //Does key exist? if (taggerPosition_.count(name)==0){ std::cerr << "ERROR in LauFlavTag::floatCalibParP0Delta : Tagger name not recognised please check your options" << std::endl; return; } //Find position in the vector from the tagger name Double_t position = taggerPosition_.at(name); calib_p0_delta_[position]->fixed(kFALSE); } } void LauFlavTag::floatCalibParP1Ave(const TString name){ if (!useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1Ave : Trying to set calibration parameters for average/delta not B0/B0bar" << std::endl; return; } if (name==""){ for (auto& param : calib_p1_ave_){ if (param==nullptr){continue;} param->fixed(kFALSE); } } else { //Does key exist? if (taggerPosition_.count(name)==0){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1Ave : Tagger name not recognised please check your options" << std::endl; return; } //Find position in the vector from the tagger name Double_t position = taggerPosition_.at(name); calib_p1_ave_[position]->fixed(kFALSE); } } void LauFlavTag::floatCalibParP1Delta(const TString name){ if (!useAveDelta_){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1Delta : Trying to set calibration parameters for average/delta not B0/B0bar" << std::endl; return; } if (name==""){ for (auto& param : calib_p1_delta_){ if (param==nullptr){continue;} param->fixed(kFALSE); } } else { //Does key exist? if (taggerPosition_.count(name)==0){ std::cerr << "ERROR in LauFlavTag::floatCalibParP1Delta : Tagger name not recognised please check your options" << std::endl; return; } //Find position in the vector from the tagger name Double_t position = taggerPosition_.at(name); calib_p1_delta_[position]->fixed(kFALSE); } } void LauFlavTag::floatAllCalibPars(){ if (useAveDelta_){ this->floatCalibParP0Ave(); this->floatCalibParP0Delta(); this->floatCalibParP1Ave(); this->floatCalibParP1Delta(); } else { this->floatCalibParP0B0(); this->floatCalibParP1B0(); this->floatCalibParP0B0bar(); this->floatCalibParP1B0bar(); } }