JetMET/src/HLTAlphaTFilter.cc

0001 /** \class HLTAlphaTFilter
0002  *
0003  *
0004  *  \author Bryn Mathias
0005  *  \modified Mark Baber, Adam Elwood
0006  *
0007  *  Filter for the AlphaT SUSY analysis
0008  *  Makes a trigger decision based on the event HT
0009  *  and the AlphaT value. AlphaT cut is chosen
0010  *  to reject all the QCD background
0011  *
0012  */
0013
0014 #include <vector>
0015
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "DataFormats/Math/interface/deltaPhi.h"
0023 #include "DataFormats/Common/interface/Handle.h"
0024 #include "HLTrigger/JetMET/interface/HLTAlphaTFilter.h"
0025 #include "HLTrigger/JetMET/interface/AlphaT.h"
0026 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0027
0028 typedef ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double>> LorentzV;
0029
0030 //
0031 // constructors and destructor
0032 //
0033 template <typename T>
0034 HLTAlphaTFilter<T>::HLTAlphaTFilter(const edm::ParameterSet& iConfig) : HLTFilter(iConfig) {
0035   inputJetTag_ = iConfig.getParameter<edm::InputTag>("inputJetTag");
0036   inputJetTagFastJet_ = iConfig.getParameter<edm::InputTag>("inputJetTagFastJet");
0037   minPtJet_ = iConfig.getParameter<std::vector<double>>("minPtJet");
0038   etaJet_ = iConfig.getParameter<std::vector<double>>("etaJet");
0039   maxNJets_ = iConfig.getParameter<unsigned int>("maxNJets");
0040   minHt_ = iConfig.getParameter<double>("minHt");
0041   minAlphaT_ = iConfig.getParameter<double>("minAlphaT");
0042   triggerType_ = iConfig.getParameter<int>("triggerType");
0043   dynamicAlphaT_ = iConfig.getParameter<bool>("dynamicAlphaT");
0044   setDHtZero_ = iConfig.getParameter<bool>("setDHtZero");
0045   // sanity checks
0046
0047   if ((minPtJet_.size() != etaJet_.size()) || ((minPtJet_.empty()) || (etaJet_.empty())) ||
0048       (((minPtJet_.size() < 2) || (etaJet_.size() < 2)))) {
0049     edm::LogError("HLTAlphaTFilter") << "inconsistent module configuration!";
0050   }
0051
0052   //register your products
0053   m_theRecoJetToken = consumes<std::vector<T>>(inputJetTag_);
0054   m_theFastJetToken = consumes<std::vector<T>>(inputJetTagFastJet_);
0055 }
0056
0057 template <typename T>
0058 HLTAlphaTFilter<T>::~HLTAlphaTFilter() = default;
0059
0060 template <typename T>
0061 void HLTAlphaTFilter<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0062   edm::ParameterSetDescription desc;
0063   makeHLTFilterDescription(desc);
0064   desc.add<edm::InputTag>("inputJetTag", edm::InputTag("hltMCJetCorJetIcone5HF07"));
0065   desc.add<edm::InputTag>("inputJetTagFastJet", edm::InputTag("hltMCJetCorJetIcone5HF07"));
0066
0067   {
0068     std::vector<double> temp1;
0069     temp1.reserve(2);
0070     temp1.push_back(20.0);
0071     temp1.push_back(20.0);
0072     desc.add<std::vector<double>>("minPtJet", temp1);
0073   }
0074   desc.add<int>("minNJet", 0);
0075   {
0076     std::vector<double> temp1;
0077     temp1.reserve(2);
0078     temp1.push_back(9999.0);
0079     temp1.push_back(9999.0);
0080     desc.add<std::vector<double>>("etaJet", temp1);
0081   }
0082   desc.add<unsigned int>("maxNJets", 32);
0083   desc.add<double>("minHt", 0.0);
0084   desc.add<double>("minAlphaT", 0.0);
0085   desc.add<int>("triggerType", trigger::TriggerJet);
0086   desc.add<bool>("dynamicAlphaT", true);  //Set to reproduce old behaviour
0087   desc.add<bool>("setDHtZero", false);    //Set to reproduce old behaviour
0088   descriptions.add(defaultModuleLabel<HLTAlphaTFilter<T>>(), desc);
0089 }
0090
0091 // ------------ method called to produce the data  ------------
0092 template <typename T>
0093 bool HLTAlphaTFilter<T>::hltFilter(edm::Event& iEvent,
0094                                    const edm::EventSetup& iSetup,
0095                                    trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0096   using namespace std;
0097   using namespace edm;
0098   using namespace reco;
0099   using namespace trigger;
0100
0101   typedef vector<T> TCollection;
0102   typedef Ref<TCollection> TRef;
0103
0104   // The filter object
0105   if (saveTags())
0106     filterproduct.addCollectionTag(inputJetTag_);
0107
0108   TRef ref;
0109   // Get the Candidates
0110   Handle<TCollection> recojets;
0111   iEvent.getByToken(m_theRecoJetToken, recojets);
0112
0113   // We have to also look at the L1 FastJet Corrections, at the same time we look at our other jets.
0114   // We calcualte our HT from the FastJet collection and AlphaT from the standard collection.
0115   CaloJetRef ref_FastJet;
0116   // Get the Candidates
0117   Handle<TCollection> recojetsFastJet;
0118   iEvent.getByToken(m_theFastJetToken, recojetsFastJet);
0119
0120   int n(0);
0121
0122   //OLD - DYNAMIC ALPHAT BEHAVIOUR
0123   //Produce AlphaT dynamically, first on two jets, then three etc until it does
0124   //or doesn't pass
0125   if (dynamicAlphaT_) {
0126     // look at all candidates,  check cuts and add to filter object
0127     int flag(0);
0128     double htFast = 0.;
0129     double aT = 0.;
0130     unsigned int njets(0);
0131
0132     if (recojets->size() > 1) {
0133       // events with at least two jets, needed for alphaT
0134       // Make a vector of Lorentz Jets for the AlphaT calcualtion
0135       std::vector<LorentzV> jets;
0136       typename TCollection::const_iterator ijet = recojets->begin();
0137       typename TCollection::const_iterator ijetFast = recojetsFastJet->begin();
0138       typename TCollection::const_iterator jjet = recojets->end();
0139
0140       for (; ijet != jjet; ijet++, ijetFast++) {
0141         if (flag == 1)
0142           break;
0143         // Do Some Jet selection!
0144         if (std::abs(ijet->eta()) > etaJet_.at(0))
0145           continue;
0146         if (ijet->et() < minPtJet_.at(0))
0147           continue;
0148         njets++;
0149
0150         if (njets >
0151             maxNJets_)  //to keep timing reasonable - if too many jets passing pt / eta cuts, just accept the event
0152           flag = 1;
0153
0154         else {
0155           if (std::abs(ijetFast->eta()) < etaJet_.at(1)) {
0156             if (ijetFast->et() > minPtJet_.at(1)) {
0157               // Add to HT
0158               htFast += ijetFast->et();
0159             }
0160           }
0161
0162           // Add to JetVector
0163           LorentzV JetLVec(ijet->pt(), ijet->eta(), ijet->phi(), ijet->mass());
0164           jets.push_back(JetLVec);
0165           aT = AlphaT(jets, setDHtZero_).value();
0166           if (htFast > minHt_ && aT > minAlphaT_) {
0167             // set flat to one so that we don't carry on looping though the jets
0168             flag = 1;
0169           }
0170         }
0171       }
0172
0173       //If passed, add jets to the filter product for the DQM
0174       if (flag == 1) {
0175         for (typename TCollection::const_iterator recojet = recojets->begin(); recojet != jjet; recojet++) {
0176           if (recojet->et() > minPtJet_.at(0)) {
0177             ref = TRef(recojets, distance(recojets->begin(), recojet));
0178             filterproduct.addObject(triggerType_, ref);
0179             n++;
0180           }
0181         }
0182       }
0183     }  // events with at least two jet
0184
0185     // filter decision
0186     bool accept(n > 0);
0187
0188     return accept;
0189   }
0190   // NEW - STATIC ALPHAT BEHAVIOUR
0191   // just reproduce
0192   else {
0193     // look at all candidates,  check cuts and add to filter object
0194     int flag(0);
0195     typename TCollection::const_iterator ijet = recojets->begin();
0196     typename TCollection::const_iterator jjet = recojets->end();
0197
0198     if ((recojets->size() > 1) && (recojetsFastJet->size() > 1)) {
0199       // events with at least two jets, needed for alphaT
0200       double htFast = 0.;
0201       typename TCollection::const_iterator ijetFast = recojetsFastJet->begin();
0202       typename TCollection::const_iterator jjetFast = recojetsFastJet->end();
0203       for (; ijetFast != jjetFast; ijetFast++) {
0204         if (std::abs(ijetFast->eta()) < etaJet_.at(1)) {
0205           if (ijetFast->et() > minPtJet_.at(1)) {
0206             // Add to HT
0207             htFast += ijetFast->et();
0208           }
0209         }
0210       }
0211
0212       if (htFast > minHt_) {
0213         unsigned int njets(0);
0214         // Make a vector of Lorentz Jets for the AlphaT calcualtion
0215         std::vector<LorentzV> jets;
0216         for (; ijet != jjet; ijet++) {
0217           if (std::abs(ijet->eta()) > etaJet_.at(0))
0218             continue;
0219           if (ijet->et() < minPtJet_.at(0))
0220             continue;
0221           njets++;
0222
0223           if (njets >
0224               maxNJets_) {  //to keep timing reasonable - if too many jets passing pt / eta cuts, just accept the event
0225             flag = 1;
0226             break;  //Added for efficiency
0227           }
0228
0229           // Add to JetVector
0230           LorentzV JetLVec(ijet->pt(), ijet->eta(), ijet->phi(), ijet->mass());
0231           jets.push_back(JetLVec);
0232         }
0233
0234         if (flag != 1) {  //Added for efficiency
0235           //Calculate the value for alphaT
0236           float aT = AlphaT(jets, setDHtZero_).value();
0237
0238           // Trigger decision!
0239           if (aT > minAlphaT_) {
0240             flag = 1;
0241           }
0242         }
0243
0244         //If passed, add the jets to the filterproduct for DQM
0245         if (flag == 1) {
0246           for (typename TCollection::const_iterator recojet = recojets->begin(); recojet != jjet; recojet++) {
0247             if (recojet->et() > minPtJet_.at(0)) {
0248               ref = TRef(recojets, distance(recojets->begin(), recojet));
0249               filterproduct.addObject(triggerType_, ref);
0250               n++;
0251             }
0252           }
0253         }
0254       }
0255     }  // events with at least two jet
0256
0257     // filter decision
0258     bool accept(n > 0);
0259
0260     return accept;
0261   }
0262
0263   // THIS WILL NEVER HAPPEN
0264   return true;
0265 }