HcalAlCaRecoProducers/plugins/AlCaIsoTracksFilter.cc

0001 // system include files
0002 #include <algorithm>
0003 #include <atomic>
0004 #include <memory>
0005 #include <cmath>
0006 #include <iostream>
0007 #include <sstream>
0008 #include <fstream>
0009
0010 // user include files
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/stream/EDFilter.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/Run.h"
0015 #include "FWCore/Framework/interface/LuminosityBlock.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0018 #include "FWCore/Common/interface/TriggerNames.h"
0019
0020 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0021 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0022
0023 #include "DataFormats/Common/interface/Handle.h"
0024 //Tracks
0025 #include "DataFormats/TrackReco/interface/Track.h"
0026 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0027 #include "DataFormats/TrackReco/interface/HitPattern.h"
0028 #include "DataFormats/TrackReco/interface/TrackBase.h"
0029 // Vertices
0030 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0031 #include "DataFormats/VertexReco/interface/Vertex.h"
0032 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0033 // RecHits
0034 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0035 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0036 //Triggers
0037 #include "DataFormats/Common/interface/TriggerResults.h"
0038 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0039
0040 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0041
0042 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0043 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0044 #include "Calibration/IsolatedParticles/interface/eCone.h"
0045 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0046
0047 #include "MagneticField/Engine/interface/MagneticField.h"
0048 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0049 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0050 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0051 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0052 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0053 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0054 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0055 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0056 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0057 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0058 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0059
0060 //#define EDM_ML_DEBUG
0061 //
0062 // class declaration
0063 //
0064
0065 namespace alCaIsoTracksFilter {
0066   struct Counters {
0067     Counters() : nAll_(0), nGood_(0), nRange_(0), nHigh_(0) {}
0068     mutable std::atomic<unsigned int> nAll_, nGood_, nRange_, nHigh_;
0069   };
0070 }  // namespace alCaIsoTracksFilter
0071
0072 class AlCaIsoTracksFilter : public edm::stream::EDFilter<edm::GlobalCache<alCaIsoTracksFilter::Counters>> {
0073 public:
0074   explicit AlCaIsoTracksFilter(edm::ParameterSet const&, const alCaIsoTracksFilter::Counters* count);
0075   ~AlCaIsoTracksFilter() override = default;
0076
0077   static std::unique_ptr<alCaIsoTracksFilter::Counters> initializeGlobalCache(edm::ParameterSet const& iConfig) {
0078     return std::make_unique<alCaIsoTracksFilter::Counters>();
0079   }
0080
0081   bool filter(edm::Event&, edm::EventSetup const&) override;
0082   void endStream() override;
0083   static void globalEndJob(const alCaIsoTracksFilter::Counters* counters);
0084   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0085
0086 private:
0087   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0088   void endRun(edm::Run const&, edm::EventSetup const&) override;
0089
0090   // ----------member data ---------------------------
0091   HLTConfigProvider hltConfig_;
0092   const std::vector<std::string> trigNames_;
0093   const edm::InputTag labelGenTrack_, labelRecVtx_;
0094   const edm::InputTag labelEB_, labelEE_, labelHBHE_;
0095   const edm::InputTag triggerEvent_, theTriggerResultsLabel_;
0096   const std::string processName_;
0097   const double a_coneR_, a_mipR_, pTrackMin_, eEcalMax_;
0098   const double maxRestrictionP_, slopeRestrictionP_;
0099   const double eIsolate_;
0100   const double hitEthrEB_, hitEthrEE0_, hitEthrEE1_;
0101   const double hitEthrEE2_, hitEthrEE3_;
0102   const double hitEthrEELo_, hitEthrEEHi_;
0103   const double pTrackLow_, pTrackHigh_, pTrackH_;
0104   const int preScale_, preScaleH_;
0105   const std::string theTrackQuality_;
0106   const std::vector<int> debEvents_;
0107   const bool usePFThresh_;
0108   spr::trackSelectionParameters selectionParameter_;
0109   double a_charIsoR_;
0110   unsigned int nRun_, nAll_, nGood_, nRange_, nHigh_;
0111   edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
0112   edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0113   edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0114   edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0115   edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0116   edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0117   edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0118   edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0119   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0120   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0121   edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> tok_ecalPFRecHitThresholds_;
0122
0123   const EcalPFRecHitThresholds* eThresholds_;
0124 };
0125
0126 //
0127 // constants, enums and typedefs
0128 //
0129
0130 //
0131 // static data member definitions
0132 //
0133
0134 //
0135 // constructors and destructor
0136 //
0137 AlCaIsoTracksFilter::AlCaIsoTracksFilter(const edm::ParameterSet& iConfig, const alCaIsoTracksFilter::Counters* count)
0138     : trigNames_(iConfig.getParameter<std::vector<std::string>>("triggers")),
0139       labelGenTrack_(iConfig.getParameter<edm::InputTag>("labelTrack")),
0140       labelRecVtx_(iConfig.getParameter<edm::InputTag>("labelVertex")),
0141       labelEB_(iConfig.getParameter<edm::InputTag>("labelEBRecHit")),
0142       labelEE_(iConfig.getParameter<edm::InputTag>("labelEERecHit")),
0143       labelHBHE_(iConfig.getParameter<edm::InputTag>("labelHBHERecHit")),
0144       triggerEvent_(iConfig.getParameter<edm::InputTag>("labelTriggerEvent")),
0145       theTriggerResultsLabel_(iConfig.getParameter<edm::InputTag>("labelTriggerResult")),
0146       processName_(iConfig.getParameter<std::string>("processName")),
0147       a_coneR_(iConfig.getParameter<double>("coneRadius")),
0148       a_mipR_(iConfig.getParameter<double>("coneRadiusMIP")),
0149       pTrackMin_(iConfig.getParameter<double>("minimumTrackP")),
0150       eEcalMax_(iConfig.getParameter<double>("maximumEcalEnergy")),
0151       maxRestrictionP_(iConfig.getParameter<double>("maxTrackP")),
0152       slopeRestrictionP_(iConfig.getParameter<double>("slopeTrackP")),
0153       eIsolate_(iConfig.getParameter<double>("isolationEnergy")),
0154       hitEthrEB_(iConfig.getParameter<double>("EBHitEnergyThreshold")),
0155       hitEthrEE0_(iConfig.getParameter<double>("EEHitEnergyThreshold0")),
0156       hitEthrEE1_(iConfig.getParameter<double>("EEHitEnergyThreshold1")),
0157       hitEthrEE2_(iConfig.getParameter<double>("EEHitEnergyThreshold2")),
0158       hitEthrEE3_(iConfig.getParameter<double>("EEHitEnergyThreshold3")),
0159       hitEthrEELo_(iConfig.getParameter<double>("EEHitEnergyThresholdLow")),
0160       hitEthrEEHi_(iConfig.getParameter<double>("EEHitEnergyThresholdHigh")),
0161       pTrackLow_(iConfig.getParameter<double>("momentumRangeLow")),
0162       pTrackHigh_(iConfig.getParameter<double>("momentumRangeHigh")),
0163       pTrackH_(iConfig.getParameter<double>("momentumHigh")),
0164       preScale_(iConfig.getParameter<int>("preScaleFactor")),
0165       preScaleH_(iConfig.getParameter<int>("preScaleHigh")),
0166       theTrackQuality_(iConfig.getParameter<std::string>("trackQuality")),
0167       debEvents_(iConfig.getParameter<std::vector<int>>("debugEvents")),
0168       usePFThresh_(iConfig.getParameter<bool>("usePFThreshold")),
0169       nRun_(0),
0170       nAll_(0),
0171       nGood_(0),
0172       nRange_(0),
0173       nHigh_(0) {
0174   //now do what ever initialization is needed
0175   const double isolationRadius(28.9);
0176   // Different isolation cuts are described in DN-2016/029
0177   // Tight cut uses 2 GeV; Loose cut uses 10 GeV
0178   // Eta dependent cut uses (maxRestrictionP_ * exp(|ieta|*log(2.5)/18))
0179   // with the factor for exponential slopeRestrictionP_ = log(2.5)/18
0180   // maxRestrictionP_ = 8 GeV as came from a study
0181   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality_);
0182   selectionParameter_.minPt = iConfig.getParameter<double>("minTrackPt");
0183   ;
0184   selectionParameter_.minQuality = trackQuality_;
0185   selectionParameter_.maxDxyPV = iConfig.getParameter<double>("maxDxyPV");
0186   selectionParameter_.maxDzPV = iConfig.getParameter<double>("maxDzPV");
0187   selectionParameter_.maxChi2 = iConfig.getParameter<double>("maxChi2");
0188   selectionParameter_.maxDpOverP = iConfig.getParameter<double>("maxDpOverP");
0189   selectionParameter_.minOuterHit = iConfig.getParameter<int>("minOuterHit");
0190   selectionParameter_.minLayerCrossed = iConfig.getParameter<int>("minLayerCrossed");
0191   selectionParameter_.maxInMiss = iConfig.getParameter<int>("maxInMiss");
0192   selectionParameter_.maxOutMiss = iConfig.getParameter<int>("maxOutMiss");
0193   a_charIsoR_ = a_coneR_ + isolationRadius;
0194
0195   // define tokens for access
0196   tok_trigEvt_ = consumes<trigger::TriggerEvent>(triggerEvent_);
0197   tok_trigRes_ = consumes<edm::TriggerResults>(theTriggerResultsLabel_);
0198   tok_genTrack_ = consumes<reco::TrackCollection>(labelGenTrack_);
0199   tok_recVtx_ = consumes<reco::VertexCollection>(labelRecVtx_);
0200   tok_bs_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("labelBeamSpot"));
0201
0202   tok_EB_ = consumes<EcalRecHitCollection>(labelEB_);
0203   tok_EE_ = consumes<EcalRecHitCollection>(labelEE_);
0204   tok_hbhe_ = consumes<HBHERecHitCollection>(labelHBHE_);
0205
0206   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0207   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0208   tok_ecalPFRecHitThresholds_ = esConsumes<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd>();
0209
0210   edm::LogVerbatim("HcalIsoTrack") << "Parameters read from config file \n"
0211                                    << "\t minPt " << selectionParameter_.minPt << "\t theTrackQuality "
0212                                    << theTrackQuality_ << "\t minQuality " << selectionParameter_.minQuality
0213                                    << "\t maxDxyPV " << selectionParameter_.maxDxyPV << "\t maxDzPV "
0214                                    << selectionParameter_.maxDzPV << "\t maxChi2 " << selectionParameter_.maxChi2
0215                                    << "\t maxDpOverP " << selectionParameter_.maxDpOverP << "\t minOuterHit "
0216                                    << selectionParameter_.minOuterHit << "\t minLayerCrossed "
0217                                    << selectionParameter_.minLayerCrossed << "\t maxInMiss "
0218                                    << selectionParameter_.maxInMiss << "\t maxOutMiss "
0219                                    << selectionParameter_.maxOutMiss << "\n"
0220                                    << "\t a_coneR " << a_coneR_ << "\t a_charIsoR " << a_charIsoR_ << "\t a_mipR "
0221                                    << a_mipR_ << "\t maxRestrictionP_ " << maxRestrictionP_ << "\t slopeRestrictionP_ "
0222                                    << slopeRestrictionP_ << "\t eIsolate_ " << eIsolate_ << "\n"
0223                                    << "\t Precale factor " << preScale_ << "\t in momentum range " << pTrackLow_ << ":"
0224                                    << pTrackHigh_ << " and prescale factor " << preScaleH_ << " for p > " << pTrackH_
0225                                    << " Threshold flag used " << usePFThresh_ << " value for EB " << hitEthrEB_
0226                                    << " EE " << hitEthrEE0_ << ":" << hitEthrEE1_ << ":" << hitEthrEE2_ << ":"
0227                                    << hitEthrEE3_ << ":" << hitEthrEELo_ << ":" << hitEthrEEHi_ << " and "
0228                                    << debEvents_.size() << " events to be debugged";
0229
0230   for (unsigned int k = 0; k < trigNames_.size(); ++k)
0231     edm::LogVerbatim("HcalIsoTrack") << "Trigger[" << k << "] " << trigNames_[k];
0232 }  // AlCaIsoTracksFilter::AlCaIsoTracksFilter  constructor
0233
0234 //
0235 // member functions
0236 //
0237
0238 // ------------ method called on each new Event  ------------
0239 bool AlCaIsoTracksFilter::filter(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0240   bool accept(false);
0241   ++nAll_;
0242 #ifdef EDM_ML_DEBUG
0243   bool debug = (debEvents_.empty())
0244                    ? true
0245                    : (std::find(debEvents_.begin(), debEvents_.end(), iEvent.id().event()) != debEvents_.end());
0246   if (debug)
0247     edm::LogVerbatim("HcalIsoTrack") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event()
0248                                      << " Luminosity " << iEvent.luminosityBlock() << " Bunch "
0249                                      << iEvent.bunchCrossing();
0250 #endif
0251
0252   // get Ecal Thresholds
0253   eThresholds_ = &iSetup.getData(tok_ecalPFRecHitThresholds_);
0254
0255   //Step1: Find if the event passes one of the chosen triggers
0256   bool triggerSatisfied(false);
0257   if (trigNames_.empty()) {
0258     triggerSatisfied = true;
0259   } else {
0260     trigger::TriggerEvent triggerEvent;
0261     auto const& triggerEventHandle = iEvent.getHandle(tok_trigEvt_);
0262     if (!triggerEventHandle.isValid()) {
0263       edm::LogWarning("HcalIsoTrack") << "Error! Can't get the product " << triggerEvent_.label();
0264     } else {
0265       triggerEvent = *(triggerEventHandle.product());
0266
0267       /////////////////////////////TriggerResults
0268       auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
0269       if (triggerResults.isValid()) {
0270         std::vector<std::string> modules;
0271         const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0272         const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
0273         for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0274           int hlt = triggerResults->accept(iHLT);
0275           for (unsigned int i = 0; i < trigNames_.size(); ++i) {
0276             if (triggerNames_[iHLT].find(trigNames_[i]) != std::string::npos) {
0277               if (hlt > 0)
0278                 triggerSatisfied = true;
0279 #ifdef EDM_ML_DEBUG
0280               if (debug)
0281                 edm::LogVerbatim("HcalIsoTrack")
0282                     << triggerNames_[iHLT] << " has got HLT flag " << hlt << ":" << triggerSatisfied;
0283 #endif
0284               if (triggerSatisfied)
0285                 break;
0286             }
0287           }
0288         }
0289       }
0290     }
0291   }
0292 #ifdef EDM_ML_DEBUG
0293   if (debug)
0294     edm::LogVerbatim("HcalIsoTrack") << "AlCaIsoTracksFilter:: triggerSatisfied: " << triggerSatisfied;
0295 #endif
0296
0297   //Step2: Get geometry/B-field information
0298   if (triggerSatisfied) {
0299     //Get magnetic field
0300     const MagneticField* bField = &(iSetup.getData(tok_magField_));
0301     const CaloGeometry* geo = &(iSetup.getData(tok_geom_));
0302
0303     //Also relevant information to extrapolate tracks to Hcal surface
0304     bool foundCollections(true);
0305     //Get track collection
0306     auto trkCollection = iEvent.getHandle(tok_genTrack_);
0307     if (!trkCollection.isValid()) {
0308       edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
0309       foundCollections = false;
0310     }
0311
0312     //Define the best vertex and the beamspot
0313     auto const& recVtxs = iEvent.getHandle(tok_recVtx_);
0314     auto const& beamSpotH = iEvent.getHandle(tok_bs_);
0315     math::XYZPoint leadPV(0, 0, 0);
0316     if (!recVtxs->empty() && !((*recVtxs)[0].isFake())) {
0317       leadPV = math::XYZPoint((*recVtxs)[0].x(), (*recVtxs)[0].y(), (*recVtxs)[0].z());
0318     } else if (beamSpotH.isValid()) {
0319       leadPV = beamSpotH->position();
0320     }
0321 #ifdef EDM_ML_DEBUG
0322     if (debug)
0323       edm::LogVerbatim("HcalIsoTrack") << "Primary Vertex " << leadPV;
0324 #endif
0325
0326     // RecHits
0327     auto barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0328     if (!barrelRecHitsHandle.isValid()) {
0329       edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEB_;
0330       foundCollections = false;
0331     }
0332     auto endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0333     if (!endcapRecHitsHandle.isValid()) {
0334       edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEE_;
0335       foundCollections = false;
0336     }
0337     auto hbhe = iEvent.getHandle(tok_hbhe_);
0338     if (!hbhe.isValid()) {
0339       edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelHBHE_;
0340       foundCollections = false;
0341     }
0342 #ifdef EDM_ML_DEBUG
0343     if (debug)
0344       edm::LogVerbatim("HcalIsoTrack") << "AlCaIsoTracksFilter:: foundCollections: " << foundCollections;
0345 #endif
0346
0347     //Step3 propagate the tracks to calorimeter surface and find
0348     // candidates for isolated tracks
0349     if (foundCollections) {
0350       //Propagate tracks to calorimeter surface)
0351       std::vector<spr::propagatedTrackDirection> trkCaloDirections;
0352       spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, false);
0353
0354       std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
0355 #ifdef EDM_ML_DEBUG
0356       if (debug)
0357         edm::LogVerbatim("HcalIsoTrack") << "AlCaIsoTracksFilter:: Has " << trkCaloDirections.size()
0358                                          << " propagated tracks from a total of " << trkCollection->size();
0359 #endif
0360       unsigned int nTracks(0), ntrin(0), ntrout(0), ntrH(0);
0361       for (trkDetItr = trkCaloDirections.begin(), nTracks = 0; trkDetItr != trkCaloDirections.end();
0362            trkDetItr++, nTracks++) {
0363         const reco::Track* pTrack = &(*(trkDetItr->trkItr));
0364         math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
0365 #ifdef EDM_ML_DEBUG
0366         if (debug)
0367           edm::LogVerbatim("HcalIsoTrack")
0368               << "This track : " << nTracks << " (pt|eta|phi|p) : " << pTrack->pt() << "|" << pTrack->eta() << "|"
0369               << pTrack->phi() << "|" << pTrack->p() << " OK HCAL " << trkDetItr->okHCAL;
0370 #endif
0371         //Selection of good track
0372         int ieta(0);
0373         if (trkDetItr->okHCAL) {
0374           HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
0375           ieta = detId.ietaAbs();
0376         }
0377         bool qltyFlag = spr::goodTrack(pTrack, leadPV, selectionParameter_, false);
0378 #ifdef EDM_ML_DEBUG
0379         if (debug)
0380           edm::LogVerbatim("HcalIsoTrack")
0381               << "qltyFlag|okECAL|okHCAL : " << qltyFlag << "|" << trkDetItr->okECAL << "|" << trkDetItr->okHCAL;
0382 #endif
0383         if (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL) {
0384           double t_p = pTrack->p();
0385           int nNearTRKs(0);
0386           std::vector<DetId> eIds;
0387           std::vector<double> eHit;
0388 #ifdef EDM_ML_DEBUG
0389           double eEcal =
0390 #endif
0391               spr::eCone_ecal(geo,
0392                               barrelRecHitsHandle,
0393                               endcapRecHitsHandle,
0394                               trkDetItr->pointHCAL,
0395                               trkDetItr->pointECAL,
0396                               a_mipR_,
0397                               trkDetItr->directionECAL,
0398                               eIds,
0399                               eHit);
0400           double eMipDR(0);
0401           for (unsigned int k = 0; k < eIds.size(); ++k) {
0402             double eThr(hitEthrEB_);
0403             if (usePFThresh_) {
0404               eThr = static_cast<double>((*eThresholds_)[eIds[k]]);
0405             } else {
0406               const GlobalPoint& pos = geo->getPosition(eIds[k]);
0407               double eta = std::abs(pos.eta());
0408               if (eIds[k].subdetId() != EcalBarrel) {
0409                 eThr = (((eta * hitEthrEE3_ + hitEthrEE2_) * eta + hitEthrEE1_) * eta + hitEthrEE0_);
0410                 if (eThr < hitEthrEELo_)
0411                   eThr = hitEthrEELo_;
0412                 else if (eThr > hitEthrEEHi_)
0413                   eThr = hitEthrEEHi_;
0414               }
0415             }
0416             if (eHit[k] > eThr)
0417               eMipDR += eHit[k];
0418           }
0419           double hmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, false);
0420           double eIsolation = (maxRestrictionP_ * exp(slopeRestrictionP_ * ((double)(ieta))));
0421           if (eIsolation < eIsolate_)
0422             eIsolation = eIsolate_;
0423 #ifdef EDM_ML_DEBUG
0424           std::string ctype =
0425               (t_p > pTrackMin_ && eMipDR < eEcalMax_ && hmaxNearP < eIsolation) ? " ***** ACCEPT *****" : "";
0426           if (debug)
0427             edm::LogVerbatim("HcalIsoTrack")
0428                 << "This track : " << nTracks << " (pt|eta|phi|p) : " << pTrack->pt() << "|" << pTrack->eta() << "|"
0429                 << pTrack->phi() << "|" << t_p << " e_MIP " << eMipDR << ":" << eEcal << " Chg Isolation " << hmaxNearP
0430                 << ":" << eIsolation << ctype;
0431 #endif
0432           if (t_p > pTrackMin_ && eMipDR < eEcalMax_ && hmaxNearP < eIsolation) {
0433             if (t_p > pTrackLow_ && t_p < pTrackHigh_)
0434               ntrin++;
0435             else if (t_p > pTrackH_)
0436               ntrH++;
0437             else
0438               ntrout++;
0439           }
0440         }
0441       }
0442       accept = (ntrout > 0);
0443       if (!accept && ntrin > 0) {
0444         ++nRange_;
0445         if (preScale_ <= 1)
0446           accept = true;
0447         else if (nRange_ % preScale_ == 1)
0448           accept = true;
0449       }
0450       if (!accept && ntrH > 0) {
0451         ++nHigh_;
0452         if (preScaleH_ <= 1)
0453           accept = true;
0454         else if (nHigh_ % preScaleH_ == 1)
0455           accept = true;
0456       }
0457 #ifdef EDM_ML_DEBUG
0458       if (debug)
0459         edm::LogVerbatim("HcalIsoTrack") << "Summary Range " << ntrout << " Low " << ntrin << " High " << ntrH
0460                                          << " Accept " << accept;
0461 #endif
0462     }
0463   }
0464   // Step 4:  Return the acceptance flag
0465   if (accept) {
0466     ++nGood_;
0467     edm::LogVerbatim("HcalIsoTrackX") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event();
0468   }
0469   return accept;
0470
0471 }  // AlCaIsoTracksFilter::filter
0472
0473 // ------------ method called once each job just after ending the event loop  ------------
0474 void AlCaIsoTracksFilter::endStream() {
0475   globalCache()->nAll_ += nAll_;
0476   globalCache()->nGood_ += nGood_;
0477   globalCache()->nRange_ += nRange_;
0478   globalCache()->nHigh_ += nHigh_;
0479 }
0480
0481 void AlCaIsoTracksFilter::globalEndJob(const alCaIsoTracksFilter::Counters* count) {
0482   edm::LogVerbatim("HcalIsoTrack") << "Selects " << count->nGood_ << " in " << count->nAll_ << " events and with "
0483                                    << count->nRange_ << " events in the p-range" << count->nHigh_
0484                                    << " events with high p";
0485 }
0486
0487 // ------------ method called when starting to processes a run  ------------
0488 void AlCaIsoTracksFilter::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0489   bool changed(false);
0490   edm::LogVerbatim("HcalIsoTrack") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init "
0491                                    << hltConfig_.init(iRun, iSetup, processName_, changed);
0492 }
0493
0494 // ------------ method called when ending the processing of a run  ------------
0495 void AlCaIsoTracksFilter::endRun(edm::Run const& iRun, edm::EventSetup const&) {
0496   ++nRun_;
0497   edm::LogVerbatim("HcalIsoTrack") << "endRun[" << nRun_ << "] " << iRun.run();
0498 }
0499
0500 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0501 void AlCaIsoTracksFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0502   edm::ParameterSetDescription desc;
0503   desc.add<edm::InputTag>("labelTrack", edm::InputTag("generalTracks"));
0504   desc.add<edm::InputTag>("labelVertex", edm::InputTag("offlinePrimaryVertices"));
0505   desc.add<edm::InputTag>("labelBeamSpot", edm::InputTag("offlineBeamSpot"));
0506   desc.add<edm::InputTag>("labelEBRecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0507   desc.add<edm::InputTag>("labelEERecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0508   desc.add<edm::InputTag>("labelHBHERecHit", edm::InputTag("hbhereco"));
0509   desc.add<edm::InputTag>("labelTriggerEvent", edm::InputTag("hltTriggerSummaryAOD", "", "HLT"));
0510   desc.add<edm::InputTag>("labelTriggerResult", edm::InputTag("TriggerResults", "", "HLT"));
0511   std::vector<std::string> trigger;
0512   desc.add<std::vector<std::string>>("triggers", trigger);
0513   desc.add<std::string>("processName", "HLT");
0514   // following 10 parameters are parameters to select good tracks
0515   desc.add<std::string>("trackQuality", "highPurity");
0516   desc.add<double>("minTrackPt", 1.0);
0517   desc.add<double>("maxDxyPV", 10.0);
0518   desc.add<double>("maxDzPV", 100.0);
0519   desc.add<double>("maxChi2", 5.0);
0520   desc.add<double>("maxDpOverP", 0.1);
0521   desc.add<int>("minOuterHit", 4);
0522   desc.add<int>("minLayerCrossed", 8);
0523   desc.add<int>("maxInMiss", 2);
0524   desc.add<int>("maxOutMiss", 2);
0525   // Minimum momentum of selected isolated track and signal zone
0526   desc.add<double>("coneRadius", 34.98);
0527   desc.add<double>("minimumTrackP", 20.0);
0528   // signal zone in ECAL and MIP energy cutoff
0529   desc.add<double>("coneRadiusMIP", 14.0);
0530   desc.add<double>("maximumEcalEnergy", 100.0);
0531   // following 3 parameters are for isolation cuts and described in the code
0532   desc.add<double>("maxTrackP", 8.0);
0533   desc.add<double>("slopeTrackP", 0.05090504066);
0534   desc.add<double>("isolationEnergy", 10.0);
0535   // energy thershold for ECAL (from Egamma group)
0536   desc.add<double>("EBHitEnergyThreshold", 0.08);
0537   desc.add<double>("EEHitEnergyThreshold0", 0.30);
0538   desc.add<double>("EEHitEnergyThreshold1", 0.00);
0539   desc.add<double>("EEHitEnergyThreshold2", 0.00);
0540   desc.add<double>("EEHitEnergyThreshold3", 0.00);
0541   desc.add<double>("EEHitEnergyThresholdLow", 0.30);
0542   desc.add<double>("EEHitEnergyThresholdHigh", 0.30);
0543   // Prescale events only containing isolated tracks in the range
0544   desc.add<double>("momentumRangeLow", 20.0);
0545   desc.add<double>("momentumRangeHigh", 40.0);
0546   desc.add<int>("preScaleFactor", 10);
0547   desc.add<double>("momentumHigh", 60.0);
0548   desc.add<int>("preScaleHigh", 5);
0549   std::vector<int> events;
0550   desc.add<std::vector<int>>("debugEvents", events);
0551   desc.add<bool>("usePFThreshold", true);
0552   descriptions.add("alcaIsoTracksFilter", desc);
0553 }
0554
0555 //define this as a plug-in
0556 #include "FWCore/Framework/interface/MakerMacros.h"
0557 DEFINE_FWK_MODULE(AlCaIsoTracksFilter);