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