TopKinFitter/interface/TtSemiLepKinFitter.h

0001 #ifndef TtSemiLepKinFitter_h
0002 #define TtSemiLepKinFitter_h
0003
0004 #include <vector>
0005
0006 #include "TLorentzVector.h"
0007
0008 #include "DataFormats/PatCandidates/interface/Lepton.h"
0009
0010 #include "AnalysisDataFormats/TopObjects/interface/TtSemiEvtSolution.h"
0011
0012 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
0013
0014 #include "TopQuarkAnalysis/TopKinFitter/interface/TopKinFitter.h"
0015 #include "TopQuarkAnalysis/TopKinFitter/interface/CovarianceMatrix.h"
0016
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018
0019 class TAbsFitParticle;
0020 class TFitConstraintM;
0021 class TFitConstraintEp;
0022
0023 /*
0024   \class   TtSemiLepKinFitter TtSemiLepKinFitter.h "TopQuarkAnalysis/TopKinFitter/interface/TtSemiLepKinFitter.h"
0025
0026   \brief   one line description to be added here...
0027
0028   text to be added here...
0029
0030 **/
0031
0032 class TtSemiLepKinFitter : public TopKinFitter {
0033 public:
0034   /// supported constraints
0035   enum Constraint { kWHadMass = 1, kWLepMass, kTopHadMass, kTopLepMass, kNeutrinoMass, kEqualTopMasses, kSumPt };
0036
0037 public:
0038   /// default constructor
0039   explicit TtSemiLepKinFitter();
0040   /// constructor initialized with built-in types and class enum's custom parameters
0041   explicit TtSemiLepKinFitter(Param jetParam,
0042                               Param lepParam,
0043                               Param metParam,
0044                               int maxNrIter,
0045                               double maxDeltaS,
0046                               double maxF,
0047                               const std::vector<Constraint>& constraints,
0048                               double mW = 80.4,
0049                               double mTop = 173.,
0050                               const std::vector<edm::ParameterSet>* udscResolutions = nullptr,
0051                               const std::vector<edm::ParameterSet>* bResolutions = nullptr,
0052                               const std::vector<edm::ParameterSet>* lepResolutions = nullptr,
0053                               const std::vector<edm::ParameterSet>* metResolutions = nullptr,
0054                               const std::vector<double>* jetEnergyResolutionScaleFactors = nullptr,
0055                               const std::vector<double>* jetEnergyResolutionEtaBinning = nullptr);
0056   /// default destructor
0057   ~TtSemiLepKinFitter();
0058
0059   /// kinematic fit interface for PAT objects
0060   template <class LeptonType>
0061   int fit(const std::vector<pat::Jet>& jets, const pat::Lepton<LeptonType>& leps, const pat::MET& met);
0062   /// kinematic fit interface for plain 4-vecs
0063   int fit(const TLorentzVector& p4HadP,
0064           const TLorentzVector& p4HadQ,
0065           const TLorentzVector& p4HadB,
0066           const TLorentzVector& p4LepB,
0067           const TLorentzVector& p4Lepton,
0068           const TLorentzVector& p4Neutrino,
0069           const int leptonCharge,
0070           const CovarianceMatrix::ObjectType leptonType);
0071   /// common core of the fit interface
0072   int fit(const TLorentzVector& p4HadP,
0073           const TLorentzVector& p4HadQ,
0074           const TLorentzVector& p4HadB,
0075           const TLorentzVector& p4LepB,
0076           const TLorentzVector& p4Lepton,
0077           const TLorentzVector& p4Neutrino,
0078           const TMatrixD& covHadP,
0079           const TMatrixD& covHadQ,
0080           const TMatrixD& covHadB,
0081           const TMatrixD& covLepB,
0082           const TMatrixD& covLepton,
0083           const TMatrixD& covNeutrino,
0084           const int leptonCharge);
0085   /// return hadronic b quark candidate
0086   const pat::Particle fittedHadB() const { return (fitter_->getStatus() == 0 ? fittedHadB_ : pat::Particle()); };
0087   /// return hadronic light quark candidate
0088   const pat::Particle fittedHadP() const { return (fitter_->getStatus() == 0 ? fittedHadP_ : pat::Particle()); };
0089   /// return hadronic light quark candidate
0090   const pat::Particle fittedHadQ() const { return (fitter_->getStatus() == 0 ? fittedHadQ_ : pat::Particle()); };
0091   /// return leptonic b quark candidate
0092   const pat::Particle fittedLepB() const { return (fitter_->getStatus() == 0 ? fittedLepB_ : pat::Particle()); };
0093   /// return lepton candidate
0094   const pat::Particle fittedLepton() const { return (fitter_->getStatus() == 0 ? fittedLepton_ : pat::Particle()); };
0095   /// return neutrino candidate
0096   const pat::Particle fittedNeutrino() const {
0097     return (fitter_->getStatus() == 0 ? fittedNeutrino_ : pat::Particle());
0098   };
0099   /// add kin fit information to the old event solution (in for legacy reasons)
0100   TtSemiEvtSolution addKinFitInfo(TtSemiEvtSolution* asol);
0101
0102 private:
0103   /// print fitter setup
0104   void printSetup() const;
0105   /// setup fitter
0106   void setupFitter();
0107   /// initialize jet inputs
0108   void setupJets();
0109   /// initialize lepton inputs
0110   void setupLeptons();
0111   /// initialize constraints
0112   void setupConstraints();
0113
0114 private:
0115   /// input particles
0116   std::unique_ptr<TAbsFitParticle> hadB_;
0117   std::unique_ptr<TAbsFitParticle> hadP_;
0118   std::unique_ptr<TAbsFitParticle> hadQ_;
0119   std::unique_ptr<TAbsFitParticle> lepB_;
0120   std::unique_ptr<TAbsFitParticle> lepton_;
0121   std::unique_ptr<TAbsFitParticle> neutrino_;
0122   /// resolutions
0123   const std::vector<edm::ParameterSet>* udscResolutions_ = nullptr;
0124   const std::vector<edm::ParameterSet>* bResolutions_ = nullptr;
0125   const std::vector<edm::ParameterSet>* lepResolutions_ = nullptr;
0126   const std::vector<edm::ParameterSet>* metResolutions_ = nullptr;
0127   /// scale factors for the jet energy resolution
0128   const std::vector<double>* jetEnergyResolutionScaleFactors_ = nullptr;
0129   const std::vector<double>* jetEnergyResolutionEtaBinning_ = nullptr;
0130   /// object used to construct the covariance matrices for the individual particles
0131   std::unique_ptr<CovarianceMatrix> covM_;
0132   /// supported constraints
0133   std::map<Constraint, std::unique_ptr<TFitConstraintM>> massConstr_;
0134   std::unique_ptr<TFitConstraintEp> sumPxConstr_;
0135   std::unique_ptr<TFitConstraintEp> sumPyConstr_;
0136   /// output particles
0137   pat::Particle fittedHadB_;
0138   pat::Particle fittedHadP_;
0139   pat::Particle fittedHadQ_;
0140   pat::Particle fittedLepB_;
0141   pat::Particle fittedLepton_;
0142   pat::Particle fittedNeutrino_;
0143   /// jet parametrization
0144   Param jetParam_;
0145   /// lepton parametrization
0146   Param lepParam_;
0147   /// met parametrization
0148   Param metParam_;
0149   /// vector of constraints to be used
0150   std::vector<Constraint> constrList_;
0151   /// internally use simple boolean for this constraint to reduce the per-event computing time
0152   bool constrainSumPt_;
0153 };
0154
0155 template <class LeptonType>
0156 int TtSemiLepKinFitter::fit(const std::vector<pat::Jet>& jets,
0157                             const pat::Lepton<LeptonType>& lepton,
0158                             const pat::MET& neutrino) {
0159   if (jets.size() < 4)
0160     throw edm::Exception(edm::errors::Configuration, "Cannot run the TtSemiLepKinFitter with less than 4 jets");
0161
0162   // get jets in right order
0163   const pat::Jet& hadP = jets[TtSemiLepEvtPartons::LightQ];
0164   const pat::Jet& hadQ = jets[TtSemiLepEvtPartons::LightQBar];
0165   const pat::Jet& hadB = jets[TtSemiLepEvtPartons::HadB];
0166   const pat::Jet& lepB = jets[TtSemiLepEvtPartons::LepB];
0167
0168   // initialize particles
0169   const TLorentzVector p4HadP(hadP.px(), hadP.py(), hadP.pz(), hadP.energy());
0170   const TLorentzVector p4HadQ(hadQ.px(), hadQ.py(), hadQ.pz(), hadQ.energy());
0171   const TLorentzVector p4HadB(hadB.px(), hadB.py(), hadB.pz(), hadB.energy());
0172   const TLorentzVector p4LepB(lepB.px(), lepB.py(), lepB.pz(), lepB.energy());
0173   const TLorentzVector p4Lepton(lepton.px(), lepton.py(), lepton.pz(), lepton.energy());
0174   const TLorentzVector p4Neutrino(neutrino.px(), neutrino.py(), 0, neutrino.et());
0175
0176   // initialize covariance matrices
0177   TMatrixD covHadP = covM_->setupMatrix(hadP, jetParam_);
0178   TMatrixD covHadQ = covM_->setupMatrix(hadQ, jetParam_);
0179   TMatrixD covHadB = covM_->setupMatrix(hadB, jetParam_, "bjets");
0180   TMatrixD covLepB = covM_->setupMatrix(lepB, jetParam_, "bjets");
0181   TMatrixD covLepton = covM_->setupMatrix(lepton, lepParam_);
0182   TMatrixD covNeutrino = covM_->setupMatrix(neutrino, metParam_);
0183
0184   // now do the part that is fully independent of PAT features
0185   return fit(p4HadP,
0186              p4HadQ,
0187              p4HadB,
0188              p4LepB,
0189              p4Lepton,
0190              p4Neutrino,
0191              covHadP,
0192              covHadQ,
0193              covHadB,
0194              covLepB,
0195              covLepton,
0196              covNeutrino,
0197              lepton.charge());
0198 }
0199
0200 #endif