JetMET/plugins/HLTRHemisphere.cc

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 // constructors and destructor
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   //register your products
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 // member functions
0064 //
0065
0066 // ------------ method called to produce the data  ------------
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   // get hold of collection of objects
0077   //   Handle<CaloJetCollection> jets;
0078   Handle<View<Jet>> jets;
0079   iEvent.getByToken(m_theJetToken, jets);
0080
0081   // get hold of the muons, if necessary
0082   Handle<vector<reco::RecoChargedCandidate>> muons;
0083   if (doMuonCorrection_)
0084     iEvent.getByToken(m_theMuonToken, muons);
0085
0086   // The output Collection
0087   std::unique_ptr<vector<math::XYZTLorentzVector>> Hemispheres(new vector<math::XYZTLorentzVector>);
0088
0089   // look at all objects, check cuts and add to filter object
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_;  // too many jets, accept for timing
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;                            // skip muons out of eta range or too low pT
0113       if (nPassMu >= 2) {                    // if we have already accepted two muons, accept the event
0114         iEvent.put(std::move(Hemispheres));  // too many muons, accept for timing
0115         return true;
0116       }
0117       muonIndex[nPassMu++] = index;
0118     }
0119     //muons as MET
0120     this->ComputeHemispheres(Hemispheres, JETS);
0121     //lead muon as jet
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);  // lead muon as jet
0128       if (nPassMu > 1) {                                       // two passing muons
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);  // lead muon as v, second muon as jet
0134         muonJets.push_back(leadMu.p4());
0135         this->ComputeHemispheres(Hemispheres, JETS, &muonJets);  // both muon as jets
0136       }
0137     }
0138   } else {  // do MuonCorrection==false
0139     if (n < 2)
0140       return false;                               // not enough jets and not adding in muons
0141     this->ComputeHemispheres(Hemispheres, JETS);  // don't do the muon isolation, just run once and done
0142   }
0143   //Format:
0144   // 0 muon: 2 hemispheres (2)
0145   // 1 muon: 2 hemisheress + leadMuP4 + 2 hemispheres (5)
0146   // 2 muon: 2 hemispheres + leadMuP4 + 2 hemispheres + 2ndMuP4 + 4 Hemispheres (10)
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) {  // put empty hemispheres if not enough jets
0163     hlist->push_back(j1R);
0164     hlist->push_back(j2R);
0165     return;
0166   }
0167   unsigned int N_comb = pow(2, nJets);  // compute the number of combinations of jets possible
0168   //Make the hemispheres
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);  // to silence LLVM analyzer warning
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);  // to silence LLVM analyzer warning
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);