diff --git a/inc/Lau1DHistPdf.hh b/inc/Lau1DHistPdf.hh index e075124..488f13d 100644 --- a/inc/Lau1DHistPdf.hh +++ b/inc/Lau1DHistPdf.hh @@ -1,136 +1,139 @@ /* 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.hh \brief File containing declaration of Lau1DHistPdf class. */ /*! \class Lau1DHistPdf \brief Class for defining a 1D histogram PDF. Class for defining a 1D histogram PDF. Employs linear interpolation to get the PDF value based on how far away a point is from nearby bin centres. The returned values are normalised to the total area. */ #ifndef LAU_1DHIST_PDF #define LAU_1DHIST_PDF +#include "TH1.h" #include "LauAbsPdf.hh" class TH1; class Lau1DHistPdf : public LauAbsPdf { public: //! Constructor /*! \param [in] theVarName the name of the abscissa variable \param [in] hist the 1D histogram from which the PDF should be constructed \param [in] minAbscissa the minimum value of the abscissa \param [in] maxAbscissa the maximum value of the abscissa \param [in] useInterpolation boolean flag to determine whether linear interpolation between bins should be used or simply the raw bin values \param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors. The seed for the random number generator used to fluctuate the bins should first be set using LauRandom::setSeed. */ Lau1DHistPdf(const TString& theVarName, const TH1* hist, Double_t minAbscissa, Double_t maxAbscissa, Bool_t useInterpolation = kTRUE, Bool_t fluctuateBins = kFALSE); //! Destructor virtual ~Lau1DHistPdf(); //! Calculate the likelihood (and intermediate info) for a given abscissa /*! \param [in] abscissas the values of the abscissa(s) */ virtual void calcLikelihoodInfo(const LauAbscissas& abscissas); using LauAbsPdf::calcLikelihoodInfo; //! Calculate the normalisation virtual void calcNorm(); //! Calculate the PDF height /*! \param [in] kinematics the current DP kinematics (not strictly required in this case since PDF has no DP dependence) */ virtual void calcPDFHeight( const LauKinematics* kinematics ); + Double_t getMean(){return hist_->GetMean();}; + protected: //! Fluctuate the histogram bin contents in accorance with their errors void doBinFluctuation(); //! Check the normalisation calculation void checkNormalisation(); //! Get the bin content from the histogram /*! \param [in] bin the bin number \return the bin content */ Double_t getBinHistValue(Int_t bin) const; //! Perform the interpolation (unnormalised) /*! \param [in] x the abscissa value \return the unnormalised PDF value */ Double_t interpolate(Double_t x) const; //! Perform the interpolation and divide by the normalisation /*! \param [in] x the abscissa value \return the normalised PDF value */ Double_t interpolateNorm(Double_t x) const; private: //! Copy constructor (not implemented) Lau1DHistPdf(const Lau1DHistPdf& rhs); //! Copy assignment operator (not implemented) Lau1DHistPdf& operator=(const Lau1DHistPdf& rhs); //! The underlying histogram TH1* hist_; //! Control boolean for using the linear interpolation Bool_t useInterpolation_; //! Control boolean for performing the fluctuation of the histogram bin contents Bool_t fluctuateBins_; //! The number of bins in the histogram Int_t nBins_; //! The histogram axis minimum Double_t axisMin_; //! The histogram axis maximum Double_t axisMax_; //! The histogram axis range Double_t axisRange_; ClassDef(Lau1DHistPdf,0) // 1D histogram pdf class }; #endif diff --git a/src/LauFlavTag.cc b/src/LauFlavTag.cc index c03ed62..78fd7ec 100644 --- a/src/LauFlavTag.cc +++ b/src/LauFlavTag.cc @@ -1,454 +1,453 @@ /* 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 #include #include "TFile.h" #include "TMinuit.h" #include "TRandom.h" #include "TSystem.h" #include "TVirtualFitter.h" #include "TH1D.h" #include "LauAbsBkgndDPModel.hh" #include "LauAbsCoeffSet.hh" #include "LauAbsPdf.hh" #include "LauAsymmCalc.hh" #include "LauComplex.hh" #include "LauConstants.hh" #include "LauDPPartialIntegralInfo.hh" #include "LauDaughters.hh" #include "LauDecayTimePdf.hh" #include "LauFitNtuple.hh" #include "LauGenNtuple.hh" #include "LauIsobarDynamics.hh" #include "LauKinematics.hh" #include "LauPrint.hh" #include "LauRandom.hh" #include "LauScfMap.hh" #include "LauFlavTag.hh" #include "Lau1DHistPdf.hh" ClassImp(LauFlavTag) LauFlavTag::LauFlavTag() : taggerPosition_(), tagVarName_(), mistagVarName_(), curEvtTagFlv_(), curEvtMistag_(), curEvtTrueTagFlv_(), perEvtAvgMistag_(), calib_p0_B0_(), calib_p0_B0bar_(), calib_p1_B0_(), calib_p1_B0bar_(), calib_p0_ave_(), calib_p0_delta_(), calib_p1_ave_(), calib_p1_delta_(), useAveDelta_(kFALSE), etaPdfs_() { } LauFlavTag::~LauFlavTag() { } void LauFlavTag::initialise() { } void LauFlavTag::addTagger(const TString name, const TString tagVarName, const TString mistagVarName, LauAbsPdf* etapdf, const Double_t tagEff_b0, const Double_t calib_p0_b0, const Double_t calib_p1_b0, const Double_t tagEff_b0bar, const Double_t calib_p0_b0bar, const Double_t calib_p1_b0bar) { // Find how many taggers have already been added Int_t position = tagVarName_.size(); // Update map to relate tagger name and position in the vectors taggerPosition_[position]=name; // Fill vectors tagVarName_.push_back(tagVarName); mistagVarName_.push_back(mistagVarName); etaPdfs_.push_back(etapdf); - //TODO UPDATE THIS - TEMP THING TO TEST THE REST - perEvtAvgMistag_.push_back(0.4); + perEvtAvgMistag_.push_back(static_cast(etapdf)->getMean()); // LauParameters TString tagEff_b0Name("tagEff_b0_"+name); TString tagEff_b0barName("tagEff_b0bar_"+name); 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* tageffb0 = new LauParameter(tagEff_b0Name,tagEff_b0,0.0,1.0,kTRUE); tagEff_B0_.push_back(tageffb0); tagEff_B0_[position]->initValue(tagEff_b0); tagEff_B0_[position]->genValue(tagEff_b0); tagEff_B0_[position]->fixed(kTRUE); //Update once full code in place LauParameter* calibp0b0 = new LauParameter(calib_p0_b0Name,calib_p0_b0,0.0,1.0,kTRUE); calib_p0_B0_.push_back(calibp0b0); calib_p0_B0_[position]->initValue(calib_p0_b0); calib_p0_B0_[position]->genValue(calib_p0_b0); calib_p0_B0_[position]->fixed(kTRUE); //Update once full code in place LauParameter* calibp1b0 = new LauParameter(calib_p1_b0Name,calib_p1_b0,0.0,1.5,kTRUE); calib_p1_B0_.push_back(calibp1b0); calib_p1_B0_[position]->initValue(calib_p1_b0); calib_p1_B0_[position]->genValue(calib_p1_b0); calib_p1_B0_[position]->fixed(kTRUE); //Update once full code in place if (tagEff_b0bar==-1.0 && calib_p0_b0bar==-1.0 && calib_p1_b0bar==-1.0){ tagEff_B0bar_.push_back(tagEff_B0_[position]->createClone(tagEff_b0barName)); 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* tageffb0bar = new LauParameter(tagEff_b0barName,tagEff_b0bar,0.0,1.0,kTRUE); tagEff_B0bar_.push_back(tageffb0bar); tagEff_B0bar_[position]->initValue(tagEff_b0bar); tagEff_B0bar_[position]->genValue(tagEff_b0bar); tagEff_B0bar_[position]->fixed(kTRUE); //Update once full code in place LauParameter* calibp0b0bar = new LauParameter(calib_p0_b0barName,calib_p0_b0bar,0.0,1.0,kTRUE); calib_p0_B0bar_.push_back(calibp0b0bar); calib_p0_B0bar_[position]->initValue(calib_p0_b0bar); calib_p0_B0bar_[position]->genValue(calib_p0_b0bar); calib_p0_B0bar_[position]->fixed(kTRUE); //Update once full code in place LauParameter* calibp1b0bar = new LauParameter(calib_p1_b0barName,calib_p1_b0bar,0.0,1.5,kTRUE); calib_p1_B0bar_.push_back(calibp1b0bar); calib_p1_B0bar_[position]->initValue(calib_p1_b0bar); calib_p1_B0bar_[position]->genValue(calib_p1_b0bar); calib_p1_B0bar_[position]->fixed(kTRUE); //Update once full code in place } // Use average and delta variables if required if (useAveDelta_){ this->useAveDeltaPars(position); } std::cout<<"INFO in LauFlavTag::addTagger : Added tagger with name "<< name << std::endl; } void LauFlavTag::useAveDeltaPars(Int_t position) { TString tagEff_ave("tagEff_ave_"+taggerPosition_[position]); TString tagEff_delta("tagEff_delta_"+taggerPosition_[position]); TString calib_p0_ave("calib_p0_ave_"+taggerPosition_[position]); TString calib_p0_delta("calib_p0_delta_"+taggerPosition_[position]); TString calib_p1_ave("calib_p1_ave_"+taggerPosition_[position]); TString calib_p1_delta("calib_p1_delta_"+taggerPosition_[position]); Double_t ave = ((tagEff_B0_[position]->unblindValue() + tagEff_B0bar_[position]->unblindValue())/2); LauParameter* tageffave = new LauParameter(tagEff_ave,ave,0.0,1.0,kTRUE); tagEff_ave_.push_back(tageffave); tagEff_ave_[position]->initValue(ave); tagEff_ave_[position]->genValue(ave); tagEff_ave_[position]->fixed(tagEff_B0_[position]->fixed()); //Update once full code in place Double_t delta = (tagEff_B0_[position]->unblindValue() - tagEff_B0bar_[position]->unblindValue()); LauParameter* tageffdelta = new LauParameter(tagEff_delta,delta,0.0,1.0,kTRUE); tagEff_delta_.push_back(tageffdelta); tagEff_delta_[position]->initValue(delta); tagEff_delta_[position]->genValue(delta); tagEff_delta_[position]->fixed(tagEff_B0_[position]->fixed()); //Update once full code in place ave = ((calib_p0_B0_[position]->unblindValue() + calib_p0_B0bar_[position]->unblindValue())/2); LauParameter* calibp0ave = new LauParameter(calib_p0_ave,ave,0.0,1.0,kTRUE); calib_p0_ave_.push_back(calibp0ave); calib_p0_ave_[position]->initValue(ave); calib_p0_ave_[position]->genValue(ave); calib_p0_ave_[position]->fixed(calib_p0_B0_[position]->fixed()); //Update once full code in place delta = (calib_p0_B0_[position]->unblindValue() - calib_p0_B0bar_[position]->unblindValue()); LauParameter* calibp0delta = new LauParameter(calib_p0_delta,delta,0.0,1.0,kTRUE); calib_p0_delta_.push_back(calibp0delta); calib_p0_delta_[position]->initValue(delta); calib_p0_delta_[position]->genValue(delta); calib_p0_delta_[position]->fixed(calib_p0_B0_[position]->fixed()); //Update once full code in place ave = ((calib_p1_B0_[position]->unblindValue() + calib_p1_B0bar_[position]->unblindValue())/2); LauParameter* calibp1ave = new LauParameter(calib_p1_ave,ave,0.0,1.0,kTRUE); calib_p1_ave_.push_back(calibp1ave); calib_p1_ave_[position]->initValue(ave); calib_p1_ave_[position]->genValue(ave); calib_p1_ave_[position]->fixed(calib_p1_B0_[position]->fixed()); //Update once full code in place delta = (calib_p1_B0_[position]->unblindValue() - calib_p1_B0bar_[position]->unblindValue()); LauParameter* calibp1delta = new LauParameter(calib_p1_delta,delta,0.0,1.0,kTRUE); calib_p1_delta_.push_back(calibp1delta); calib_p1_delta_[position]->initValue(delta); calib_p1_delta_[position]->genValue(delta); calib_p1_delta_[position]->fixed(calib_p1_B0_[position]->fixed()); //Update once full code in place } void LauFlavTag::cacheInputFitVars(LauFitDataTree* inputFitData) { evtTagFlv_.clear(); evtMistag_.clear(); evtTrueTagFlv_.clear(); // Loop over the taggers to check the branches for (UInt_t i=0; i < tagVarName_.size(); ++i){ if ( ! inputFitData->haveBranch( tagVarName_[i] ) ) { std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << tagVarName_[i] << "\"." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( ! inputFitData->haveBranch( mistagVarName_[i] ) ) { std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << mistagVarName_[i] << "\"." << std::endl; gSystem->Exit(EXIT_FAILURE); } } if ( ! inputFitData->haveBranch( trueTagVarName_ ) ) { std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << trueTagVarName_ << "\"." << std::endl; gSystem->Exit(EXIT_FAILURE); } UInt_t nEvents = inputFitData->nEvents(); evtTagFlv_.reserve( nEvents ); evtMistag_.reserve( nEvents ); evtTrueTagFlv_.reserve( nEvents ); LauFitData::const_iterator fitdata_iter; for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) { const LauFitData& dataValues = inputFitData->getData(iEvt); //Clear vectors curEvtTagFlv_.clear(); curEvtMistag_.clear(); // For untagged events see if we have a truth tag for normalisation modes curEvtTrueTagFlv_ = 0; fitdata_iter = dataValues.find( trueTagVarName_ ); curEvtTrueTagFlv_ = static_cast( fitdata_iter->second ); if ( curEvtTrueTagFlv_ > 1 ) { std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid true tagging output " << curEvtTrueTagFlv_ << " for event " << iEvt << ", setting it to +1" << std::endl; curEvtTrueTagFlv_ = +1; } else if ( curEvtTrueTagFlv_ < -1 ){ std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid true tagging output " << curEvtTrueTagFlv_ << " for event " << iEvt << ", setting it to -1" << std::endl; curEvtTrueTagFlv_ = -1; } evtTrueTagFlv_.push_back(curEvtTrueTagFlv_); // TODO - surely this part should only be done if we have per-event mistagging for (UInt_t i=0; i < tagVarName_.size(); ++i){ fitdata_iter = dataValues.find( tagVarName_[i] ); curEvtTagFlv_.push_back(static_cast( fitdata_iter->second )); if ( curEvtTagFlv_[i] > 1 ) { std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv_[i] << " for event " << iEvt << ", setting it to +1" << std::endl; curEvtTagFlv_[i] = +1; } else if ( curEvtTagFlv_[i] < -1 ) { std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv_[i] << " for event " << iEvt << ", setting it to -1" << std::endl; curEvtTagFlv_[i] = -1; } fitdata_iter = dataValues.find( mistagVarName_[i]); curEvtMistag_.push_back(static_cast( fitdata_iter->second )); if (curEvtMistag_[i] > 0.5){ std::cerr<<"WARNING in LauFlavTag::cacheInputFitVars : Mistag value "<value( calib_p0_ave_[position]->unblindValue() + 0.5*calib_p0_delta_[position]->unblindValue() ); calib_p0_B0bar_[position]->value( calib_p0_ave_[position]->unblindValue() - 0.5*calib_p0_delta_[position]->unblindValue() ); calib_p1_B0_[position]->value( calib_p1_ave_[position]->unblindValue() + 0.5*calib_p1_delta_[position]->unblindValue() ); calib_p1_B0bar_[position]->value( calib_p1_ave_[position]->unblindValue() - 0.5*calib_p1_delta_[position]->unblindValue() ); } if (flag == 1){ return calib_p0_B0_[position]->unblindValue() + calib_p1_B0_[position]->unblindValue() * (curEvtMistag_[position] - perEvtAvgMistag_[position]); } else{ return calib_p0_B0bar_[position]->unblindValue() + calib_p1_B0bar_[position]->unblindValue() * (curEvtMistag_[position] - perEvtAvgMistag_[position]); } std::cerr << "ERROR in LauFlavTag::getOmega : Current event flavour tag not defined" << std::endl; return 0.0; } Double_t LauFlavTag::getCapitalOmega(const Int_t position, const Int_t flag) //const { if (TMath::Abs(flag) != 1){ 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] == -1){ deltam1 = 1; } else if(curEvtTagFlv_[position] == 1){ deltap1 = 1; } else{ delta0 = 1; } //Efficiency Double_t eff=0.0; if (useAveDelta_){ tagEff_B0_[position]->value( tagEff_ave_[position]->unblindValue() + 0.5*tagEff_delta_[position]->unblindValue() ); tagEff_B0bar_[position]->value( tagEff_ave_[position]->unblindValue() - 0.5*tagEff_delta_[position]->unblindValue() ); } if (flag==1){ eff = tagEff_B0_[position]->unblindValue(); } else { eff = tagEff_B0bar_[position]->unblindValue(); } //Little omega Double_t omega = this->getOmega(position, flag); //eta PDF value std::vector abs; abs.push_back(curEvtMistag_[position]); etaPdfs_[position]->calcLikelihoodInfo(abs); Double_t h = etaPdfs_[position]->getLikelihood(); //std::cout << deltap1 << " " << deltam1 << " " << delta0 << " " << eff << " " << omega << " " << h << std::endl; //Put it together if (flag == 1){ //Particle return (deltap1*eff*(1-omega) + deltam1*eff*omega)*h + delta0*(1-eff); } else { return (deltam1*eff*(1-omega) + deltap1*eff*omega)*h + delta0*(1-eff); } } Double_t LauFlavTag::getOmegaGen(const Int_t position, const Double_t eta, const Int_t flag) //const { if (TMath::Abs(flag) != 1){ std::cerr << "ERROR in LauFlavTag::getOmegaGen : Invalid flag, you must request either omega (-1) or omega bar (+1) to be returned" << std::endl; return 0.0; } //If we are floating average omega and delta omega we need to use those parameters instead if (useAveDelta_){ calib_p0_B0_[position]->value( calib_p0_ave_[position]->unblindValue() + 0.5*calib_p0_delta_[position]->unblindValue() ); calib_p0_B0bar_[position]->value( calib_p0_ave_[position]->unblindValue() - 0.5*calib_p0_delta_[position]->unblindValue() ); calib_p1_B0_[position]->value( calib_p1_ave_[position]->unblindValue() + 0.5*calib_p1_delta_[position]->unblindValue() ); calib_p1_B0bar_[position]->value( calib_p1_ave_[position]->unblindValue() - 0.5*calib_p1_delta_[position]->unblindValue() ); } if (flag == 1){ return calib_p0_B0_[position]->unblindValue() + calib_p1_B0_[position]->unblindValue() * (eta - perEvtAvgMistag_[position]); } else{ return calib_p0_B0bar_[position]->unblindValue() + calib_p1_B0bar_[position]->unblindValue() * (eta - perEvtAvgMistag_[position]); } std::cerr << "ERROR in LauFlavTag::getOmegaGen : Current event flavour tag not defined" << std::endl; return 0.0; } Double_t LauFlavTag::getCapitalOmegaGen(const Int_t position, const Double_t eta, const Int_t curEvtTagFlv, const Int_t flag) { if (TMath::Abs(flag) != 1){ std::cerr << "ERROR in LauFlavTag::getCapitalOmegaGen : 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 == -1){ deltam1 = 1; } else if(curEvtTagFlv == 1){ deltap1 = 1; } else{ delta0 = 1; } //Efficiency Double_t eff=0.0; if (useAveDelta_){ tagEff_B0_[position]->value( tagEff_ave_[position]->unblindValue() + 0.5*tagEff_delta_[position]->unblindValue() ); tagEff_B0bar_[position]->value( tagEff_ave_[position]->unblindValue() - 0.5*tagEff_delta_[position]->unblindValue() ); } if (flag==1){ eff = tagEff_B0_[position]->unblindValue(); } else { eff = tagEff_B0bar_[position]->unblindValue(); } //Little omega Double_t omega = this->getOmegaGen(position, eta, flag); //eta PDF value std::vector abs; abs.push_back(eta); etaPdfs_[position]->calcLikelihoodInfo(abs); Double_t h = etaPdfs_[position]->getLikelihood(); //Put it together if (flag == -1){ return (deltam1*eff*(1-omega) + deltap1*eff*omega)*h + delta0*(1-eff)*0.5; } else { return (deltap1*eff*(1-omega) + deltam1*eff*omega)*h + delta0*(1-eff)*0.5; } //return(0.); } Double_t LauFlavTag::getEtaGen(Int_t position) const { LauFitData data = etaPdfs_[position]->generate(0); //TODO Add DP dependence? Double_t etagen = data.find(etaPdfs_[position]->varName())->second; if (etagen > 0.5){etagen = 0.5;} if (etagen < 0.0){etagen = 0.0;} return etagen; } void LauFlavTag::setTrueTagVarName(TString trueTagVarName){ trueTagVarName_ = trueTagVarName; }