TopJetCombination/src/TtSemiLepHypothesis.cc

0001 #include "CommonTools/CandUtils/interface/AddFourMomenta.h"
0002 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
0003 #include "DataFormats/PatCandidates/interface/Electron.h"
0004 #include "DataFormats/PatCandidates/interface/Muon.h"
0005 #include "TopQuarkAnalysis/TopTools/interface/MEzCalculator.h"
0006
0007 #include "TopQuarkAnalysis/TopJetCombination/interface/TtSemiLepHypothesis.h"
0008
0009 /// default constructor
0010 TtSemiLepHypothesis::TtSemiLepHypothesis(const edm::ParameterSet& cfg)
0011     : jetsToken_(consumes<std::vector<pat::Jet> >(cfg.getParameter<edm::InputTag>("jets"))),
0012       lepsToken_(consumes<edm::View<reco::RecoCandidate> >(cfg.getParameter<edm::InputTag>("leps"))),
0013       metsToken_(consumes<std::vector<pat::MET> >(cfg.getParameter<edm::InputTag>("mets"))),
0014       nJetsConsideredToken_(consumes<int>(cfg.getParameter<edm::InputTag>("nJetsConsidered"))),
0015       numberOfRealNeutrinoSolutions_(-1) {
0016   getMatch_ = false;
0017   if (cfg.exists("match")) {
0018     getMatch_ = true;
0019     matchToken_ = consumes<std::vector<std::vector<int> > >(cfg.getParameter<edm::InputTag>("match"));
0020   }
0021   if (cfg.exists("neutrinoSolutionType"))
0022     neutrinoSolutionType_ = cfg.getParameter<int>("neutrinoSolutionType");
0023   else
0024     neutrinoSolutionType_ = -1;
0025   if (cfg.exists("jetCorrectionLevel")) {
0026     jetCorrectionLevel_ = cfg.getParameter<std::string>("jetCorrectionLevel");
0027   }
0028   produces<std::vector<std::pair<reco::CompositeCandidate, std::vector<int> > > >();
0029   produces<int>("Key");
0030   produces<int>("NumberOfRealNeutrinoSolutions");
0031   produces<int>("NumberOfConsideredJets");
0032 }
0033
0034 /// produce the event hypothesis as CompositeCandidate and Key
0035 void TtSemiLepHypothesis::produce(edm::Event& evt, const edm::EventSetup& setup) {
0036   edm::Handle<std::vector<pat::Jet> > jets;
0037   evt.getByToken(jetsToken_, jets);
0038
0039   edm::Handle<edm::View<reco::RecoCandidate> > leps;
0040   evt.getByToken(lepsToken_, leps);
0041
0042   edm::Handle<std::vector<pat::MET> > mets;
0043   evt.getByToken(metsToken_, mets);
0044
0045   edm::Handle<int> nJetsConsidered;
0046   evt.getByToken(nJetsConsideredToken_, nJetsConsidered);
0047
0048   std::vector<std::vector<int> > matchVec;
0049   if (getMatch_) {
0050     edm::Handle<std::vector<std::vector<int> > > matchHandle;
0051     evt.getByToken(matchToken_, matchHandle);
0052     matchVec = *matchHandle;
0053   } else {
0054     std::vector<int> dummyMatch;
0055     for (unsigned int i = 0; i < 4; ++i)
0056       dummyMatch.push_back(-1);
0057     matchVec.push_back(dummyMatch);
0058   }
0059
0060   // declare unique_ptr for products
0061   std::unique_ptr<std::vector<std::pair<reco::CompositeCandidate, std::vector<int> > > > pOut(
0062       new std::vector<std::pair<reco::CompositeCandidate, std::vector<int> > >);
0063   std::unique_ptr<int> pKey(new int);
0064   std::unique_ptr<int> pNeutrinoSolutions(new int);
0065   std::unique_ptr<int> pJetsConsidered(new int);
0066
0067   // go through given vector of jet combinations
0068   unsigned int idMatch = 0;
0069   typedef std::vector<std::vector<int> >::iterator MatchVecIterator;
0070   for (MatchVecIterator match = matchVec.begin(); match != matchVec.end(); ++match) {
0071     // reset pointers
0072     resetCandidates();
0073     // build hypothesis
0074     buildHypo(evt, leps, mets, jets, *match, idMatch++);
0075     pOut->push_back(std::make_pair(hypo(), *match));
0076   }
0077   // feed out hyps and matches
0078   evt.put(std::move(pOut));
0079
0080   // build and feed out key
0081   buildKey();
0082   *pKey = key();
0083   evt.put(std::move(pKey), "Key");
0084
0085   // feed out number of real neutrino solutions
0086   *pNeutrinoSolutions = numberOfRealNeutrinoSolutions_;
0087   evt.put(std::move(pNeutrinoSolutions), "NumberOfRealNeutrinoSolutions");
0088
0089   // feed out number of considered jets
0090   *pJetsConsidered = *nJetsConsidered;
0091   evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0092 }
0093
0094 /// reset candidate pointers before hypo build process
0095 void TtSemiLepHypothesis::resetCandidates() {
0096   numberOfRealNeutrinoSolutions_ = -1;
0097   lightQ_.reset();
0098   lightQBar_.reset();
0099   hadronicB_.reset();
0100   leptonicB_.reset();
0101   neutrino_.reset();
0102   lepton_.reset();
0103 }
0104
0105 /// return event hypothesis
0106 reco::CompositeCandidate TtSemiLepHypothesis::hypo() {
0107   // check for sanity of the hypothesis
0108   if (!lightQ_ || !lightQBar_ || !hadronicB_ || !leptonicB_ || !neutrino_ || !lepton_)
0109     return reco::CompositeCandidate();
0110
0111   // setup transient references
0112   reco::CompositeCandidate hyp, hadTop, hadW, lepTop, lepW;
0113
0114   AddFourMomenta addFourMomenta;
0115   // build up the top branch that decays leptonically
0116   lepW.addDaughter(std::move(lepton_), TtSemiLepDaughter::Lep);
0117   lepW.addDaughter(std::move(neutrino_), TtSemiLepDaughter::Nu);
0118   addFourMomenta.set(lepW);
0119   lepTop.addDaughter(lepW, TtSemiLepDaughter::LepW);
0120   lepTop.addDaughter(std::move(leptonicB_), TtSemiLepDaughter::LepB);
0121   addFourMomenta.set(lepTop);
0122
0123   // build up the top branch that decays hadronically
0124   hadW.addDaughter(std::move(lightQ_), TtSemiLepDaughter::HadP);
0125   hadW.addDaughter(std::move(lightQBar_), TtSemiLepDaughter::HadQ);
0126   addFourMomenta.set(hadW);
0127   hadTop.addDaughter(hadW, TtSemiLepDaughter::HadW);
0128   hadTop.addDaughter(std::move(hadronicB_), TtSemiLepDaughter::HadB);
0129   addFourMomenta.set(hadTop);
0130
0131   // build ttbar hypotheses
0132   hyp.addDaughter(lepTop, TtSemiLepDaughter::LepTop);
0133   hyp.addDaughter(hadTop, TtSemiLepDaughter::HadTop);
0134   addFourMomenta.set(hyp);
0135
0136   return hyp;
0137 }
0138
0139 /// determine lepton type of reco candidate and return a corresponding WDecay::LeptonType; the type is kNone if it is whether a muon nor an electron
0140 WDecay::LeptonType TtSemiLepHypothesis::leptonType(const reco::RecoCandidate* cand) {
0141   // check whetherwe are dealing with a reco muon or a reco electron
0142   WDecay::LeptonType type = WDecay::kNone;
0143   if (dynamic_cast<const reco::Muon*>(cand)) {
0144     type = WDecay::kMuon;
0145   } else if (dynamic_cast<const reco::GsfElectron*>(cand)) {
0146     type = WDecay::kElec;
0147   }
0148   return type;
0149 }
0150
0151 /// helper function to construct the proper correction level string for corresponding quarkType
0152 std::string TtSemiLepHypothesis::jetCorrectionLevel(const std::string& quarkType) {
0153   // jetCorrectionLevel was not configured
0154   if (jetCorrectionLevel_.empty())
0155     throw cms::Exception("Configuration")
0156         << "Unconfigured jetCorrectionLevel. Please use an appropriate, non-empty string.\n";
0157
0158   // quarkType is unknown
0159   if (!(quarkType == "wQuarkMix" || quarkType == "udsQuark" || quarkType == "cQuark" || quarkType == "bQuark"))
0160     throw cms::Exception("Configuration") << quarkType << " is unknown as a quarkType for the jetCorrectionLevel.\n";
0161
0162   // combine correction level; start with a ':' even if
0163   // there is no flavor tag to be added, as it is needed
0164   // by makeCandidate to disentangle the correction tag
0165   // from a potential flavor tag, which can be empty
0166   std::string level = jetCorrectionLevel_ + ":";
0167   if (level == "L5Flavor:" || level == "L6UE:" || level == "L7Parton:") {
0168     if (quarkType == "wQuarkMix") {
0169       level += "wMix";
0170     }
0171     if (quarkType == "udsQuark") {
0172       level += "uds";
0173     }
0174     if (quarkType == "cQuark") {
0175       level += "charm";
0176     }
0177     if (quarkType == "bQuark") {
0178       level += "bottom";
0179     }
0180   } else {
0181     level += "none";
0182   }
0183   return level;
0184 }
0185
0186 /// use one object in a jet collection to set a ShallowClonePtrCandidate with proper jet corrections
0187 std::unique_ptr<reco::ShallowClonePtrCandidate> TtSemiLepHypothesis::makeCandidate(
0188     const edm::Handle<std::vector<pat::Jet> >& handle, const int& idx, const std::string& correctionLevel) {
0189   edm::Ptr<pat::Jet> ptr = edm::Ptr<pat::Jet>(handle, idx);
0190   // disentangle the correction from the potential flavor tag
0191   // by the separating ':'; the flavor tag can be empty though
0192   std::string step = correctionLevel.substr(0, correctionLevel.find(':'));
0193   std::string flavor = correctionLevel.substr(1 + correctionLevel.find(':'));
0194   float corrFactor = 1.;
0195   if (flavor == "wMix")
0196     corrFactor = 0.75 * ptr->jecFactor(step, "uds") + 0.25 * ptr->jecFactor(step, "charm");
0197   else
0198     corrFactor = ptr->jecFactor(step, flavor);
0199   return std::make_unique<reco::ShallowClonePtrCandidate>(ptr, ptr->charge(), ptr->p4() * corrFactor, ptr->vertex());
0200 }
0201
0202 /// set neutrino, using mW = 80.4 to calculate the neutrino pz
0203 void TtSemiLepHypothesis::setNeutrino(const edm::Handle<std::vector<pat::MET> >& met,
0204                                       const edm::Handle<edm::View<reco::RecoCandidate> >& leps,
0205                                       const int& idx,
0206                                       const int& type) {
0207   edm::Ptr<pat::MET> ptr = edm::Ptr<pat::MET>(met, 0);
0208   MEzCalculator mez;
0209   mez.SetMET(*(met->begin()));
0210   if (leptonType(&(leps->front())) == WDecay::kMuon)
0211     mez.SetLepton((*leps)[idx], true);
0212   else if (leptonType(&(leps->front())) == WDecay::kElec)
0213     mez.SetLepton((*leps)[idx], false);
0214   else
0215     throw cms::Exception("UnimplementedFeature")
0216         << "Type of lepton given together with MET for solving neutrino kinematics is neither muon nor electron.\n";
0217   double pz = mez.Calculate(type);
0218   numberOfRealNeutrinoSolutions_ = mez.IsComplex() ? 0 : 2;
0219   const math::XYZTLorentzVector p4(
0220       ptr->px(), ptr->py(), pz, sqrt(ptr->px() * ptr->px() + ptr->py() * ptr->py() + pz * pz));
0221   neutrino_ = std::make_unique<reco::ShallowClonePtrCandidate>(ptr, ptr->charge(), p4, ptr->vertex());
0222 }
0223
0224 /// minimalistic build function for simple hypotheses
0225 void TtSemiLepHypothesis::buildHypo(const edm::Handle<edm::View<reco::RecoCandidate> >& leps,
0226                                     const edm::Handle<std::vector<pat::MET> >& mets,
0227                                     const edm::Handle<std::vector<pat::Jet> >& jets,
0228                                     std::vector<int>& match) {
0229   // -----------------------------------------------------
0230   // add jets
0231   // -----------------------------------------------------
0232   for (unsigned idx = 0; idx < match.size(); ++idx) {
0233     if (isValid(match[idx], jets)) {
0234       switch (idx) {
0235         case TtSemiLepEvtPartons::LightQ:
0236           lightQ_ = makeCandidate(jets, match[idx], jetCorrectionLevel("wQuarkMix"));
0237           break;
0238         case TtSemiLepEvtPartons::LightQBar:
0239           lightQBar_ = makeCandidate(jets, match[idx], jetCorrectionLevel("wQuarkMix"));
0240           break;
0241         case TtSemiLepEvtPartons::HadB:
0242           hadronicB_ = makeCandidate(jets, match[idx], jetCorrectionLevel("bQuark"));
0243           break;
0244         case TtSemiLepEvtPartons::LepB:
0245           leptonicB_ = makeCandidate(jets, match[idx], jetCorrectionLevel("bQuark"));
0246           break;
0247       }
0248     }
0249   }
0250
0251   // -----------------------------------------------------
0252   // add lepton
0253   // -----------------------------------------------------
0254   if (leps->empty())
0255     return;
0256   lepton_ = makeCandidate(leps, 0);
0257   match.push_back(0);
0258
0259   // -----------------------------------------------------
0260   // add neutrino
0261   // -----------------------------------------------------
0262   if (mets->empty())
0263     return;
0264   if (neutrinoSolutionType_ == -1)
0265     neutrino_ = makeCandidate(mets, 0);
0266   else
0267     setNeutrino(mets, leps, 0, neutrinoSolutionType_);
0268 }