File indexing completed on 2021-10-01 22:40:16
0001
0002
0003
0004
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
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
0075
0076
0077
0078
0079 namespace AlCaIsoTracks {
0080 struct Counters {
0081 Counters() : nAll_(0), nGood_(0), nRange_(0) {}
0082 mutable std::atomic<unsigned int> nAll_, nGood_, nRange_;
0083 };
0084 }
0085
0086 class AlCaIsoTracksProducer : public edm::stream::EDProducer<edm::GlobalCache<AlCaIsoTracks::Counters> > {
0087 public:
0088 explicit AlCaIsoTracksProducer(edm::ParameterSet const&, const AlCaIsoTracks::Counters* count);
0089 ~AlCaIsoTracksProducer() override;
0090
0091 static std::unique_ptr<AlCaIsoTracks::Counters> initializeGlobalCache(edm::ParameterSet const&) {
0092 return std::make_unique<AlCaIsoTracks::Counters>();
0093 }
0094
0095 void produce(edm::Event&, edm::EventSetup const&) override;
0096 void endStream() override;
0097 static void globalEndJob(const AlCaIsoTracks::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
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, const AlCaIsoTracks::Counters* counters)
0151 : nRun_(0),
0152 nAll_(0),
0153 nGood_(0),
0154 nRange_(0),
0155 trigNames_(iConfig.getParameter<std::vector<std::string> >("triggers")),
0156 theTrackQuality_(iConfig.getParameter<std::string>("trackQuality")),
0157 processName_(iConfig.getParameter<std::string>("processName")),
0158 a_coneR_(iConfig.getParameter<double>("coneRadius")),
0159 a_mipR_(iConfig.getParameter<double>("coneRadiusMIP")),
0160 maxRestrictionP_(iConfig.getParameter<double>("maxTrackP")),
0161 slopeRestrictionP_(iConfig.getParameter<double>("slopeTrackP")),
0162 pTrackMin_(iConfig.getParameter<double>("minimumTrackP")),
0163 eEcalMax_(iConfig.getParameter<double>("maximumEcalEnergy")),
0164 eIsolate_(iConfig.getParameter<double>("isolationEnergy")),
0165 pTrackLow_(iConfig.getParameter<double>("momentumRangeLow")),
0166 pTrackHigh_(iConfig.getParameter<double>("momentumRangeHigh")),
0167 preScale_(iConfig.getParameter<int>("preScaleFactor")),
0168 labelGenTrack_(iConfig.getParameter<edm::InputTag>("TrackLabel")),
0169 labelRecVtx_(iConfig.getParameter<edm::InputTag>("VertexLabel")),
0170 labelBS_(iConfig.getParameter<edm::InputTag>("BeamSpotLabel")),
0171 labelEB_(iConfig.getParameter<edm::InputTag>("EBRecHitLabel")),
0172 labelEE_(iConfig.getParameter<edm::InputTag>("EERecHitLabel")),
0173 labelHBHE_(iConfig.getParameter<edm::InputTag>("HBHERecHitLabel")),
0174 labelHltGT_(iConfig.getParameter<edm::InputTag>("L1GTSeedLabel")),
0175 labelTriggerEvent_(iConfig.getParameter<edm::InputTag>("TriggerEventLabel")),
0176 labelTriggerResults_(iConfig.getParameter<edm::InputTag>("TriggerResultLabel")),
0177 labelIsoTk_(iConfig.getParameter<std::string>("IsoTrackLabel")) {
0178
0179
0180
0181
0182
0183
0184 const double isolationRadius(28.9);
0185 selectionParameter_.minPt = iConfig.getParameter<double>("minTrackPt");
0186 ;
0187 selectionParameter_.minQuality = reco::TrackBase::qualityByName(theTrackQuality_);
0188 selectionParameter_.maxDxyPV = iConfig.getParameter<double>("maxDxyPV");
0189 selectionParameter_.maxDzPV = iConfig.getParameter<double>("maxDzPV");
0190 selectionParameter_.maxChi2 = iConfig.getParameter<double>("maxChi2");
0191 selectionParameter_.maxDpOverP = iConfig.getParameter<double>("maxDpOverP");
0192 selectionParameter_.minOuterHit = iConfig.getParameter<int>("minOuterHit");
0193 selectionParameter_.minLayerCrossed = iConfig.getParameter<int>("minLayerCrossed");
0194 selectionParameter_.maxInMiss = iConfig.getParameter<int>("maxInMiss");
0195 selectionParameter_.maxOutMiss = iConfig.getParameter<int>("maxOutMiss");
0196 a_charIsoR_ = a_coneR_ + isolationRadius;
0197
0198
0199 tok_hltGT_ = consumes<trigger::TriggerFilterObjectWithRefs>(labelHltGT_);
0200 tok_trigEvt_ = consumes<trigger::TriggerEvent>(labelTriggerEvent_);
0201 tok_trigRes_ = consumes<edm::TriggerResults>(labelTriggerResults_);
0202 tok_genTrack_ = consumes<reco::TrackCollection>(labelGenTrack_);
0203 tok_recVtx_ = consumes<reco::VertexCollection>(labelRecVtx_);
0204 tok_bs_ = consumes<reco::BeamSpot>(labelBS_);
0205 tok_EB_ = consumes<EcalRecHitCollection>(labelEB_);
0206 tok_EE_ = consumes<EcalRecHitCollection>(labelEE_);
0207 tok_hbhe_ = consumes<HBHERecHitCollection>(labelHBHE_);
0208
0209 tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>();
0210 tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>();
0211
0212 edm::LogVerbatim("HcalIsoTrack") << "Parameters read from config file \n"
0213 << "\t minPt " << selectionParameter_.minPt << "\t theTrackQuality "
0214 << theTrackQuality_ << "\t minQuality " << selectionParameter_.minQuality
0215 << "\t maxDxyPV " << selectionParameter_.maxDxyPV << "\t maxDzPV "
0216 << selectionParameter_.maxDzPV << "\t maxChi2 " << selectionParameter_.maxChi2
0217 << "\t maxDpOverP " << selectionParameter_.maxDpOverP << "\t minOuterHit "
0218 << selectionParameter_.minOuterHit << "\t minLayerCrossed "
0219 << selectionParameter_.minLayerCrossed << "\t maxInMiss "
0220 << selectionParameter_.maxInMiss << "\t maxOutMiss "
0221 << selectionParameter_.maxOutMiss << "\n"
0222 << "\t a_coneR " << a_coneR_ << "\t a_charIsoR " << a_charIsoR_ << "\t a_mipR "
0223 << a_mipR_ << "\t pTrackMin " << pTrackMin_ << "\t eEcalMax " << eEcalMax_
0224 << "\t maxRestrictionP_ " << maxRestrictionP_ << "\t slopeRestrictionP_ "
0225 << slopeRestrictionP_ << "\t eIsolate_ " << eIsolate_ << "\t Process "
0226 << processName_ << "\n"
0227 << "\t Precale factor " << preScale_ << "\t in momentum range " << pTrackLow_ << ":"
0228 << pTrackHigh_;
0229 for (unsigned int k = 0; k < trigNames_.size(); ++k)
0230 edm::LogVerbatim("HcalIsoTrack") << "Trigger[" << k << "] " << trigNames_[k];
0231
0232
0233 produces<reco::HcalIsolatedTrackCandidateCollection>(labelIsoTk_);
0234 produces<reco::VertexCollection>(labelRecVtx_.label());
0235 produces<EcalRecHitCollection>(labelEB_.instance());
0236 produces<EcalRecHitCollection>(labelEE_.instance());
0237 produces<HBHERecHitCollection>(labelHBHE_.label());
0238
0239 edm::LogVerbatim("HcalIsoTrack") << " Expected to produce the collections:\n"
0240 << "reco::HcalIsolatedTrackCandidateCollection "
0241 << " with label HcalIsolatedTrackCollection\n"
0242 << "reco::VertexCollection with label " << labelRecVtx_.label() << "\n"
0243 << "EcalRecHitCollection with label EcalRecHitsEB\n"
0244 << "EcalRecHitCollection with label EcalRecHitsEE\n"
0245 << "HBHERecHitCollection with label " << labelHBHE_.label();
0246 }
0247
0248 AlCaIsoTracksProducer::~AlCaIsoTracksProducer() {}
0249
0250 void AlCaIsoTracksProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0251 edm::ParameterSetDescription desc;
0252
0253 desc.add<edm::InputTag>("TrackLabel", edm::InputTag("generalTracks"));
0254 desc.add<edm::InputTag>("VertexLabel", edm::InputTag("offlinePrimaryVertices"));
0255 desc.add<edm::InputTag>("BeamSpotLabel", edm::InputTag("offlineBeamSpot"));
0256 desc.add<edm::InputTag>("EBRecHitLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0257 desc.add<edm::InputTag>("EERecHitLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0258 desc.add<edm::InputTag>("HBHERecHitLabel", edm::InputTag("hbhereco"));
0259 desc.add<edm::InputTag>("L1GTSeedLabel", edm::InputTag("hltL1sV0SingleJet60"));
0260 desc.add<edm::InputTag>("TriggerEventLabel", edm::InputTag("hltTriggerSummaryAOD", "", "HLT"));
0261 desc.add<edm::InputTag>("TriggerResultLabel", edm::InputTag("TriggerResults", "", "HLT"));
0262 desc.add<std::string>("IsoTrackLabel", "HcalIsolatedTrackCollection");
0263 std::vector<std::string> triggers = {"HLT_IsoTrackHB", "HLT_IsoTrackHE"};
0264 desc.add<std::vector<std::string> >("triggers", triggers);
0265 desc.add<std::string>("processName", "HLT");
0266
0267 desc.add<std::string>("trackQuality", "highPurity");
0268 desc.add<double>("minTrackPt", 1.0);
0269 desc.add<double>("maxDxyPV", 10.0);
0270 desc.add<double>("maxDzPV", 100.0);
0271 desc.add<double>("maxChi2", 5.0);
0272 desc.add<double>("maxDpOverP", 0.1);
0273 desc.add<int>("minOuterHit", 4);
0274 desc.add<int>("minLayerCrossed", 8);
0275 desc.add<int>("maxInMiss", 2);
0276 desc.add<int>("maxOutMiss", 2);
0277
0278 desc.add<double>("coneRadius", 34.98);
0279 desc.add<double>("minimumTrackP", 20.0);
0280
0281 desc.add<double>("coneRadiusMIP", 14.0);
0282 desc.add<double>("maximumEcalEnergy", 2.0);
0283
0284 desc.add<double>("maxTrackP", 8.0);
0285 desc.add<double>("slopeTrackP", 0.05090504066);
0286 desc.add<double>("isolationEnergy", 10.0);
0287
0288 desc.add<double>("momentumRangeLow", 20.0);
0289 desc.add<double>("momentumRangeHigh", 40.0);
0290 desc.add<int>("preScaleFactor", 10);
0291 descriptions.add("alcaisotrk", desc);
0292 }
0293
0294 void AlCaIsoTracksProducer::produce(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0295 nAll_++;
0296 #ifdef EDM_ML_DEBUG
0297 edm::LogVerbatim("HcalIsoTrack") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event() << " Luminosity "
0298 << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing();
0299 #endif
0300 bool valid(true);
0301
0302 trigger::TriggerEvent triggerEvent;
0303 auto const& triggerEventHandle = iEvent.getHandle(tok_trigEvt_);
0304 if (!triggerEventHandle.isValid()) {
0305 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelTriggerEvent_;
0306 valid = false;
0307 }
0308 auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
0309 if (!triggerResults.isValid()) {
0310 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelTriggerResults_;
0311 valid = false;
0312 }
0313
0314 auto trkCollection = iEvent.getHandle(tok_genTrack_);
0315 if (!trkCollection.isValid()) {
0316 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
0317 valid = false;
0318 }
0319
0320 auto const& recVtxs = iEvent.getHandle(tok_recVtx_);
0321 if (!recVtxs.isValid()) {
0322 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
0323 valid = false;
0324 }
0325
0326 auto const& beamSpotH = iEvent.getHandle(tok_bs_);
0327 math::XYZPoint leadPV(0, 0, 0);
0328 if (valid) {
0329 if (!recVtxs->empty() && !((*recVtxs)[0].isFake())) {
0330 leadPV = math::XYZPoint((*recVtxs)[0].x(), (*recVtxs)[0].y(), (*recVtxs)[0].z());
0331 } else if (beamSpotH.isValid()) {
0332 leadPV = beamSpotH->position();
0333 }
0334 }
0335 #ifdef EDM_ML_DEBUG
0336 edm::LogVerbatim("HcalIsoTrack") << "Primary Vertex " << leadPV;
0337 #endif
0338
0339 auto barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0340 if (!barrelRecHitsHandle.isValid()) {
0341 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEB_;
0342 valid = false;
0343 }
0344 auto endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0345 if (!endcapRecHitsHandle.isValid()) {
0346 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEE_;
0347 valid = false;
0348 }
0349 auto hbhe = iEvent.getHandle(tok_hbhe_);
0350 if (!hbhe.isValid()) {
0351 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelHBHE_;
0352 valid = false;
0353 }
0354
0355
0356 double ptL1(0), etaL1(0), phiL1(0);
0357 auto const& l1trigobj = iEvent.getHandle(tok_hltGT_);
0358
0359 if (l1trigobj.isValid()) {
0360 std::vector<edm::Ref<l1extra::L1JetParticleCollection> > l1tauobjref;
0361 l1trigobj->getObjects(trigger::TriggerL1TauJet, l1tauobjref);
0362 setPtEtaPhi(l1tauobjref, ptL1, etaL1, phiL1);
0363
0364 std::vector<edm::Ref<l1extra::L1JetParticleCollection> > l1jetobjref;
0365 l1trigobj->getObjects(trigger::TriggerL1CenJet, l1jetobjref);
0366 setPtEtaPhi(l1jetobjref, ptL1, etaL1, phiL1);
0367
0368 std::vector<edm::Ref<l1extra::L1JetParticleCollection> > l1forjetobjref;
0369 l1trigobj->getObjects(trigger::TriggerL1ForJet, l1forjetobjref);
0370 setPtEtaPhi(l1forjetobjref, ptL1, etaL1, phiL1);
0371 } else {
0372 valid = false;
0373 }
0374
0375 auto outputHcalIsoTrackColl = std::make_unique<reco::HcalIsolatedTrackCandidateCollection>();
0376 auto outputVColl = std::make_unique<reco::VertexCollection>();
0377 auto outputEBColl = std::make_unique<EBRecHitCollection>();
0378 auto outputEEColl = std::make_unique<EERecHitCollection>();
0379 auto outputHBHEColl = std::make_unique<HBHERecHitCollection>();
0380
0381
0382 if (!valid) {
0383 edm::LogWarning("HcalIsoTrack") << "Error! Can't get some of the products";
0384 } else {
0385 trigger::TriggerEvent triggerEvent = *(triggerEventHandle.product());
0386 if (triggerResults.isValid()) {
0387 const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0388 const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
0389 reco::HcalIsolatedTrackCandidateCollection* isotk = select(triggerResults,
0390 triggerNames_,
0391 trkCollection,
0392 leadPV,
0393 barrelRecHitsHandle,
0394 endcapRecHitsHandle,
0395 hbhe,
0396 ptL1,
0397 etaL1,
0398 phiL1);
0399 #ifdef EDM_ML_DEBUG
0400 edm::LogVerbatim("HcalIsoTrack") << "AlCaIsoTracksProducer::select returns " << isotk->size()
0401 << " isolated tracks";
0402 #endif
0403
0404 if (!isotk->empty()) {
0405 int ntrin(0), ntrout(0);
0406 for (reco::HcalIsolatedTrackCandidateCollection::const_iterator itr = isotk->begin(); itr != isotk->end();
0407 ++itr) {
0408 if (itr->p() > pTrackLow_ && itr->p() < pTrackHigh_)
0409 ntrin++;
0410 else
0411 ntrout++;
0412 }
0413 bool selectEvent = ntrout > 0;
0414 if (!selectEvent && ntrin > 0) {
0415 ++nRange_;
0416 if (preScale_ <= 1)
0417 selectEvent = true;
0418 else if (nRange_ % preScale_ == 1)
0419 selectEvent = true;
0420 }
0421 if (selectEvent) {
0422 for (reco::HcalIsolatedTrackCandidateCollection::const_iterator itr = isotk->begin(); itr != isotk->end();
0423 ++itr)
0424 outputHcalIsoTrackColl->push_back(*itr);
0425
0426 for (reco::VertexCollection::const_iterator vtx = recVtxs->begin(); vtx != recVtxs->end(); ++vtx)
0427 outputVColl->push_back(*vtx);
0428
0429 for (edm::SortedCollection<EcalRecHit>::const_iterator ehit = barrelRecHitsHandle->begin();
0430 ehit != barrelRecHitsHandle->end();
0431 ++ehit)
0432 outputEBColl->push_back(*ehit);
0433
0434 for (edm::SortedCollection<EcalRecHit>::const_iterator ehit = endcapRecHitsHandle->begin();
0435 ehit != endcapRecHitsHandle->end();
0436 ++ehit)
0437 outputEEColl->push_back(*ehit);
0438
0439 for (std::vector<HBHERecHit>::const_iterator hhit = hbhe->begin(); hhit != hbhe->end(); ++hhit)
0440 outputHBHEColl->push_back(*hhit);
0441 ++nGood_;
0442 }
0443 }
0444 }
0445 }
0446 iEvent.put(std::move(outputHcalIsoTrackColl), labelIsoTk_);
0447 iEvent.put(std::move(outputVColl), labelRecVtx_.label());
0448 iEvent.put(std::move(outputEBColl), labelEB_.instance());
0449 iEvent.put(std::move(outputEEColl), labelEE_.instance());
0450 iEvent.put(std::move(outputHBHEColl), labelHBHE_.label());
0451 }
0452
0453 void AlCaIsoTracksProducer::endStream() {
0454 globalCache()->nAll_ += nAll_;
0455 globalCache()->nGood_ += nGood_;
0456 globalCache()->nRange_ += nRange_;
0457 }
0458
0459 void AlCaIsoTracksProducer::globalEndJob(const AlCaIsoTracks::Counters* count) {
0460 edm::LogVerbatim("HcalIsoTrack") << "Finds " << count->nGood_ << " good tracks in " << count->nAll_ << " events and "
0461 << count->nRange_ << " events in the momentum raange";
0462 }
0463
0464 void AlCaIsoTracksProducer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0465 bool changed(false);
0466 edm::LogVerbatim("HcalIsoTrack") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init "
0467 << hltConfig_.init(iRun, iSetup, processName_, changed);
0468
0469 bField = &(iSetup.getData(tok_magField_));
0470 geo = &(iSetup.getData(tok_geom_));
0471 }
0472
0473 void AlCaIsoTracksProducer::endRun(edm::Run const& iRun, edm::EventSetup const&) {
0474 edm::LogVerbatim("HcalIsoTrack") << "endRun [" << nRun_ << "] " << iRun.run();
0475 ++nRun_;
0476 }
0477
0478 reco::HcalIsolatedTrackCandidateCollection* AlCaIsoTracksProducer::select(
0479 edm::Handle<edm::TriggerResults> const& triggerResults,
0480 const std::vector<std::string>& triggerNames_,
0481 edm::Handle<reco::TrackCollection>& trkCollection,
0482 math::XYZPoint& leadPV,
0483 edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
0484 edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
0485 edm::Handle<HBHERecHitCollection>& hbhe,
0486 double ptL1,
0487 double etaL1,
0488 double phiL1) {
0489 reco::HcalIsolatedTrackCandidateCollection* trackCollection = new reco::HcalIsolatedTrackCandidateCollection;
0490 bool ok(false);
0491
0492
0493 for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0494 int hlt = triggerResults->accept(iHLT);
0495 for (unsigned int i = 0; i < trigNames_.size(); ++i) {
0496 if (triggerNames_[iHLT].find(trigNames_[i]) != std::string::npos) {
0497 if (hlt > 0) {
0498 ok = true;
0499 }
0500 edm::LogVerbatim("HcalIsoTrack") << "The trigger we are looking for " << triggerNames_[iHLT] << " Flag " << hlt
0501 << ":" << ok;
0502 }
0503 }
0504 }
0505
0506
0507 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
0508 spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, false);
0509
0510 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
0511 unsigned int nTracks(0), nselTracks(0);
0512 for (trkDetItr = trkCaloDirections.begin(), nTracks = 0; trkDetItr != trkCaloDirections.end();
0513 trkDetItr++, nTracks++) {
0514 const reco::Track* pTrack = &(*(trkDetItr->trkItr));
0515 math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
0516 #ifdef EDM_ML_DEBUG
0517 edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) :" << pTrack->pt() << "|"
0518 << pTrack->eta() << "|" << pTrack->phi() << "|" << pTrack->p();
0519 #endif
0520
0521 bool qltyFlag = spr::goodTrack(pTrack, leadPV, selectionParameter_, false);
0522 #ifdef EDM_ML_DEBUG
0523 edm::LogVerbatim("HcalIsoTrack") << "qltyFlag|okECAL|okHCAL : " << qltyFlag << "|" << trkDetItr->okECAL << "|"
0524 << trkDetItr->okHCAL;
0525 #endif
0526 if (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL) {
0527 double t_p = pTrack->p();
0528 nselTracks++;
0529 int nRH_eMipDR(0), nNearTRKs(0);
0530 double eMipDR = spr::eCone_ecal(geo,
0531 barrelRecHitsHandle,
0532 endcapRecHitsHandle,
0533 trkDetItr->pointHCAL,
0534 trkDetItr->pointECAL,
0535 a_mipR_,
0536 trkDetItr->directionECAL,
0537 nRH_eMipDR);
0538 double hmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, false);
0539 HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
0540 int ieta = detId.ietaAbs();
0541 double eIsolation = (maxRestrictionP_ * exp(slopeRestrictionP_ * ((double)(ieta))));
0542 if (eIsolation < eIsolate_)
0543 eIsolation = eIsolate_;
0544 #ifdef EDM_ML_DEBUG
0545 edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) :" << pTrack->pt() << "|"
0546 << pTrack->eta() << "|" << pTrack->phi() << "|" << t_p << " e_MIP " << eMipDR
0547 << " Chg Isolation " << hmaxNearP << ":" << eIsolation;
0548 #endif
0549 if (t_p > pTrackMin_ && eMipDR < eEcalMax_ && hmaxNearP < eIsolation) {
0550 reco::HcalIsolatedTrackCandidate newCandidate(v4);
0551 newCandidate.SetMaxP(hmaxNearP);
0552 newCandidate.SetEnergyEcal(eMipDR);
0553 newCandidate.setL1(ptL1, etaL1, phiL1);
0554 newCandidate.SetEtaPhiEcal((trkDetItr->pointECAL).eta(), (trkDetItr->pointECAL).phi());
0555 HcalDetId detId = HcalDetId(trkDetItr->detIdHCAL);
0556 newCandidate.SetEtaPhiHcal(
0557 (trkDetItr->pointHCAL).eta(), (trkDetItr->pointHCAL).phi(), detId.ieta(), detId.iphi());
0558 int indx(0);
0559 for (reco::TrackCollection::const_iterator trkItr1 = trkCollection->begin(); trkItr1 != trkCollection->end();
0560 ++trkItr1, ++indx) {
0561 const reco::Track* pTrack1 = &(*trkItr1);
0562 if (pTrack1 == pTrack) {
0563 reco::TrackRef tRef = reco::TrackRef(trkCollection, indx);
0564 newCandidate.setTrack(tRef);
0565 break;
0566 }
0567 }
0568 trackCollection->push_back(newCandidate);
0569 }
0570 }
0571 }
0572 return trackCollection;
0573 }
0574
0575 void AlCaIsoTracksProducer::setPtEtaPhi(std::vector<edm::Ref<l1extra::L1JetParticleCollection> >& objref,
0576 double& ptL1,
0577 double& etaL1,
0578 double& phiL1) {
0579 for (unsigned int p = 0; p < objref.size(); p++) {
0580 if (objref[p]->pt() > ptL1) {
0581 ptL1 = objref[p]->pt();
0582 phiL1 = objref[p]->phi();
0583 etaL1 = objref[p]->eta();
0584 }
0585 }
0586 }
0587
0588 #include "FWCore/Framework/interface/MakerMacros.h"
0589
0590 DEFINE_FWK_MODULE(AlCaIsoTracksProducer);