Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:20:29

0001 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0002 
0003 #include "DataFormats/Common/interface/Handle.h"
0004 
0005 #include "DataFormats/Common/interface/Ref.h"
0006 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0007 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0008 
0009 #include "DataFormats/JetReco/interface/CaloJet.h"
0010 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0011 
0012 #include "DataFormats/MuonReco/interface/Muon.h"
0013 
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 
0017 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0019 #include "FWCore/Utilities/interface/InputTag.h"
0020 
0021 #include "HLTrigger/JetMET/interface/HLTRFilter.h"
0022 #include "DataFormats/Common/interface/Ref.h"
0023 
0024 //
0025 // constructors and destructor
0026 //
0027 HLTRFilter::HLTRFilter(const edm::ParameterSet& iConfig)
0028     : HLTFilter(iConfig),
0029       inputTag_(iConfig.getParameter<edm::InputTag>("inputTag")),
0030       inputMetTag_(iConfig.getParameter<edm::InputTag>("inputMetTag")),
0031       doMuonCorrection_(iConfig.getParameter<bool>("doMuonCorrection")),
0032       min_R_(iConfig.getParameter<double>("minR")),
0033       min_MR_(iConfig.getParameter<double>("minMR")),
0034       DoRPrime_(iConfig.getParameter<bool>("doRPrime")),
0035       accept_NJ_(iConfig.getParameter<bool>("acceptNJ")),
0036       R_offset_((iConfig.existsAs<double>("R2Offset") ? iConfig.getParameter<double>("R2Offset") : 0)),
0037       MR_offset_((iConfig.existsAs<double>("MROffset") ? iConfig.getParameter<double>("MROffset") : 0)),
0038       R_MR_cut_((iConfig.existsAs<double>("RMRCut") ? iConfig.getParameter<double>("RMRCut") : -999999.))
0039 
0040 {
0041   m_theInputToken = consumes<std::vector<math::XYZTLorentzVector>>(inputTag_);
0042   m_theMETToken = consumes<edm::View<reco::MET>>(inputMetTag_);
0043   LogDebug("") << "Inputs/minR/minMR/doRPrime/acceptNJ/R2Offset/MROffset/RMRCut : " << inputTag_.encode() << " "
0044                << inputMetTag_.encode() << " " << min_R_ << " " << min_MR_ << " " << accept_NJ_ << R_offset_ << " "
0045                << MR_offset_ << " " << R_MR_cut_ << ".";
0046 
0047   //put a dummy METCollection into the event, holding values for MR and Rsq
0048   produces<reco::METCollection>();
0049 }
0050 
0051 HLTRFilter::~HLTRFilter() = default;
0052 
0053 void HLTRFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0054   edm::ParameterSetDescription desc;
0055   makeHLTFilterDescription(desc);
0056   desc.add<edm::InputTag>("inputTag", edm::InputTag("hltRHemisphere"));
0057   desc.add<edm::InputTag>("inputMetTag", edm::InputTag("hltMet"));
0058   desc.add<bool>("doMuonCorrection", false);
0059   desc.add<double>("minR", 0.3);
0060   desc.add<double>("minMR", 100.0);
0061   desc.add<bool>("doRPrime", false);
0062   desc.add<bool>("acceptNJ", true);
0063   desc.add<double>("R2Offset", 0.0);
0064   desc.add<double>("MROffset", 0.0);
0065   desc.add<double>("RMRCut", -999999.0);
0066   descriptions.add("hltRFilter", desc);
0067 }
0068 
0069 //
0070 // member functions
0071 //
0072 
0073 // ------------ method called to produce the data  ------------
0074 bool HLTRFilter::hltFilter(edm::Event& iEvent,
0075                            const edm::EventSetup& iSetup,
0076                            trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0077   using namespace std;
0078   using namespace edm;
0079   using namespace reco;
0080 
0081   // get hold of collection of objects
0082   Handle<vector<math::XYZTLorentzVector>> hemispheres;
0083   iEvent.getByToken(m_theInputToken, hemispheres);
0084 
0085   // get hold of the MET Collection
0086   edm::Handle<edm::View<reco::MET>> inputMet;
0087   iEvent.getByToken(m_theMETToken, inputMet);
0088 
0089   // check the the input collections are available
0090   if (not hemispheres.isValid() or not inputMet.isValid())
0091     return false;
0092 
0093   if (hemispheres
0094           ->empty()) {  // the Hemisphere Maker will produce an empty collection of hemispheres if the number of jets in the
0095     return accept_NJ_;  // event is greater than the maximum number of jets
0096   }
0097 
0098   //***********************************
0099   //Calculate R
0100 
0101   int nMuons;
0102   switch (hemispheres->size()) {
0103     case 2:
0104       nMuons = 0;
0105       break;
0106     case 5:
0107       nMuons = 1;
0108       break;
0109     case 10:
0110       nMuons = 2;
0111       break;
0112     default:
0113       return false;  //invalid hemisphere collection
0114   }
0115 
0116   //muons as MET
0117   TLorentzVector ja(hemispheres->at(0).x(), hemispheres->at(0).y(), hemispheres->at(0).z(), hemispheres->at(0).t());
0118   TLorentzVector jb(hemispheres->at(1).x(), hemispheres->at(1).y(), hemispheres->at(1).z(), hemispheres->at(1).t());
0119 
0120   std::vector<math::XYZTLorentzVector> muonVec;
0121 
0122   double MR = CalcMR(ja, jb);
0123   double R = CalcR(MR, ja, jb, inputMet, muonVec);
0124 
0125   if (MR >= min_MR_ && R >= min_R_ && ((R * R - R_offset_) * (MR - MR_offset_)) >= R_MR_cut_) {
0126     addObjects(iEvent, filterproduct, MR, R * R);
0127     return true;
0128   }
0129   if (nMuons == 0) {
0130     addObjects(iEvent, filterproduct, MR, R * R);
0131     return false;  // if no muons and we get here, reject event
0132   }
0133   //Lead Muon as Jet
0134   ja.SetXYZT(hemispheres->at(3).x(), hemispheres->at(3).y(), hemispheres->at(3).z(), hemispheres->at(3).t());
0135   jb.SetXYZT(hemispheres->at(4).x(), hemispheres->at(4).y(), hemispheres->at(4).z(), hemispheres->at(4).t());
0136   muonVec.push_back(hemispheres->at(2));  // lead muon at position 2
0137 
0138   MR = CalcMR(ja, jb);
0139   R = CalcR(MR, ja, jb, inputMet, muonVec);
0140   if (MR >= min_MR_ && R >= min_R_ && ((R * R - R_offset_) * (MR - MR_offset_)) >= R_MR_cut_) {
0141     addObjects(iEvent, filterproduct, MR, R * R);
0142     return true;
0143   }
0144 
0145   if (nMuons == 1) {
0146     addObjects(iEvent, filterproduct, MR, R * R);
0147     return false;  // if no muons and we get here, reject event
0148   }
0149 
0150   muonVec.pop_back();
0151   //Second Muon as Jet
0152   ja.SetXYZT(hemispheres->at(6).x(), hemispheres->at(6).y(), hemispheres->at(6).z(), hemispheres->at(6).t());
0153   jb.SetXYZT(hemispheres->at(7).x(), hemispheres->at(7).y(), hemispheres->at(7).z(), hemispheres->at(7).t());
0154   muonVec.push_back(hemispheres->at(5));  // sublead muon at position 5
0155 
0156   MR = CalcMR(ja, jb);
0157   R = CalcR(MR, ja, jb, inputMet, muonVec);
0158   if (MR >= min_MR_ && R >= min_R_ && ((R * R - R_offset_) * (MR - MR_offset_)) >= R_MR_cut_) {
0159     addObjects(iEvent, filterproduct, MR, R * R);
0160     return true;
0161   }
0162 
0163   ja.SetXYZT(hemispheres->at(8).x(), hemispheres->at(8).y(), hemispheres->at(8).z(), hemispheres->at(8).t());
0164   jb.SetXYZT(hemispheres->at(9).x(), hemispheres->at(9).y(), hemispheres->at(9).z(), hemispheres->at(9).t());
0165   muonVec.push_back(hemispheres->at(2));  // lead muon at position 2
0166 
0167   MR = CalcMR(ja, jb);
0168   R = CalcR(MR, ja, jb, inputMet, muonVec);
0169   if (MR >= min_MR_ && R >= min_R_ && ((R * R - R_offset_) * (MR - MR_offset_)) >= R_MR_cut_) {
0170     addObjects(iEvent, filterproduct, MR, R * R);
0171     return true;
0172   }
0173 
0174   // filter decision
0175   addObjects(iEvent, filterproduct, MR, R * R);
0176   return false;
0177 }
0178 
0179 void HLTRFilter::addObjects(edm::Event& iEvent,
0180                             trigger::TriggerFilterObjectWithRefs& filterproduct,
0181                             double MR,
0182                             double Rsq) const {
0183   //create METCollection for storing MR and Rsq results
0184   std::unique_ptr<reco::METCollection> razorObject(new reco::METCollection());
0185 
0186   reco::MET::LorentzVector mrRsqP4(MR, Rsq, 0, 0);
0187   reco::MET::Point vtx(0, 0, 0);
0188 
0189   reco::MET mrRsq(mrRsqP4, vtx);
0190   razorObject->push_back(mrRsq);
0191 
0192   edm::RefProd<reco::METCollection> ref_before_put = iEvent.getRefBeforePut<reco::METCollection>();
0193 
0194   //put the METCollection into the event (necessary because of how addCollectionTag works...)
0195   iEvent.put(std::move(razorObject));
0196   edm::Ref<reco::METCollection> mrRsqRef(ref_before_put, 0);
0197 
0198   //add it to the trigger object collection
0199   if (saveTags())
0200     filterproduct.addCollectionTag(edm::InputTag(*moduleLabel()));
0201   filterproduct.addObject(trigger::TriggerMET, mrRsqRef);  //give it the ID of a MET object
0202 }
0203 
0204 double HLTRFilter::CalcMR(TLorentzVector ja, TLorentzVector jb) {
0205   if (ja.Pt() <= 0.1)
0206     return -1;
0207 
0208   ja.SetPtEtaPhiM(ja.Pt(), ja.Eta(), ja.Phi(), 0.0);
0209   jb.SetPtEtaPhiM(jb.Pt(), jb.Eta(), jb.Phi(), 0.0);
0210 
0211   if (ja.Pt() > jb.Pt()) {
0212     TLorentzVector temp = ja;
0213     ja = jb;
0214     jb = temp;
0215   }
0216 
0217   double A = ja.P();
0218   double B = jb.P();
0219   double az = ja.Pz();
0220   double bz = jb.Pz();
0221   TVector3 jaT, jbT;
0222   jaT.SetXYZ(ja.Px(), ja.Py(), 0.0);
0223   jbT.SetXYZ(jb.Px(), jb.Py(), 0.0);
0224   double ATBT = (jaT + jbT).Mag2();
0225 
0226   double MR = sqrt((A + B) * (A + B) - (az + bz) * (az + bz) -
0227                    (jbT.Dot(jbT) - jaT.Dot(jaT)) * (jbT.Dot(jbT) - jaT.Dot(jaT)) / (jaT + jbT).Mag2());
0228 
0229   double mybeta = (jbT.Dot(jbT) - jaT.Dot(jaT)) / sqrt(ATBT * ((A + B) * (A + B) - (az + bz) * (az + bz)));
0230 
0231   double mygamma = 1. / sqrt(1. - mybeta * mybeta);
0232 
0233   //use gamma times MRstar
0234   return MR * mygamma;
0235 }
0236 
0237 double HLTRFilter::CalcR(double MR,
0238                          TLorentzVector ja,
0239                          TLorentzVector jb,
0240                          edm::Handle<edm::View<reco::MET>> inputMet,
0241                          const std::vector<math::XYZTLorentzVector>& muons) {
0242   //now we can calculate MTR
0243   TVector3 met;
0244   met.SetPtEtaPhi((inputMet->front()).pt(), 0.0, (inputMet->front()).phi());
0245 
0246   std::vector<math::XYZTLorentzVector>::const_iterator muonIt;
0247   for (muonIt = muons.begin(); muonIt != muons.end(); muonIt++) {
0248     TVector3 tmp;
0249     tmp.SetPtEtaPhi(muonIt->pt(), 0, muonIt->phi());
0250     met -= tmp;
0251   }
0252 
0253   double MTR = sqrt(0.5 * (met.Mag() * (ja.Pt() + jb.Pt()) - met.Dot(ja.Vect() + jb.Vect())));
0254 
0255   //filter events
0256   return float(MTR) / float(MR);  //R
0257 }
0258 DEFINE_FWK_MODULE(HLTRFilter);