Egamma/plugins/HLTEgammaDoubleLegCombFilter.cc

0001 /** \class HLTEgammaDoubleLegCombFilter
0002  *
0003  * $Id: HLTEgammaDoubleLegCombFilter.cc,
0004  *
0005  *  \author Sam Harper (RAL)
0006  *
0007  */
0008
0009 #include "HLTEgammaDoubleLegCombFilter.h"
0010
0011 #include "DataFormats/Common/interface/Handle.h"
0012
0013 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0018 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0019 #include "DataFormats/Math/interface/deltaR.h"
0020
0021 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0022
0023 //
0024 // constructors and destructor
0025 //
0026 HLTEgammaDoubleLegCombFilter::HLTEgammaDoubleLegCombFilter(const edm::ParameterSet& iConfig) : HLTFilter(iConfig) {
0027   firstLegLastFilterTag_ = iConfig.getParameter<edm::InputTag>("firstLegLastFilter");
0028   secondLegLastFilterTag_ = iConfig.getParameter<edm::InputTag>("secondLegLastFilter");
0029   nrRequiredFirstLeg_ = iConfig.getParameter<int>("nrRequiredFirstLeg");
0030   nrRequiredSecondLeg_ = iConfig.getParameter<int>("nrRequiredSecondLeg");
0031   nrRequiredUniqueLeg_ = iConfig.getParameter<int>("nrRequiredUniqueLeg");
0032   maxMatchDR_ = iConfig.getParameter<double>("maxMatchDR");
0033   firstLegLastFilterToken_ = consumes<trigger::TriggerFilterObjectWithRefs>(firstLegLastFilterTag_);
0034   secondLegLastFilterToken_ = consumes<trigger::TriggerFilterObjectWithRefs>(secondLegLastFilterTag_);
0035 }
0036
0037 void HLTEgammaDoubleLegCombFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0038   edm::ParameterSetDescription desc;
0039   makeHLTFilterDescription(desc);
0040   desc.add<edm::InputTag>("firstLegLastFilter", edm::InputTag("firstFilter"));
0041   desc.add<edm::InputTag>("secondLegLastFilter", edm::InputTag("secondFilter"));
0042   desc.add<int>("nrRequiredFirstLeg", 0);
0043   desc.add<int>("nrRequiredSecondLeg", 0);
0044   desc.add<int>("nrRequiredUniqueLeg", 0);
0045   desc.add<double>("maxMatchDR", 0.01);
0046   descriptions.add("hltEgammaDoubleLegCombFilter", desc);
0047 }
0048
0049 HLTEgammaDoubleLegCombFilter::~HLTEgammaDoubleLegCombFilter() = default;
0050
0051 // ------------ method called to produce the data  ------------
0052 bool HLTEgammaDoubleLegCombFilter::hltFilter(edm::Event& iEvent,
0053                                              const edm::EventSetup& iSetup,
0054                                              trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0055   //right, issue 1, we dont know if this is a TriggerElectron, TriggerPhoton, TriggerCluster (should never be a TriggerCluster btw as that implies the 4-vectors are not stored in AOD)
0056
0057   //trigger::TriggerObjectType firstLegTrigType;
0058   std::vector<math::XYZPoint> firstLegP3s;
0059
0060   //trigger::TriggerObjectType secondLegTrigType;
0061   std::vector<math::XYZPoint> secondLegP3s;
0062
0063   getP3OfLegCands(iEvent, firstLegLastFilterToken_, firstLegP3s);
0064   getP3OfLegCands(iEvent, secondLegLastFilterToken_, secondLegP3s);
0065
0066   std::vector<std::pair<int, int> > matchedCands;
0067   matchCands(firstLegP3s, secondLegP3s, matchedCands);
0068
0069   int nr1stLegOnly = 0;
0070   int nr2ndLegOnly = 0;
0071   int nrBoth = 0;
0072   ;
0073
0074   for (auto& matchedCand : matchedCands) {
0075     if (matchedCand.first >= 0) {  //we found a first leg cand
0076       if (matchedCand.second >= 0)
0077         nrBoth++;  //we also found a second leg cand
0078       else
0079         nr1stLegOnly++;  //we didnt find a second leg cand
0080     } else if (matchedCand.second >= 0)
0081       nr2ndLegOnly++;  //we found a second leg cand but we didnt find a first leg
0082   }
0083
0084   bool accept = false;
0085   if (nr1stLegOnly + nr2ndLegOnly + nrBoth >= nrRequiredUniqueLeg_) {
0086     if (nr1stLegOnly >= nrRequiredFirstLeg_ && nr2ndLegOnly >= nrRequiredSecondLeg_)
0087       accept = true;
0088     else {
0089       int nrNeededFirstLeg = std::max(0, nrRequiredFirstLeg_ - nr1stLegOnly);
0090       int nrNeededSecondLeg = std::max(0, nrRequiredSecondLeg_ - nr2ndLegOnly);
0091
0092       if (nrBoth >= nrNeededFirstLeg + nrNeededSecondLeg)
0093         accept = true;
0094     }
0095   }
0096
0097   return accept;
0098 }
0099
0100 //-1 if no match is found
0101 void HLTEgammaDoubleLegCombFilter::matchCands(const std::vector<math::XYZPoint>& firstLegP3s,
0102                                               const std::vector<math::XYZPoint>& secondLegP3s,
0103                                               std::vector<std::pair<int, int> >& matchedCands) const {
0104   std::vector<size_t> matched2ndLegs;
0105   for (size_t firstLegNr = 0; firstLegNr < firstLegP3s.size(); firstLegNr++) {
0106     int matchedNr = -1;
0107     for (size_t secondLegNr = 0; secondLegNr < secondLegP3s.size() && matchedNr == -1; secondLegNr++) {
0108       if (reco::deltaR2(firstLegP3s[firstLegNr], secondLegP3s[secondLegNr]) < maxMatchDR_ * maxMatchDR_)
0109         matchedNr = secondLegNr;
0110     }
0111     matchedCands.push_back(std::make_pair(firstLegNr, matchedNr));
0112     if (matchedNr >= 0)
0113       matched2ndLegs.push_back(static_cast<size_t>(matchedNr));
0114   }
0115   std::sort(matched2ndLegs.begin(), matched2ndLegs.end());
0116
0117   for (size_t secondLegNr = 0; secondLegNr < secondLegP3s.size(); secondLegNr++) {
0118     if (!std::binary_search(matched2ndLegs.begin(), matched2ndLegs.end(), secondLegNr)) {  //wasnt matched already
0119       matchedCands.push_back(std::make_pair<int, int>(-1, static_cast<int>(secondLegNr)));
0120     }
0121   }
0122 }
0123
0124 //we use position and p3 interchangably here, we only use eta/phi so its alright
0125 void HLTEgammaDoubleLegCombFilter::getP3OfLegCands(
0126     const edm::Event& iEvent,
0127     const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs>& filterToken,
0128     std::vector<math::XYZPoint>& p3s) {
0129   edm::Handle<trigger::TriggerFilterObjectWithRefs> filterOutput;
0130   iEvent.getByToken(filterToken, filterOutput);
0131
0132   //its easier on the if statement flow if I try everything at once, shouldnt add to timing
0133   std::vector<edm::Ref<reco::RecoEcalCandidateCollection> > phoCands;
0134   filterOutput->getObjects(trigger::TriggerPhoton, phoCands);
0135   std::vector<edm::Ref<reco::RecoEcalCandidateCollection> > clusCands;
0136   filterOutput->getObjects(trigger::TriggerCluster, clusCands);
0137   std::vector<edm::Ref<reco::ElectronCollection> > eleCands;
0138   filterOutput->getObjects(trigger::TriggerElectron, eleCands);
0139   std::vector<edm::Ref<reco::CaloJetCollection> > jetCands;
0140   filterOutput->getObjects(trigger::TriggerJet, jetCands);
0141
0142   if (!phoCands.empty()) {  //its photons
0143     for (auto& phoCand : phoCands) {
0144       p3s.push_back(phoCand->superCluster()->position());
0145     }
0146   } else if (!clusCands.empty()) {
0147     //try trigger cluster (should never be this, at the time of writing (17/1/11) this would indicate an error)
0148     for (auto& clusCand : clusCands) {
0149       p3s.push_back(clusCand->superCluster()->position());
0150     }
0151   } else if (!eleCands.empty()) {
0152     for (auto& eleCand : eleCands) {
0153       p3s.push_back(eleCand->superCluster()->position());
0154     }
0155   } else if (!jetCands.empty()) {
0156     for (auto& jetCand : jetCands) {
0157       math::XYZPoint p3;
0158       p3.SetX(jetCand->p4().x());
0159       p3.SetY(jetCand->p4().y());
0160       p3.SetZ(jetCand->p4().z());
0161       p3s.push_back(p3);
0162     }
0163   }
0164 }
0165
0166 // declare this class as a framework plugin
0167 #include "FWCore/Framework/interface/MakerMacros.h"
0168 DEFINE_FWK_MODULE(HLTEgammaDoubleLegCombFilter);