File indexing completed on 2024-10-29 06:08:28
0001 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0002 #include "FWCore/Framework/interface/ESHandle.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004
0005 #include "DataFormats/Common/interface/Handle.h"
0006 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0007
0008 #include "DataFormats/Common/interface/Ref.h"
0009
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 #include "TVector3.h"
0014 #include "TLorentzVector.h"
0015
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "FWCore/Utilities/interface/InputTag.h"
0019
0020 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0021
0022 #include "HLTrigger/JetMET/interface/HLTRHemisphere.h"
0023
0024 #include <vector>
0025
0026
0027
0028
0029 HLTRHemisphere::HLTRHemisphere(const edm::ParameterSet& iConfig)
0030 : inputTag_(iConfig.getParameter<edm::InputTag>("inputTag")),
0031 muonTag_(iConfig.getParameter<edm::InputTag>("muonTag")),
0032 doMuonCorrection_(iConfig.getParameter<bool>("doMuonCorrection")),
0033 muonEta_(iConfig.getParameter<double>("maxMuonEta")),
0034 min_Jet_Pt_(iConfig.getParameter<double>("minJetPt")),
0035 max_Eta_(iConfig.getParameter<double>("maxEta")),
0036 max_NJ_(iConfig.getParameter<int>("maxNJ")),
0037 accNJJets_(iConfig.getParameter<bool>("acceptNJ")) {
0038 LogDebug("") << "Input/minJetPt/maxEta/maxNJ/acceptNJ : " << inputTag_.encode() << " " << min_Jet_Pt_ << "/"
0039 << max_Eta_ << "/" << max_NJ_ << "/" << accNJJets_ << ".";
0040
0041 m_theJetToken = consumes<edm::View<reco::Jet>>(inputTag_);
0042 m_theMuonToken = consumes<std::vector<reco::RecoChargedCandidate>>(muonTag_);
0043
0044 produces<std::vector<math::XYZTLorentzVector>>();
0045 }
0046
0047 HLTRHemisphere::~HLTRHemisphere() = default;
0048
0049 void HLTRHemisphere::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0050 edm::ParameterSetDescription desc;
0051 desc.add<edm::InputTag>("inputTag", edm::InputTag("hltMCJetCorJetIcone5HF07"));
0052 desc.add<edm::InputTag>("muonTag", edm::InputTag(""));
0053 desc.add<bool>("doMuonCorrection", false);
0054 desc.add<double>("maxMuonEta", 2.1);
0055 desc.add<double>("minJetPt", 30.0);
0056 desc.add<double>("maxEta", 3.0);
0057 desc.add<int>("maxNJ", 7);
0058 desc.add<bool>("acceptNJ", true);
0059 descriptions.add("hltRHemisphere", desc);
0060 }
0061
0062
0063
0064
0065
0066
0067 bool HLTRHemisphere::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0068 using namespace std;
0069 using namespace edm;
0070 using namespace reco;
0071 using namespace math;
0072 using namespace trigger;
0073
0074 typedef XYZTLorentzVector LorentzVector;
0075
0076
0077
0078 Handle<View<Jet>> jets;
0079 iEvent.getByToken(m_theJetToken, jets);
0080
0081
0082 Handle<vector<reco::RecoChargedCandidate>> muons;
0083 if (doMuonCorrection_)
0084 iEvent.getByToken(m_theMuonToken, muons);
0085
0086
0087 std::unique_ptr<vector<math::XYZTLorentzVector>> Hemispheres(new vector<math::XYZTLorentzVector>);
0088
0089
0090 int n(0);
0091 vector<math::XYZTLorentzVector> JETS;
0092 for (auto const& i : *jets) {
0093 if (std::abs(i.eta()) < max_Eta_ && i.pt() >= min_Jet_Pt_) {
0094 JETS.push_back(i.p4());
0095 n++;
0096 }
0097 }
0098
0099 if (n > max_NJ_ && max_NJ_ != -1) {
0100 iEvent.put(std::move(Hemispheres));
0101 return accNJJets_;
0102 }
0103
0104 if (doMuonCorrection_) {
0105 const int nMu = 2;
0106 int muonIndex[nMu] = {-1, -1};
0107 std::vector<reco::RecoChargedCandidate>::const_iterator muonIt;
0108 int index = 0;
0109 int nPassMu = 0;
0110 for (muonIt = muons->begin(); muonIt != muons->end(); muonIt++, index++) {
0111 if (std::abs(muonIt->eta()) > muonEta_ || muonIt->pt() < min_Jet_Pt_)
0112 continue;
0113 if (nPassMu >= 2) {
0114 iEvent.put(std::move(Hemispheres));
0115 return true;
0116 }
0117 muonIndex[nPassMu++] = index;
0118 }
0119
0120 this->ComputeHemispheres(Hemispheres, JETS);
0121
0122 if (nPassMu > 0) {
0123 std::vector<math::XYZTLorentzVector> muonJets;
0124 reco::RecoChargedCandidate leadMu = muons->at(muonIndex[0]);
0125 muonJets.push_back(leadMu.p4());
0126 Hemispheres->push_back(leadMu.p4());
0127 this->ComputeHemispheres(Hemispheres, JETS, &muonJets);
0128 if (nPassMu > 1) {
0129 muonJets.pop_back();
0130 reco::RecoChargedCandidate secondMu = muons->at(muonIndex[1]);
0131 muonJets.push_back(secondMu.p4());
0132 Hemispheres->push_back(secondMu.p4());
0133 this->ComputeHemispheres(Hemispheres, JETS, &muonJets);
0134 muonJets.push_back(leadMu.p4());
0135 this->ComputeHemispheres(Hemispheres, JETS, &muonJets);
0136 }
0137 }
0138 } else {
0139 if (n < 2)
0140 return false;
0141 this->ComputeHemispheres(Hemispheres, JETS);
0142 }
0143
0144
0145
0146
0147 iEvent.put(std::move(Hemispheres));
0148 return true;
0149 }
0150
0151 void HLTRHemisphere::ComputeHemispheres(std::unique_ptr<std::vector<math::XYZTLorentzVector>>& hlist,
0152 const std::vector<math::XYZTLorentzVector>& JETS,
0153 std::vector<math::XYZTLorentzVector>* extraJets) {
0154 using namespace math;
0155 using namespace reco;
0156 XYZTLorentzVector j1R(0.1, 0., 0., 0.1);
0157 XYZTLorentzVector j2R(0.1, 0., 0., 0.1);
0158 int nJets = JETS.size();
0159 if (extraJets)
0160 nJets += extraJets->size();
0161
0162 if (nJets < 2) {
0163 hlist->push_back(j1R);
0164 hlist->push_back(j2R);
0165 return;
0166 }
0167 unsigned int N_comb = pow(2, nJets);
0168
0169 double M_minR = 9999999999.0;
0170 unsigned int j_count;
0171 for (unsigned int i = 0; i < N_comb; i++) {
0172 XYZTLorentzVector j_temp1, j_temp2;
0173 unsigned int itemp = i;
0174 j_count = N_comb / 2;
0175 unsigned int count = 0;
0176 while (j_count > 0) {
0177 if (itemp / j_count == 1) {
0178 if (count < JETS.size())
0179 j_temp1 += JETS.at(count);
0180 else {
0181 assert(extraJets);
0182 j_temp1 += extraJets->at(count - JETS.size());
0183 }
0184 } else {
0185 if (count < JETS.size())
0186 j_temp2 += JETS.at(count);
0187 else {
0188 assert(extraJets);
0189 j_temp2 += extraJets->at(count - JETS.size());
0190 }
0191 }
0192 itemp -= j_count * (itemp / j_count);
0193 j_count /= 2;
0194 count++;
0195 }
0196 double M_temp = j_temp1.M2() + j_temp2.M2();
0197 if (M_temp < M_minR) {
0198 M_minR = M_temp;
0199 j1R = j_temp1;
0200 j2R = j_temp2;
0201 }
0202 }
0203
0204 hlist->push_back(j1R);
0205 hlist->push_back(j2R);
0206 return;
0207 }
0208
0209 DEFINE_FWK_MODULE(HLTRHemisphere);