TopEventSelection/plugins/TtSemiLepSignalSelMVAComputer.cc

0001 #include "PhysicsTools/JetMCUtils/interface/combination.h"
0002
0003 #include "TopQuarkAnalysis/TopEventSelection/plugins/TtSemiLepSignalSelMVAComputer.h"
0004 #include "TopQuarkAnalysis/TopEventSelection/interface/TtSemiLepSignalSelEval.h"
0005
0006 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0007
0008 #include "DataFormats/Common/interface/TriggerResults.h"
0009 #include "FWCore/Framework/interface/TriggerNamesService.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011 #include "DataFormats/PatCandidates/interface/Flags.h"
0012
0013 TtSemiLepSignalSelMVAComputer::TtSemiLepSignalSelMVAComputer(const edm::ParameterSet& cfg)
0014     : mvaToken_(esConsumes()),
0015       muonsToken_(consumes<edm::View<pat::Muon> >(cfg.getParameter<edm::InputTag>("muons"))),
0016       jetsToken_(consumes<std::vector<pat::Jet> >(cfg.getParameter<edm::InputTag>("jets"))),
0017       METsToken_(consumes<edm::View<pat::MET> >(cfg.getParameter<edm::InputTag>("mets"))),
0018       electronsToken_(consumes<edm::View<pat::Electron> >(cfg.getParameter<edm::InputTag>("elecs"))),
0019       putToken_(produces("DiscSel")) {}
0020
0021 TtSemiLepSignalSelMVAComputer::~TtSemiLepSignalSelMVAComputer() {}
0022
0023 void TtSemiLepSignalSelMVAComputer::produce(edm::Event& evt, const edm::EventSetup& setup) {
0024   mvaComputer.update(&setup.getData(mvaToken_), "ttSemiLepSignalSelMVA");
0025
0026   //make your preselection! This must!! be the same one as in TraintreeSaver.cc
0027   edm::Handle<edm::View<pat::MET> > MET_handle;
0028   evt.getByToken(METsToken_, MET_handle);
0029   if (!MET_handle.isValid())
0030     return;
0031   const edm::View<pat::MET> MET = *MET_handle;
0032
0033   edm::Handle<std::vector<pat::Jet> > jet_handle;
0034   evt.getByToken(jetsToken_, jet_handle);
0035   if (!jet_handle.isValid())
0036     return;
0037   const std::vector<pat::Jet> jets = *jet_handle;
0038   unsigned int nJets = 0;
0039   std::vector<pat::Jet> seljets;
0040   for (std::vector<pat::Jet>::const_iterator it = jets.begin(); it != jets.end(); it++) {
0041     if (!(pat::Flags::test(*it, pat::Flags::Overlap::Electrons)))
0042       continue;
0043     if (it->pt() > 30. && fabs(it->eta()) < 2.4) {
0044       seljets.push_back(*it);
0045       nJets++;
0046     }
0047   }
0048
0049   edm::Handle<edm::View<pat::Muon> > muon_handle;
0050   evt.getByToken(muonsToken_, muon_handle);
0051   if (!muon_handle.isValid())
0052     return;
0053   const edm::View<pat::Muon> muons = *muon_handle;
0054   int nmuons = 0;
0055   std::vector<pat::Muon> selMuons;
0056   for (edm::View<pat::Muon>::const_iterator it = muons.begin(); it != muons.end(); it++) {
0057     reco::TrackRef gltr = it->track();       // global track
0058     reco::TrackRef trtr = it->innerTrack();  // tracker track
0059     if (it->pt() > 30 && fabs(it->eta()) < 2.1 && (it->pt() / (it->pt() + it->trackIso() + it->caloIso())) > 0.95 &&
0060         it->isGlobalMuon()) {
0061       if (gltr.isNull())
0062         continue;  //temporary problems with dead trackrefs
0063       if ((gltr->chi2() / gltr->ndof()) < 10 && trtr->numberOfValidHits() > 11) {
0064         double dRmin = 9999.;
0065         for (std::vector<pat::Jet>::const_iterator ajet = seljets.begin(); ajet != seljets.end(); ajet++) {
0066           math::XYZTLorentzVector jet = ajet->p4();
0067           math::XYZTLorentzVector muon = it->p4();
0068           double tmpdR = DeltaR(muon, jet);
0069           if (tmpdR < dRmin)
0070             dRmin = tmpdR;
0071         }
0072         if (dRmin > 0.3) {  //temporary problems with muon isolation
0073           nmuons++;
0074           selMuons.push_back(*it);
0075         }
0076       }
0077     }
0078   }
0079
0080   edm::Handle<edm::View<pat::Electron> > electron_handle;
0081   evt.getByToken(electronsToken_, electron_handle);
0082   if (!electron_handle.isValid())
0083     return;
0084   const edm::View<pat::Electron> electrons = *electron_handle;
0085   int nelectrons = 0;
0086   for (edm::View<pat::Electron>::const_iterator it = electrons.begin(); it != electrons.end(); it++) {
0087     if (it->pt() > 30 && fabs(it->eta()) < 2.4 && (it->pt() / (it->pt() + it->trackIso() + it->caloIso())) > 0.95 &&
0088         it->isElectronIDAvailable("eidTight")) {
0089       if (it->electronID("eidTight") == 1)
0090         nelectrons++;
0091     }
0092   }
0093
0094   double discrim;
0095   // discriminator output for events which do not pass the preselection is set to -1
0096   if (nmuons != 1 || nJets < 4 || nelectrons > 0)
0097     discrim = -1.;  //std::cout<<"nJets: "<<seljets.size()<<" numLeptons: "<<nleptons<<std::endl;}
0098   else {
0099     //check wheter a event was already selected (problem with duplicated events)
0100     math::XYZTLorentzVector muon = selMuons.begin()->p4();
0101
0102     TtSemiLepSignalSel selection(seljets, muon, MET);
0103
0104     discrim = evaluateTtSemiLepSignalSel(mvaComputer, selection);
0105   }
0106
0107   evt.emplace(putToken_, discrim);
0108 }
0109
0110 double TtSemiLepSignalSelMVAComputer::DeltaPhi(const math::XYZTLorentzVector& v1, const math::XYZTLorentzVector& v2) {
0111   double dPhi = fabs(v1.Phi() - v2.Phi());
0112   if (dPhi > TMath::Pi())
0113     dPhi = 2 * TMath::Pi() - dPhi;
0114   return dPhi;
0115 }
0116
0117 double TtSemiLepSignalSelMVAComputer::DeltaR(const math::XYZTLorentzVector& v1, const math::XYZTLorentzVector& v2) {
0118   double dPhi = DeltaPhi(v1, v2);
0119   double dR = TMath::Sqrt((v1.Eta() - v2.Eta()) * (v1.Eta() - v2.Eta()) + dPhi * dPhi);
0120   return dR;
0121 }
0122
0123 // implement the plugins for the computer container
0124 // -> register TtSemiLepSignalSelMVARcd
0125 // -> define TtSemiLepSignalSelMVAFileSource
0126 MVA_COMPUTER_CONTAINER_IMPLEMENT(TtSemiLepSignalSelMVA);