HcalAlCaRecoProducers/plugins/AlCaIsoTracksProducer.cc

0001 // -*- C++ -*-
0002 //#define EDM_ML_DEBUG
0003
0004 // system include files
0005 #include <atomic>
0006 #include <memory>
0007 #include <string>
0008 #include <cmath>
0009 #include <iostream>
0010 #include <sstream>
0011 #include <fstream>
0012 #include <vector>
0013 #include <boost/regex.hpp>
0014
0015 // user include files
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/stream/EDProducer.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023 #include "FWCore/Common/interface/TriggerNames.h"
0024 #include "FWCore/Utilities/interface/Exception.h"
0025
0026 #include "DataFormats/Common/interface/Handle.h"
0027 #include "DataFormats/Common/interface/Ref.h"
0028 #include "DataFormats/DetId/interface/DetId.h"
0029 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0030 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0031 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0032 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0033 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0034 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0035 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0036 #include "DataFormats/HcalIsolatedTrack/interface/HcalIsolatedTrackCandidate.h"
0037 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0038 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0039 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
0040 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
0041 #include "DataFormats/Provenance/interface/ProductID.h"
0042 #include "DataFormats/TrackReco/interface/Track.h"
0043 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0044 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0045 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0046 #include "DataFormats/TrackReco/interface/HitPattern.h"
0047 #include "DataFormats/TrackReco/interface/TrackBase.h"
0048 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0049 #include "DataFormats/VertexReco/interface/Vertex.h"
0050 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0051 #include "DataFormats/Common/interface/TriggerResults.h"
0052
0053 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0054
0055 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0056 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0057 #include "Calibration/IsolatedParticles/interface/eCone.h"
0058 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0059
0060 #include "MagneticField/Engine/interface/MagneticField.h"
0061 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0062 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0063 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0064 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0065 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0066 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0067 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0068 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0069 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0070 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0071 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0072 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0073
0074 //#define EDM_ML_DEBUG
0075 //
0076 // class declaration
0077 //
0078
0079 namespace alCaIsoTracksProducer {
0080   struct Counters {
0081     Counters() : nAll_(0), nGood_(0), nRange_(0) {}
0082     mutable std::atomic<unsigned int> nAll_, nGood_, nRange_;
0083   };
0084 }  // namespace alCaIsoTracksProducer
0085
0086 class AlCaIsoTracksProducer : public edm::stream::EDProducer<edm::GlobalCache<alCaIsoTracksProducer::Counters> > {
0087 public:
0088   explicit AlCaIsoTracksProducer(edm::ParameterSet const&, const alCaIsoTracksProducer::Counters* count);
0089   ~AlCaIsoTracksProducer() override = default;
0090
0091   static std::unique_ptr<alCaIsoTracksProducer::Counters> initializeGlobalCache(edm::ParameterSet const&) {
0092     return std::make_unique<alCaIsoTracksProducer::Counters>();
0093   }
0094
0095   void produce(edm::Event&, edm::EventSetup const&) override;
0096   void endStream() override;
0097   static void globalEndJob(const alCaIsoTracksProducer::Counters* counters);
0098   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0099
0100 private:
0101   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0102   void endRun(edm::Run const&, edm::EventSetup const&) override;
0103   reco::HcalIsolatedTrackCandidateCollection* select(edm::Handle<edm::TriggerResults> const& triggerResults,
0104                                                      const std::vector<std::string>& triggerNames_,
0105                                                      edm::Handle<reco::TrackCollection>& trkCollection,
0106                                                      math::XYZPoint& leadPV,
0107                                                      edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
0108                                                      edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
0109                                                      edm::Handle<HBHERecHitCollection>& hbhe,
0110                                                      double ptL1,
0111                                                      double etaL1,
0112                                                      double phiL1);
0113   void setPtEtaPhi(std::vector<edm::Ref<l1extra::L1JetParticleCollection> >& objref,
0114                    double& ptL1,
0115                    double& etaL1,
0116                    double& phiL1);
0117
0118   // ----------member data ---------------------------
0119   HLTConfigProvider hltConfig_;
0120   unsigned int nRun_, nAll_, nGood_, nRange_;
0121   spr::trackSelectionParameters selectionParameter_;
0122   const std::vector<std::string> trigNames_;
0123   const std::string theTrackQuality_, processName_;
0124   const double a_coneR_, a_mipR_;
0125   const double maxRestrictionP_, slopeRestrictionP_;
0126   const double pTrackMin_, eEcalMax_, eIsolate_;
0127   const double pTrackLow_, pTrackHigh_;
0128   const int preScale_;
0129   const edm::InputTag labelGenTrack_, labelRecVtx_, labelBS_;
0130   const edm::InputTag labelEB_, labelEE_, labelHBHE_, labelHltGT_;
0131   const edm::InputTag labelTriggerEvent_, labelTriggerResults_;
0132   const std::string labelIsoTk_;
0133   double a_charIsoR_;
0134   const MagneticField* bField;
0135   const CaloGeometry* geo;
0136
0137   edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_hltGT_;
0138   edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
0139   edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0140   edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0141   edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0142   edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0143   edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0144   edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0145   edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0146   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0147   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0148 };
0149
0150 AlCaIsoTracksProducer::AlCaIsoTracksProducer(edm::ParameterSet const& iConfig,
0151                                              const alCaIsoTracksProducer::Counters* counters)
0152     : nRun_(0),
0153       nAll_(0),
0154       nGood_(0),
0155       nRange_(0),
0156       trigNames_(iConfig.getParameter<std::vector<std::string> >("triggers")),
0157       theTrackQuality_(iConfig.getParameter<std::string>("trackQuality")),
0158       processName_(iConfig.getParameter<std::string>("processName")),
0159       a_coneR_(iConfig.getParameter<double>("coneRadius")),
0160       a_mipR_(iConfig.getParameter<double>("coneRadiusMIP")),
0161       maxRestrictionP_(iConfig.getParameter<double>("maxTrackP")),
0162       slopeRestrictionP_(iConfig.getParameter<double>("slopeTrackP")),
0163       pTrackMin_(iConfig.getParameter<double>("minimumTrackP")),
0164       eEcalMax_(iConfig.getParameter<double>("maximumEcalEnergy")),
0165       eIsolate_(iConfig.getParameter<double>("isolationEnergy")),
0166       pTrackLow_(iConfig.getParameter<double>("momentumRangeLow")),
0167       pTrackHigh_(iConfig.getParameter<double>("momentumRangeHigh")),
0168       preScale_(iConfig.getParameter<int>("preScaleFactor")),
0169       labelGenTrack_(iConfig.getParameter<edm::InputTag>("TrackLabel")),
0170       labelRecVtx_(iConfig.getParameter<edm::InputTag>("VertexLabel")),
0171       labelBS_(iConfig.getParameter<edm::InputTag>("BeamSpotLabel")),
0172       labelEB_(iConfig.getParameter<edm::InputTag>("EBRecHitLabel")),
0173       labelEE_(iConfig.getParameter<edm::InputTag>("EERecHitLabel")),
0174       labelHBHE_(iConfig.getParameter<edm::InputTag>("HBHERecHitLabel")),
0175       labelHltGT_(iConfig.getParameter<edm::InputTag>("L1GTSeedLabel")),
0176       labelTriggerEvent_(iConfig.getParameter<edm::InputTag>("TriggerEventLabel")),
0177       labelTriggerResults_(iConfig.getParameter<edm::InputTag>("TriggerResultLabel")),
0178       labelIsoTk_(iConfig.getParameter<std::string>("IsoTrackLabel")) {
0179   // Get the run parameters
0180   // Different isolation cuts are described in DN-2016/029
0181   // Tight cut uses 2 GeV; Loose cut uses 10 GeV
0182   // Eta dependent cut uses (maxRestrictionP_ * exp(|ieta|*log(2.5)/18))
0183   // with the factor for exponential slopeRestrictionP_ = log(2.5)/18
0184   // maxRestrictionP_ = 8 GeV as came from a study
0185   const double isolationRadius(28.9);
0186   selectionParameter_.minPt = iConfig.getParameter<double>("minTrackPt");
0187   ;
0188   selectionParameter_.minQuality = reco::TrackBase::qualityByName(theTrackQuality_);
0189   selectionParameter_.maxDxyPV = iConfig.getParameter<double>("maxDxyPV");
0190   selectionParameter_.maxDzPV = iConfig.getParameter<double>("maxDzPV");
0191   selectionParameter_.maxChi2 = iConfig.getParameter<double>("maxChi2");
0192   selectionParameter_.maxDpOverP = iConfig.getParameter<double>("maxDpOverP");
0193   selectionParameter_.minOuterHit = iConfig.getParameter<int>("minOuterHit");
0194   selectionParameter_.minLayerCrossed = iConfig.getParameter<int>("minLayerCrossed");
0195   selectionParameter_.maxInMiss = iConfig.getParameter<int>("maxInMiss");
0196   selectionParameter_.maxOutMiss = iConfig.getParameter<int>("maxOutMiss");
0197   a_charIsoR_ = a_coneR_ + isolationRadius;
0198
0199   // define tokens for access
0200   tok_hltGT_ = consumes<trigger::TriggerFilterObjectWithRefs>(labelHltGT_);
0201   tok_trigEvt_ = consumes<trigger::TriggerEvent>(labelTriggerEvent_);
0202   tok_trigRes_ = consumes<edm::TriggerResults>(labelTriggerResults_);
0203   tok_genTrack_ = consumes<reco::TrackCollection>(labelGenTrack_);
0204   tok_recVtx_ = consumes<reco::VertexCollection>(labelRecVtx_);
0205   tok_bs_ = consumes<reco::BeamSpot>(labelBS_);
0206   tok_EB_ = consumes<EcalRecHitCollection>(labelEB_);
0207   tok_EE_ = consumes<EcalRecHitCollection>(labelEE_);
0208   tok_hbhe_ = consumes<HBHERecHitCollection>(labelHBHE_);
0209   // for event setup
0210   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>();
0211   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>();
0212
0213   edm::LogVerbatim("HcalIsoTrack") << "Parameters read from config file \n"
0214                                    << "\t minPt " << selectionParameter_.minPt << "\t theTrackQuality "
0215                                    << theTrackQuality_ << "\t minQuality " << selectionParameter_.minQuality
0216                                    << "\t maxDxyPV " << selectionParameter_.maxDxyPV << "\t maxDzPV "
0217                                    << selectionParameter_.maxDzPV << "\t maxChi2 " << selectionParameter_.maxChi2
0218                                    << "\t maxDpOverP " << selectionParameter_.maxDpOverP << "\t minOuterHit "
0219                                    << selectionParameter_.minOuterHit << "\t minLayerCrossed "
0220                                    << selectionParameter_.minLayerCrossed << "\t maxInMiss "
0221                                    << selectionParameter_.maxInMiss << "\t maxOutMiss "
0222                                    << selectionParameter_.maxOutMiss << "\n"
0223                                    << "\t a_coneR " << a_coneR_ << "\t a_charIsoR " << a_charIsoR_ << "\t a_mipR "
0224                                    << a_mipR_ << "\t pTrackMin " << pTrackMin_ << "\t eEcalMax " << eEcalMax_
0225                                    << "\t maxRestrictionP_ " << maxRestrictionP_ << "\t slopeRestrictionP_ "
0226                                    << slopeRestrictionP_ << "\t eIsolate_ " << eIsolate_ << "\t Process "
0227                                    << processName_ << "\n"
0228                                    << "\t Precale factor " << preScale_ << "\t in momentum range " << pTrackLow_ << ":"
0229                                    << pTrackHigh_;
0230   for (unsigned int k = 0; k < trigNames_.size(); ++k)
0231     edm::LogVerbatim("HcalIsoTrack") << "Trigger[" << k << "] " << trigNames_[k];
0232
0233   //create also IsolatedPixelTrackCandidateCollection which contains isolation info and reference to primary track
0234   produces<reco::HcalIsolatedTrackCandidateCollection>(labelIsoTk_);
0235   produces<reco::VertexCollection>(labelRecVtx_.label());
0236   produces<EcalRecHitCollection>(labelEB_.instance());
0237   produces<EcalRecHitCollection>(labelEE_.instance());
0238   produces<HBHERecHitCollection>(labelHBHE_.label());
0239
0240   edm::LogVerbatim("HcalIsoTrack") << " Expected to produce the collections:\n"
0241                                    << "reco::HcalIsolatedTrackCandidateCollection "
0242                                    << " with label HcalIsolatedTrackCollection\n"
0243                                    << "reco::VertexCollection with label " << labelRecVtx_.label() << "\n"
0244                                    << "EcalRecHitCollection with label EcalRecHitsEB\n"
0245                                    << "EcalRecHitCollection with label EcalRecHitsEE\n"
0246                                    << "HBHERecHitCollection with label " << labelHBHE_.label();
0247 }
0248
0249 void AlCaIsoTracksProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0250   edm::ParameterSetDescription desc;
0251   // producer for  (HCAL isolated tracks)
0252   desc.add<edm::InputTag>("TrackLabel", edm::InputTag("generalTracks"));
0253   desc.add<edm::InputTag>("VertexLabel", edm::InputTag("offlinePrimaryVertices"));
0254   desc.add<edm::InputTag>("BeamSpotLabel", edm::InputTag("offlineBeamSpot"));
0255   desc.add<edm::InputTag>("EBRecHitLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0256   desc.add<edm::InputTag>("EERecHitLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0257   desc.add<edm::InputTag>("HBHERecHitLabel", edm::InputTag("hbhereco"));
0258   desc.add<edm::InputTag>("L1GTSeedLabel", edm::InputTag("hltL1sV0SingleJet60"));
0259   desc.add<edm::InputTag>("TriggerEventLabel", edm::InputTag("hltTriggerSummaryAOD", "", "HLT"));
0260   desc.add<edm::InputTag>("TriggerResultLabel", edm::InputTag("TriggerResults", "", "HLT"));
0261   desc.add<std::string>("IsoTrackLabel", "HcalIsolatedTrackCollection");
0262   std::vector<std::string> triggers = {"HLT_IsoTrackHB", "HLT_IsoTrackHE"};
0263   desc.add<std::vector<std::string> >("triggers", triggers);
0264   desc.add<std::string>("processName", "HLT");
0265   // following 10 parameters are parameters to select good tracks
0266   desc.add<std::string>("trackQuality", "highPurity");
0267   desc.add<double>("minTrackPt", 1.0);
0268   desc.add<double>("maxDxyPV", 10.0);
0269   desc.add<double>("maxDzPV", 100.0);
0270   desc.add<double>("maxChi2", 5.0);
0271   desc.add<double>("maxDpOverP", 0.1);
0272   desc.add<int>("minOuterHit", 4);
0273   desc.add<int>("minLayerCrossed", 8);
0274   desc.add<int>("maxInMiss", 2);
0275   desc.add<int>("maxOutMiss", 2);
0276   // Minimum momentum of selected isolated track and signal zone
0277   desc.add<double>("coneRadius", 34.98);
0278   desc.add<double>("minimumTrackP", 20.0);
0279   // signal zone in ECAL and MIP energy cutoff
0280   desc.add<double>("coneRadiusMIP", 14.0);
0281   desc.add<double>("maximumEcalEnergy", 2.0);
0282   // following 3 parameters are for isolation cuts and described in the code
0283   desc.add<double>("maxTrackP", 8.0);
0284   desc.add<double>("slopeTrackP", 0.05090504066);
0285   desc.add<double>("isolationEnergy", 10.0);
0286   // Prescale events only containing isolated tracks in the range
0287   desc.add<double>("momentumRangeLow", 20.0);
0288   desc.add<double>("momentumRangeHigh", 40.0);
0289   desc.add<int>("preScaleFactor", 10);
0290   descriptions.add("alcaisotrk", desc);
0291 }
0292
0293 void AlCaIsoTracksProducer::produce(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0294   nAll_++;
0295 #ifdef EDM_ML_DEBUG
0296   edm::LogVerbatim("HcalIsoTrack") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event() << " Luminosity "
0297                                    << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing();
0298 #endif
0299   bool valid(true);
0300   //Step1: Get all the relevant containers
0301   trigger::TriggerEvent triggerEvent;
0302   auto const& triggerEventHandle = iEvent.getHandle(tok_trigEvt_);
0303   if (!triggerEventHandle.isValid()) {
0304     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelTriggerEvent_;
0305     valid = false;
0306   }
0307   auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
0308   if (!triggerResults.isValid()) {
0309     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelTriggerResults_;
0310     valid = false;
0311   }
0312
0313   auto trkCollection = iEvent.getHandle(tok_genTrack_);
0314   if (!trkCollection.isValid()) {
0315     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
0316     valid = false;
0317   }
0318
0319   auto const& recVtxs = iEvent.getHandle(tok_recVtx_);
0320   if (!recVtxs.isValid()) {
0321     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
0322     valid = false;
0323   }
0324
0325   auto const& beamSpotH = iEvent.getHandle(tok_bs_);
0326   math::XYZPoint leadPV(0, 0, 0);
0327   if (valid) {
0328     if (!recVtxs->empty() && !((*recVtxs)[0].isFake())) {
0329       leadPV = math::XYZPoint((*recVtxs)[0].x(), (*recVtxs)[0].y(), (*recVtxs)[0].z());
0330     } else if (beamSpotH.isValid()) {
0331       leadPV = beamSpotH->position();
0332     }
0333   }
0334 #ifdef EDM_ML_DEBUG
0335   edm::LogVerbatim("HcalIsoTrack") << "Primary Vertex " << leadPV;
0336 #endif
0337
0338   auto barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0339   if (!barrelRecHitsHandle.isValid()) {
0340     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEB_;
0341     valid = false;
0342   }
0343   auto endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0344   if (!endcapRecHitsHandle.isValid()) {
0345     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEE_;
0346     valid = false;
0347   }
0348   auto hbhe = iEvent.getHandle(tok_hbhe_);
0349   if (!hbhe.isValid()) {
0350     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelHBHE_;
0351     valid = false;
0352   }
0353
0354   //Get L1 trigger object
0355   double ptL1(0), etaL1(0), phiL1(0);
0356   auto const& l1trigobj = iEvent.getHandle(tok_hltGT_);
0357
0358   if (l1trigobj.isValid()) {
0359     std::vector<edm::Ref<l1extra::L1JetParticleCollection> > l1tauobjref;
0360     l1trigobj->getObjects(trigger::TriggerL1TauJet, l1tauobjref);
0361     setPtEtaPhi(l1tauobjref, ptL1, etaL1, phiL1);
0362
0363     std::vector<edm::Ref<l1extra::L1JetParticleCollection> > l1jetobjref;
0364     l1trigobj->getObjects(trigger::TriggerL1CenJet, l1jetobjref);
0365     setPtEtaPhi(l1jetobjref, ptL1, etaL1, phiL1);
0366
0367     std::vector<edm::Ref<l1extra::L1JetParticleCollection> > l1forjetobjref;
0368     l1trigobj->getObjects(trigger::TriggerL1ForJet, l1forjetobjref);
0369     setPtEtaPhi(l1forjetobjref, ptL1, etaL1, phiL1);
0370   } else {
0371     valid = false;
0372   }
0373
0374   auto outputHcalIsoTrackColl = std::make_unique<reco::HcalIsolatedTrackCandidateCollection>();
0375   auto outputVColl = std::make_unique<reco::VertexCollection>();
0376   auto outputEBColl = std::make_unique<EBRecHitCollection>();
0377   auto outputEEColl = std::make_unique<EERecHitCollection>();
0378   auto outputHBHEColl = std::make_unique<HBHERecHitCollection>();
0379
0380   //For valid HLT record
0381   if (!valid) {
0382     edm::LogWarning("HcalIsoTrack") << "Error! Can't get some of the products";
0383   } else {
0384     trigger::TriggerEvent triggerEvent = *(triggerEventHandle.product());
0385     if (triggerResults.isValid()) {
0386       const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0387       const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
0388       reco::HcalIsolatedTrackCandidateCollection* isotk = select(triggerResults,
0389                                                                  triggerNames_,
0390                                                                  trkCollection,
0391                                                                  leadPV,
0392                                                                  barrelRecHitsHandle,
0393                                                                  endcapRecHitsHandle,
0394                                                                  hbhe,
0395                                                                  ptL1,
0396                                                                  etaL1,
0397                                                                  phiL1);
0398 #ifdef EDM_ML_DEBUG
0399       edm::LogVerbatim("HcalIsoTrack") << "AlCaIsoTracksProducer::select returns " << isotk->size()
0400                                        << " isolated tracks";
0401 #endif
0402
0403       if (!isotk->empty()) {
0404         int ntrin(0), ntrout(0);
0405         for (reco::HcalIsolatedTrackCandidateCollection::const_iterator itr = isotk->begin(); itr != isotk->end();
0406              ++itr) {
0407           if (itr->p() > pTrackLow_ && itr->p() < pTrackHigh_)
0408             ntrin++;
0409           else
0410             ntrout++;
0411         }
0412         bool selectEvent = ntrout > 0;
0413         if (!selectEvent && ntrin > 0) {
0414           ++nRange_;
0415           if (preScale_ <= 1)
0416             selectEvent = true;
0417           else if (nRange_ % preScale_ == 1)
0418             selectEvent = true;
0419         }
0420         if (selectEvent) {
0421           for (reco::HcalIsolatedTrackCandidateCollection::const_iterator itr = isotk->begin(); itr != isotk->end();
0422                ++itr)
0423             outputHcalIsoTrackColl->push_back(*itr);
0424
0425           for (reco::VertexCollection::const_iterator vtx = recVtxs->begin(); vtx != recVtxs->end(); ++vtx)
0426             outputVColl->push_back(*vtx);
0427
0428           for (edm::SortedCollection<EcalRecHit>::const_iterator ehit = barrelRecHitsHandle->begin();
0429                ehit != barrelRecHitsHandle->end();
0430                ++ehit)
0431             outputEBColl->push_back(*ehit);
0432
0433           for (edm::SortedCollection<EcalRecHit>::const_iterator ehit = endcapRecHitsHandle->begin();
0434                ehit != endcapRecHitsHandle->end();
0435                ++ehit)
0436             outputEEColl->push_back(*ehit);
0437
0438           for (std::vector<HBHERecHit>::const_iterator hhit = hbhe->begin(); hhit != hbhe->end(); ++hhit)
0439             outputHBHEColl->push_back(*hhit);
0440           ++nGood_;
0441         }
0442       }
0443     }
0444   }
0445   iEvent.put(std::move(outputHcalIsoTrackColl), labelIsoTk_);
0446   iEvent.put(std::move(outputVColl), labelRecVtx_.label());
0447   iEvent.put(std::move(outputEBColl), labelEB_.instance());
0448   iEvent.put(std::move(outputEEColl), labelEE_.instance());
0449   iEvent.put(std::move(outputHBHEColl), labelHBHE_.label());
0450 }
0451
0452 void AlCaIsoTracksProducer::endStream() {
0453   globalCache()->nAll_ += nAll_;
0454   globalCache()->nGood_ += nGood_;
0455   globalCache()->nRange_ += nRange_;
0456 }
0457
0458 void AlCaIsoTracksProducer::globalEndJob(const alCaIsoTracksProducer::Counters* count) {
0459   edm::LogVerbatim("HcalIsoTrack") << "Finds " << count->nGood_ << " good tracks in " << count->nAll_ << " events and "
0460                                    << count->nRange_ << " events in the momentum raange";
0461 }
0462
0463 void AlCaIsoTracksProducer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0464   bool changed(false);
0465   edm::LogVerbatim("HcalIsoTrack") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init "
0466                                    << hltConfig_.init(iRun, iSetup, processName_, changed);
0467
0468   bField = &(iSetup.getData(tok_magField_));
0469   geo = &(iSetup.getData(tok_geom_));
0470 }
0471
0472 void AlCaIsoTracksProducer::endRun(edm::Run const& iRun, edm::EventSetup const&) {
0473   edm::LogVerbatim("HcalIsoTrack") << "endRun [" << nRun_ << "] " << iRun.run();
0474   ++nRun_;
0475 }
0476
0477 reco::HcalIsolatedTrackCandidateCollection* AlCaIsoTracksProducer::select(
0478     edm::Handle<edm::TriggerResults> const& triggerResults,
0479     const std::vector<std::string>& triggerNames_,
0480     edm::Handle<reco::TrackCollection>& trkCollection,
0481     math::XYZPoint& leadPV,
0482     edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
0483     edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
0484     edm::Handle<HBHERecHitCollection>& hbhe,
0485     double ptL1,
0486     double etaL1,
0487     double phiL1) {
0488   reco::HcalIsolatedTrackCandidateCollection* trackCollection = new reco::HcalIsolatedTrackCandidateCollection;
0489   bool ok(false);
0490
0491   // Find a good HLT trigger
0492   for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0493     int hlt = triggerResults->accept(iHLT);
0494     for (unsigned int i = 0; i < trigNames_.size(); ++i) {
0495       if (triggerNames_[iHLT].find(trigNames_[i]) != std::string::npos) {
0496         if (hlt > 0) {
0497           ok = true;
0498         }
0499         edm::LogVerbatim("HcalIsoTrack") << "The trigger we are looking for " << triggerNames_[iHLT] << " Flag " << hlt
0500                                          << ":" << ok;
0501       }
0502     }
0503   }
0504
0505   //Propagate tracks to calorimeter surface)
0506   std::vector<spr::propagatedTrackDirection> trkCaloDirections;
0507   spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, false);
0508
0509   std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
0510   unsigned int nTracks(0);
0511   for (trkDetItr = trkCaloDirections.begin(), nTracks = 0; trkDetItr != trkCaloDirections.end();
0512        trkDetItr++, nTracks++) {
0513     const reco::Track* pTrack = &(*(trkDetItr->trkItr));
0514     math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
0515 #ifdef EDM_ML_DEBUG
0516     edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) :" << pTrack->pt() << "|"
0517                                      << pTrack->eta() << "|" << pTrack->phi() << "|" << pTrack->p();
0518 #endif
0519     //Selection of good track
0520     bool qltyFlag = spr::goodTrack(pTrack, leadPV, selectionParameter_, false);
0521 #ifdef EDM_ML_DEBUG
0522     edm::LogVerbatim("HcalIsoTrack") << "qltyFlag|okECAL|okHCAL : " << qltyFlag << "|" << trkDetItr->okECAL << "|"
0523                                      << trkDetItr->okHCAL;
0524 #endif
0525     if (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL) {
0526       double t_p = pTrack->p();
0527       int nRH_eMipDR(0), nNearTRKs(0);
0528       double eMipDR = spr::eCone_ecal(geo,
0529                                       barrelRecHitsHandle,
0530                                       endcapRecHitsHandle,
0531                                       trkDetItr->pointHCAL,
0532                                       trkDetItr->pointECAL,
0533                                       a_mipR_,
0534                                       trkDetItr->directionECAL,
0535                                       nRH_eMipDR);
0536       double hmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, false);
0537       HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
0538       int ieta = detId.ietaAbs();
0539       double eIsolation = (maxRestrictionP_ * exp(slopeRestrictionP_ * ((double)(ieta))));
0540       if (eIsolation < eIsolate_)
0541         eIsolation = eIsolate_;
0542 #ifdef EDM_ML_DEBUG
0543       edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) :" << pTrack->pt() << "|"
0544                                        << pTrack->eta() << "|" << pTrack->phi() << "|" << t_p << " e_MIP " << eMipDR
0545                                        << " Chg Isolation " << hmaxNearP << ":" << eIsolation;
0546 #endif
0547       if (t_p > pTrackMin_ && eMipDR < eEcalMax_ && hmaxNearP < eIsolation) {
0548         reco::HcalIsolatedTrackCandidate newCandidate(v4);
0549         newCandidate.SetMaxP(hmaxNearP);
0550         newCandidate.SetEnergyEcal(eMipDR);
0551         newCandidate.setL1(ptL1, etaL1, phiL1);
0552         newCandidate.SetEtaPhiEcal((trkDetItr->pointECAL).eta(), (trkDetItr->pointECAL).phi());
0553         HcalDetId detId = HcalDetId(trkDetItr->detIdHCAL);
0554         newCandidate.SetEtaPhiHcal(
0555             (trkDetItr->pointHCAL).eta(), (trkDetItr->pointHCAL).phi(), detId.ieta(), detId.iphi());
0556         int indx(0);
0557         for (reco::TrackCollection::const_iterator trkItr1 = trkCollection->begin(); trkItr1 != trkCollection->end();
0558              ++trkItr1, ++indx) {
0559           const reco::Track* pTrack1 = &(*trkItr1);
0560           if (pTrack1 == pTrack) {
0561             reco::TrackRef tRef = reco::TrackRef(trkCollection, indx);
0562             newCandidate.setTrack(tRef);
0563             break;
0564           }
0565         }
0566         trackCollection->push_back(newCandidate);
0567       }
0568     }
0569   }
0570   return trackCollection;
0571 }
0572
0573 void AlCaIsoTracksProducer::setPtEtaPhi(std::vector<edm::Ref<l1extra::L1JetParticleCollection> >& objref,
0574                                         double& ptL1,
0575                                         double& etaL1,
0576                                         double& phiL1) {
0577   for (unsigned int p = 0; p < objref.size(); p++) {
0578     if (objref[p]->pt() > ptL1) {
0579       ptL1 = objref[p]->pt();
0580       phiL1 = objref[p]->phi();
0581       etaL1 = objref[p]->eta();
0582     }
0583   }
0584 }
0585
0586 #include "FWCore/Framework/interface/MakerMacros.h"
0587
0588 DEFINE_FWK_MODULE(AlCaIsoTracksProducer);