DeDx/plugins/HLTDeDxFilter.cc

0001 /** \class HLTDeDxFilter
0002 *
0003 *
0004 *  \author Claude Nuttens
0005 *
0006 */
0007
0008 #include "RecoTracker/DeDx/plugins/HLTDeDxFilter.h"
0009
0010 #include "DataFormats/Common/interface/Handle.h"
0011
0012 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 //#include "DataFormats/JetReco/interface/CaloJetCollection.h"
0017 #include "DataFormats/TrackReco/interface/Track.h"
0018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0019
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "DataFormats/TrackReco/interface/DeDxData.h"
0023 //#include "DataFormats/Math/interface/deltaPhi.h"
0024 #include "DataFormats/Math/interface/deltaR.h"
0025 #include "DataFormats/Common/interface/ValueMap.h"
0026 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0028 #include "FWCore/Utilities/interface/InputTag.h"
0029 #include <vector>
0030 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0031 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0032
0033 //
0034 // constructors and destructor
0035 //
0036 HLTDeDxFilter::HLTDeDxFilter(const edm::ParameterSet& iConfig) : HLTFilter(iConfig) {
0037   minDEDx_ = iConfig.getParameter<double>("minDEDx");
0038   minPT_ = iConfig.getParameter<double>("minPT");
0039   minNOM_ = iConfig.getParameter<double>("minNOM");
0040   maxETA_ = iConfig.getParameter<double>("maxETA");
0041   minNumValidHits_ = iConfig.getParameter<double>("minNumValidHits");
0042   maxNHitMissIn_ = iConfig.getParameter<double>("maxNHitMissIn");
0043   maxNHitMissMid_ = iConfig.getParameter<double>("maxNHitMissMid");
0044   maxRelTrkIsoDeltaRp3_ = iConfig.getParameter<double>("maxRelTrkIsoDeltaRp3");
0045   relTrkIsoDeltaRSize_ = iConfig.getParameter<double>("relTrkIsoDeltaRSize");
0046   maxAssocCaloE_ = iConfig.getParameter<double>("maxAssocCaloE");
0047   maxAssocCaloEDeltaRSize_ = iConfig.getParameter<double>("maxAssocCaloEDeltaRSize");
0048   inputTracksTag_ = iConfig.getParameter<edm::InputTag>("inputTracksTag");
0049   inputdedxTag_ = iConfig.getParameter<edm::InputTag>("inputDeDxTag");
0050   caloTowersTag_ = iConfig.getParameter<edm::InputTag>("caloTowersTag");
0051
0052   if (maxAssocCaloE_ >= 0)
0053     caloTowersToken_ = consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("caloTowersTag"));
0054   inputTracksToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("inputTracksTag"));
0055   inputdedxToken_ = consumes<edm::ValueMap<reco::DeDxData>>(iConfig.getParameter<edm::InputTag>("inputDeDxTag"));
0056
0057   thisModuleTag_ = edm::InputTag(iConfig.getParameter<std::string>("@module_label"));
0058
0059   //register your products
0060   produces<reco::RecoChargedCandidateCollection>();
0061 }
0062
0063 HLTDeDxFilter::~HLTDeDxFilter() {}
0064
0065 void HLTDeDxFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0066   edm::ParameterSetDescription desc;
0067   desc.add<bool>("saveTags", false);
0068   desc.add<double>("minDEDx", 0.0);
0069   desc.add<double>("minPT", 0.0);
0070   desc.add<double>("minNOM", 0.0);
0071   desc.add<double>("maxETA", 5.5);
0072   desc.add<double>("minNumValidHits", 0);
0073   desc.add<double>("maxNHitMissIn", 99);
0074   desc.add<double>("maxNHitMissMid", 99);
0075   desc.add<double>("maxRelTrkIsoDeltaRp3", -1);
0076   desc.add<double>("relTrkIsoDeltaRSize", 0.3);
0077   desc.add<double>("maxAssocCaloE", -99);
0078   desc.add<double>("maxAssocCaloEDeltaRSize", 0.5);
0079   desc.add<edm::InputTag>("caloTowersTag", edm::InputTag("hltTowerMakerForAll"));
0080   desc.add<edm::InputTag>("inputTracksTag", edm::InputTag("hltL3Mouns"));
0081   desc.add<edm::InputTag>("inputDeDxTag", edm::InputTag("HLTdedxHarm2"));
0082   descriptions.add("hltDeDxFilter", desc);
0083 }
0084
0085 // ------------ method called to produce the data  ------------
0086 bool HLTDeDxFilter::hltFilter(edm::Event& iEvent,
0087                               const edm::EventSetup& iSetup,
0088                               trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0089   using namespace std;
0090   using namespace edm;
0091   using namespace reco;
0092   using namespace trigger;
0093
0094   auto chargedCandidates = std::make_unique<std::vector<RecoChargedCandidate>>();
0095
0096   ModuleDescription moduleDesc_;
0097
0098   if (saveTags()) {
0099     filterproduct.addCollectionTag(thisModuleTag_);
0100     filterproduct.addCollectionTag(inputTracksTag_);
0101     filterproduct.addCollectionTag(inputdedxTag_);
0102   }
0103
0104   edm::Handle<reco::TrackCollection> trackCollectionHandle;
0105   iEvent.getByToken(inputTracksToken_, trackCollectionHandle);
0106   const reco::TrackCollection& trackCollection = *trackCollectionHandle.product();
0107
0108   edm::Handle<edm::ValueMap<reco::DeDxData>> dEdxTrackHandle;
0109   iEvent.getByToken(inputdedxToken_, dEdxTrackHandle);
0110   const edm::ValueMap<reco::DeDxData>& dEdxTrack = *dEdxTrackHandle.product();
0111
0112   edm::Handle<CaloTowerCollection> caloTowersHandle;
0113   if (maxAssocCaloE_ >= 0)
0114     iEvent.getByToken(caloTowersToken_, caloTowersHandle);
0115
0116   bool accept = false;
0117
0118   // early return
0119   if (trackCollection.empty())
0120     return accept;
0121
0122   //fill local arrays for eta, phi, and pt
0123   float eta[trackCollection.size()], phi[trackCollection.size()], pt[trackCollection.size()];
0124   for (unsigned int i = 0; i < trackCollection.size(); i++) {
0125     eta[i] = trackCollection[i].eta();
0126     phi[i] = trackCollection[i].phi();
0127     pt[i] = trackCollection[i].pt();
0128   }
0129   for (unsigned int i = 0; i < trackCollection.size(); i++) {
0130     reco::TrackRef track = reco::TrackRef(trackCollectionHandle, i);
0131     if (pt[i] > minPT_ && fabs(eta[i]) < maxETA_ && dEdxTrack[track].numberOfMeasurements() > minNOM_ &&
0132         dEdxTrack[track].dEdx() > minDEDx_) {
0133       if (track->numberOfValidHits() < minNumValidHits_)
0134         continue;
0135       if (track->hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS) > maxNHitMissIn_)
0136         continue;
0137       if (track->hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS) > maxNHitMissMid_)
0138         continue;
0139       if (saveTags()) {
0140         Particle::Charge q = track->charge();
0141         //SAVE DEDX INFORMATION AS IF IT WAS THE MASS OF THE PARTICLE
0142         Particle::LorentzVector p4(
0143             track->px(), track->py(), track->pz(), sqrt(pow(track->p(), 2) + pow(dEdxTrack[track].dEdx(), 2)));
0144         Particle::Point vtx(track->vx(), track->vy(), track->vz());
0145         //SAVE NOH, NOM, NOS INFORMATION AS IF IT WAS THE PDGID OF THE PARTICLE
0146         int Hits = ((dEdxTrack[track].numberOfSaturatedMeasurements() & 0xFF) << 16) |
0147                    ((dEdxTrack[track].numberOfMeasurements() & 0xFF) << 8) | (track->found() & 0xFF);
0148         RecoChargedCandidate cand(q, p4, vtx, Hits, 0);
0149         cand.setTrack(track);
0150         chargedCandidates->push_back(cand);
0151       }
0152
0153       //calculate relative trk isolation only if parameter maxRelTrkIsoDeltaRp3 is greater than 0
0154       if (maxRelTrkIsoDeltaRp3_ >= 0) {
0155         auto ptCone = trackCollection[i].pt();
0156         for (unsigned int j = 0; j < trackCollection.size(); j++) {
0157           if (i == j)
0158             continue;  // do not compare track to itself
0159           auto trkDeltaR2 = deltaR2(eta[i], phi[i], eta[j], phi[j]);
0160           if (trkDeltaR2 < relTrkIsoDeltaRSize_ * relTrkIsoDeltaRSize_) {
0161             ptCone += pt[j];
0162           }
0163         }
0164         double relTrkIso = (ptCone - pt[i]) / (pt[i]);
0165         if (relTrkIso > maxRelTrkIsoDeltaRp3_)
0166           continue;
0167       }
0168
0169       //calculate the calorimeter energy associated with the track if maxAssocCaloE_ >= 0
0170       if (maxAssocCaloE_ >= 0) {
0171         //Access info about Calo Towers
0172         double caloEMDeltaRp5 = 0;
0173         double caloHadDeltaRp5 = 0;
0174         const CaloTowerCollection& caloTower = *caloTowersHandle.product();
0175         for (CaloTowerCollection::const_iterator j = caloTower.begin(); j != caloTower.end(); j++) {
0176           auto caloDeltaR2 = deltaR2(eta[i], phi[i], j->eta(), j->phi());
0177           double Eem = j->emEnergy();
0178           double Ehad = j->hadEnergy();
0179
0180           if (caloDeltaR2 < (maxAssocCaloEDeltaRSize_ * maxAssocCaloEDeltaRSize_)) {
0181             caloEMDeltaRp5 += Eem;
0182             caloHadDeltaRp5 += Ehad;
0183           }
0184         }
0185         if (caloEMDeltaRp5 + caloHadDeltaRp5 > maxAssocCaloE_)
0186           continue;
0187       }
0188
0189       accept = true;
0190     }
0191   }
0192
0193   // put filter object into the Event
0194   if (saveTags()) {
0195     auto const chargedCandidatesHandle = iEvent.put(std::move(chargedCandidates));
0196     for (unsigned int i = 0; i < chargedCandidatesHandle->size(); ++i) {
0197       filterproduct.addObject(trigger::TriggerMuon, reco::RecoChargedCandidateRef(chargedCandidatesHandle, i));
0198     }
0199   }
0200
0201   return accept;
0202 }
0203
0204 //define this as a plug-in
0205 #include "FWCore/Framework/interface/MakerMacros.h"
0206 DEFINE_FWK_MODULE(HLTDeDxFilter);