L1TNtuples/src/L1AnalysisRecoMuon2.cc

0001 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoMuon2.h"
0002 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0003 #include <DataFormats/PatCandidates/interface/Muon.h>
0004 #include "L1Trigger/L1TNtuples/interface/MuonID.h"
0005 #include <TLorentzVector.h>
0006
0007 using namespace std;
0008 using namespace muon;
0009
0010 L1Analysis::L1AnalysisRecoMuon2::L1AnalysisRecoMuon2(const edm::ParameterSet& pset, edm::ConsumesCollector iC)
0011     : muPropagatorSetup1st_(pset.getParameter<edm::ParameterSet>("muProp1st"), iC),
0012       muPropagatorSetup2nd_(pset.getParameter<edm::ParameterSet>("muProp2nd"), iC) {}
0013
0014 L1Analysis::L1AnalysisRecoMuon2::~L1AnalysisRecoMuon2() {}
0015
0016 void L1Analysis::L1AnalysisRecoMuon2::SetMuon(const edm::Event& event,
0017                                               const edm::EventSetup& setup,
0018                                               edm::Handle<reco::MuonCollection> muons,
0019                                               edm::Handle<reco::VertexCollection> vertices,
0020                                               double METx,
0021                                               double METy,
0022                                               unsigned maxMuon) {
0023   recoMuon_.nMuons = 0;
0024
0025   for (reco::MuonCollection::const_iterator it = muons->begin(); it != muons->end() && recoMuon_.nMuons < maxMuon;
0026        ++it) {
0027     recoMuon_.e.push_back(it->energy());
0028     recoMuon_.pt.push_back(it->pt());
0029     recoMuon_.et.push_back(it->et());
0030     recoMuon_.eta.push_back(it->eta());
0031     recoMuon_.phi.push_back(it->phi());
0032     recoMuon_.charge.push_back(it->charge());
0033
0034     //check isLooseMuon
0035     bool flagLoose = isLooseMuonCustom(*it);
0036     recoMuon_.isLooseMuon.push_back(flagLoose);
0037
0038     //check isMediumMuon
0039     bool flagMedium = isMediumMuonCustom(*it);
0040     recoMuon_.isMediumMuon.push_back(flagMedium);
0041
0042     //check isTightMuon
0043     bool flagTight = false;
0044     if (vertices.isValid())
0045       flagTight = isTightMuonCustom(*it, (*vertices)[0]);
0046     recoMuon_.isTightMuon.push_back(flagTight);
0047
0048     double iso = (it->pfIsolationR04().sumChargedHadronPt +
0049                   max(0.,
0050                       it->pfIsolationR04().sumNeutralHadronEt + it->pfIsolationR04().sumPhotonEt -
0051                           0.5 * it->pfIsolationR04().sumPUPt)) /
0052                  it->pt();
0053     recoMuon_.iso.push_back(iso);
0054
0055     double MET_local = TMath::Sqrt(METx * METx + METy * METy);
0056     recoMuon_.met.push_back(MET_local);
0057
0058     TLorentzVector METP4;
0059     METP4.SetPxPyPzE(METx, METy, 0, MET_local);
0060
0061     TLorentzVector Muon;
0062     Muon.SetPtEtaPhiE(it->pt(), it->eta(), it->phi(), it->energy());
0063
0064     double scalSum = MET_local + Muon.Pt();
0065     TLorentzVector vecSum(Muon);
0066     vecSum += METP4;
0067     double vecSumPt = vecSum.Pt();
0068
0069     recoMuon_.mt.push_back(TMath::Sqrt(scalSum * scalSum - vecSumPt * vecSumPt));
0070
0071     recoMuon_.nMuons++;
0072
0073     // extrapolation of track coordinates
0074     TrajectoryStateOnSurface stateAtMuSt1 = muPropagator1st_.extrapolate(*it);
0075     if (stateAtMuSt1.isValid()) {
0076       recoMuon_.etaSt1.push_back(stateAtMuSt1.globalPosition().eta());
0077       recoMuon_.phiSt1.push_back(stateAtMuSt1.globalPosition().phi());
0078     } else {
0079       recoMuon_.etaSt1.push_back(-9999);
0080       recoMuon_.phiSt1.push_back(-9999);
0081     }
0082
0083     TrajectoryStateOnSurface stateAtMuSt2 = muPropagator2nd_.extrapolate(*it);
0084     if (stateAtMuSt2.isValid()) {
0085       recoMuon_.etaSt2.push_back(stateAtMuSt2.globalPosition().eta());
0086       recoMuon_.phiSt2.push_back(stateAtMuSt2.globalPosition().phi());
0087     } else {
0088       recoMuon_.etaSt2.push_back(-9999);
0089       recoMuon_.phiSt2.push_back(-9999);
0090     }
0091   }
0092 }
0093
0094 void L1Analysis::L1AnalysisRecoMuon2::init(const edm::EventSetup& eventSetup) {
0095   muPropagator1st_ = muPropagatorSetup1st_.init(eventSetup);
0096   muPropagator2nd_ = muPropagatorSetup2nd_.init(eventSetup);
0097 }