File indexing completed on 2022-01-18 23:28:32
0001
0002
0003
0004 #include <algorithm>
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
0014
0015 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0016 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0017 #include "Calibration/IsolatedParticles/interface/eCone.h"
0018 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0019 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0020 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0021
0022 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0023 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0024 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0025 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0026 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0027
0028 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0029 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0030 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0031 #include "DataFormats/HcalCalibObjects/interface/HcalIsoTrkCalibVariables.h"
0032 #include "DataFormats/HcalCalibObjects/interface/HcalIsoTrkEventVariables.h"
0033 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0034 #include "DataFormats/Math/interface/LorentzVector.h"
0035 #include "DataFormats/Math/interface/deltaR.h"
0036 #include "DataFormats/MuonReco/interface/Muon.h"
0037 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0038 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0039 #include "DataFormats/TrackReco/interface/Track.h"
0040 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0041 #include "DataFormats/TrackReco/interface/HitPattern.h"
0042 #include "DataFormats/TrackReco/interface/TrackBase.h"
0043 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0044 #include "DataFormats/VertexReco/interface/Vertex.h"
0045 #include "DataFormats/Common/interface/TriggerResults.h"
0046 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0047 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0048 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0049 #include "DataFormats/L1TGlobal/interface/GlobalAlgBlk.h"
0050 #include "DataFormats/L1TGlobal/interface/GlobalExtBlk.h"
0051 #include "DataFormats/L1Trigger/interface/BXVector.h"
0052 #include "L1Trigger/L1TGlobal/interface/L1TGlobalUtil.h"
0053 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0054 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0055 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0056 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0057
0058 #include "FWCore/Common/interface/TriggerNames.h"
0059 #include "FWCore/Framework/interface/Frameworkfwd.h"
0060 #include "FWCore/Framework/interface/stream/EDProducer.h"
0061 #include "FWCore/Framework/interface/Event.h"
0062 #include "FWCore/Framework/interface/EventSetup.h"
0063 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0064 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0065 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0066
0067 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0068 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0069 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0070 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0071 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0072 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0073 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0074 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0075
0076 #include "MagneticField/Engine/interface/MagneticField.h"
0077 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0078 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0079 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0080
0081
0082
0083
0084
0085
0086 namespace alcaHcalIsoTracks {
0087 struct Counters {
0088 Counters() : nAll_(0), nGood_(0), nRange_(0) {}
0089 mutable std::atomic<unsigned int> nAll_, nGood_, nRange_;
0090 };
0091 }
0092
0093 class AlCaHcalIsotrkProducer : public edm::stream::EDProducer<edm::GlobalCache<alcaHcalIsoTracks::Counters>> {
0094 public:
0095 explicit AlCaHcalIsotrkProducer(edm::ParameterSet const&, const alcaHcalIsoTracks::Counters*);
0096 ~AlCaHcalIsotrkProducer() override = default;
0097
0098 static std::unique_ptr<alcaHcalIsoTracks::Counters> initializeGlobalCache(edm::ParameterSet const&) {
0099 return std::make_unique<alcaHcalIsoTracks::Counters>();
0100 }
0101
0102 void produce(edm::Event&, edm::EventSetup const&) override;
0103 void endStream() override;
0104 static void globalEndJob(const alcaHcalIsoTracks::Counters* counters);
0105 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0106
0107 private:
0108 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0109 void endRun(edm::Run const&, edm::EventSetup const&) override;
0110 std::array<int, 3> getProducts(int goodPV,
0111 double eventWeight,
0112 std::vector<math::XYZTLorentzVector>& vecL1,
0113 std::vector<math::XYZTLorentzVector>& vecL3,
0114 math::XYZPoint& leadPV,
0115 std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
0116 std::vector<spr::propagatedTrackID>& trkCaloDets,
0117 const CaloGeometry* geo,
0118 const CaloTopology* topo,
0119 const HcalTopology* theHBHETopology,
0120 const EcalChannelStatus* theEcalChStatus,
0121 const EcalSeverityLevelAlgo* theEcalSevlv,
0122 edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
0123 edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
0124 edm::Handle<HBHERecHitCollection>& hbhe,
0125 edm::Handle<CaloTowerCollection>& towerHandle,
0126 const edm::Handle<reco::GenParticleCollection>& genParticles,
0127 const HcalRespCorrs* respCorrs,
0128 const edm::Handle<reco::MuonCollection>& muonh,
0129 std::vector<HcalIsoTrkCalibVariables>& hocalib,
0130 HcalIsoTrkEventVariables& hocalibEvent,
0131 const edm::EventID& eventId);
0132 double dR(math::XYZTLorentzVector&, math::XYZTLorentzVector&);
0133 double trackP(const reco::Track*, const edm::Handle<reco::GenParticleCollection>&);
0134 double rhoh(const edm::Handle<CaloTowerCollection>&);
0135 double eThreshold(const DetId& id, const CaloGeometry* geo) const;
0136 DetId newId(const DetId&);
0137 void storeEnergy(const HcalRespCorrs* respCorrs,
0138 const std::vector<DetId>& ids,
0139 std::vector<double>& edet,
0140 double& eHcal,
0141 std::vector<unsigned int>& detIds,
0142 std::vector<double>& hitEnergies);
0143 std::pair<double, double> storeEnergy(const HcalRespCorrs* respCorrs,
0144 edm::Handle<HBHERecHitCollection>& hbhe,
0145 const std::vector<DetId>& ids,
0146 std::vector<double>& hitEnergy1,
0147 std::vector<double>& hitEnergy2);
0148 bool notaMuon(const reco::Track* pTrack0, const edm::Handle<reco::MuonCollection>& muonh);
0149
0150
0151 l1t::L1TGlobalUtil* l1GtUtils_;
0152 HLTConfigProvider hltConfig_;
0153 unsigned int nRun_, nAll_, nGood_, nRange_;
0154 const std::vector<std::string> trigNames_;
0155 spr::trackSelectionParameters selectionParameter_;
0156 const std::string theTrackQuality_;
0157 const std::string processName_, l1Filter_;
0158 const std::string l2Filter_, l3Filter_;
0159 const double a_coneR_, a_mipR_, a_mipR2_, a_mipR3_;
0160 const double a_mipR4_, a_mipR5_, pTrackMin_, eEcalMax_;
0161 const double maxRestrictionP_, slopeRestrictionP_;
0162 const double hcalScale_, eIsolate1_, eIsolate2_;
0163 const double pTrackLow_, pTrackHigh_;
0164 const bool ignoreTrigger_, useL1Trigger_;
0165 const bool unCorrect_, collapseDepth_;
0166 const double hitEthrEB_, hitEthrEE0_, hitEthrEE1_;
0167 const double hitEthrEE2_, hitEthrEE3_;
0168 const double hitEthrEELo_, hitEthrEEHi_;
0169 const edm::InputTag triggerEvent_, theTriggerResultsLabel_;
0170 const std::string labelGenTrack_, labelRecVtx_, labelEB_;
0171 const std::string labelEE_, labelHBHE_, labelTower_, l1TrigName_;
0172 const std::vector<int> oldID_, newDepth_;
0173 const bool hep17_;
0174 const std::string labelIsoTkVar_, labelIsoTkEvtVar_;
0175 const std::vector<int> debEvents_;
0176 const bool usePFThresh_;
0177
0178 double a_charIsoR_, a_coneR1_, a_coneR2_;
0179 const HcalDDDRecConstants* hdc_;
0180 const EcalPFRecHitThresholds* eThresholds_;
0181
0182 std::vector<double> etabins_, phibins_;
0183 std::vector<int> oldDet_, oldEta_, oldDepth_;
0184 double etadist_, phidist_, etahalfdist_, phihalfdist_;
0185
0186 edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
0187 edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0188 edm::EDGetTokenT<reco::GenParticleCollection> tok_parts_;
0189 edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0190 edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0191 edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0192 edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0193 edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0194 edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0195 edm::EDGetTokenT<CaloTowerCollection> tok_cala_;
0196 edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
0197 edm::EDGetTokenT<BXVector<GlobalAlgBlk>> tok_alg_;
0198 edm::EDGetTokenT<reco::MuonCollection> tok_Muon_;
0199
0200 edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec_;
0201 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_bFieldH_;
0202 edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_ecalChStatus_;
0203 edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0204 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0205 edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_caloTopology_;
0206 edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
0207 edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_resp_;
0208 edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> tok_ecalPFRecHitThresholds_;
0209
0210 bool debug_;
0211 };
0212
0213 AlCaHcalIsotrkProducer::AlCaHcalIsotrkProducer(edm::ParameterSet const& iConfig,
0214 const alcaHcalIsoTracks::Counters* counters)
0215 : nRun_(0),
0216 nAll_(0),
0217 nGood_(0),
0218 nRange_(0),
0219 trigNames_(iConfig.getParameter<std::vector<std::string>>("triggers")),
0220 theTrackQuality_(iConfig.getParameter<std::string>("trackQuality")),
0221 processName_(iConfig.getParameter<std::string>("processName")),
0222 l1Filter_(iConfig.getParameter<std::string>("l1Filter")),
0223 l2Filter_(iConfig.getParameter<std::string>("l2Filter")),
0224 l3Filter_(iConfig.getParameter<std::string>("l3Filter")),
0225 a_coneR_(iConfig.getParameter<double>("coneRadius")),
0226 a_mipR_(iConfig.getParameter<double>("coneRadiusMIP")),
0227 a_mipR2_(iConfig.getParameter<double>("coneRadiusMIP2")),
0228 a_mipR3_(iConfig.getParameter<double>("coneRadiusMIP3")),
0229 a_mipR4_(iConfig.getParameter<double>("coneRadiusMIP4")),
0230 a_mipR5_(iConfig.getParameter<double>("coneRadiusMIP5")),
0231 pTrackMin_(iConfig.getParameter<double>("minimumTrackP")),
0232 eEcalMax_(iConfig.getParameter<double>("maximumEcalEnergy")),
0233 maxRestrictionP_(iConfig.getParameter<double>("maxTrackP")),
0234 slopeRestrictionP_(iConfig.getParameter<double>("slopeTrackP")),
0235 hcalScale_(iConfig.getUntrackedParameter<double>("hHcalScale", 1.0)),
0236 eIsolate1_(iConfig.getParameter<double>("isolationEnergyTight")),
0237 eIsolate2_(iConfig.getParameter<double>("isolationEnergyLoose")),
0238 pTrackLow_(iConfig.getParameter<double>("momentumLow")),
0239 pTrackHigh_(iConfig.getParameter<double>("momentumHigh")),
0240 ignoreTrigger_(iConfig.getUntrackedParameter<bool>("ignoreTriggers", false)),
0241 useL1Trigger_(iConfig.getUntrackedParameter<bool>("useL1Trigger", false)),
0242 unCorrect_(iConfig.getUntrackedParameter<bool>("unCorrect", false)),
0243 collapseDepth_(iConfig.getUntrackedParameter<bool>("collapseDepth", false)),
0244 hitEthrEB_(iConfig.getParameter<double>("EBHitEnergyThreshold")),
0245 hitEthrEE0_(iConfig.getParameter<double>("EEHitEnergyThreshold0")),
0246 hitEthrEE1_(iConfig.getParameter<double>("EEHitEnergyThreshold1")),
0247 hitEthrEE2_(iConfig.getParameter<double>("EEHitEnergyThreshold2")),
0248 hitEthrEE3_(iConfig.getParameter<double>("EEHitEnergyThreshold3")),
0249 hitEthrEELo_(iConfig.getParameter<double>("EEHitEnergyThresholdLow")),
0250 hitEthrEEHi_(iConfig.getParameter<double>("EEHitEnergyThresholdHigh")),
0251 triggerEvent_(iConfig.getParameter<edm::InputTag>("labelTriggerEvent")),
0252 theTriggerResultsLabel_(iConfig.getParameter<edm::InputTag>("labelTriggerResult")),
0253 labelGenTrack_(iConfig.getParameter<std::string>("labelTrack")),
0254 labelRecVtx_(iConfig.getParameter<std::string>("labelVertex")),
0255 labelEB_(iConfig.getParameter<std::string>("labelEBRecHit")),
0256 labelEE_(iConfig.getParameter<std::string>("labelEERecHit")),
0257 labelHBHE_(iConfig.getParameter<std::string>("labelHBHERecHit")),
0258 labelTower_(iConfig.getParameter<std::string>("labelCaloTower")),
0259 l1TrigName_(iConfig.getUntrackedParameter<std::string>("l1TrigName", "L1_SingleJet60")),
0260 oldID_(iConfig.getUntrackedParameter<std::vector<int>>("oldID")),
0261 newDepth_(iConfig.getUntrackedParameter<std::vector<int>>("newDepth")),
0262 hep17_(iConfig.getUntrackedParameter<bool>("hep17")),
0263 labelIsoTkVar_(iConfig.getParameter<std::string>("isoTrackLabel")),
0264 labelIsoTkEvtVar_(iConfig.getParameter<std::string>("isoTrackEventLabel")),
0265 debEvents_(iConfig.getParameter<std::vector<int>>("debugEvents")),
0266 usePFThresh_(iConfig.getParameter<bool>("usePFThreshold")) {
0267
0268 const double isolationRadius(28.9), innerR(10.0), outerR(30.0);
0269 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality_);
0270 selectionParameter_.minPt = iConfig.getParameter<double>("minTrackPt");
0271 selectionParameter_.maxChi2 = iConfig.getParameter<double>("maxChi2");
0272 selectionParameter_.maxDpOverP = iConfig.getParameter<double>("maxDpOverP");
0273 selectionParameter_.minOuterHit = iConfig.getParameter<int>("minOuterHit");
0274 selectionParameter_.minLayerCrossed = iConfig.getParameter<int>("minLayerCrossed");
0275 selectionParameter_.maxInMiss = iConfig.getParameter<int>("maxInMiss");
0276 selectionParameter_.minQuality = trackQuality_;
0277 selectionParameter_.maxDxyPV = iConfig.getParameter<double>("maxDxyPV");
0278 selectionParameter_.maxDzPV = iConfig.getParameter<double>("maxDzPV");
0279 selectionParameter_.maxChi2 = iConfig.getParameter<double>("maxChi2");
0280 selectionParameter_.maxDpOverP = iConfig.getParameter<double>("maxDpOverP");
0281 selectionParameter_.minOuterHit = iConfig.getParameter<int>("minOuterHit");
0282 selectionParameter_.minLayerCrossed = iConfig.getParameter<int>("minLayerCrossed");
0283 selectionParameter_.maxInMiss = iConfig.getParameter<int>("maxInMiss");
0284 selectionParameter_.maxOutMiss = iConfig.getParameter<int>("maxOutMiss");
0285 a_charIsoR_ = a_coneR_ + isolationRadius;
0286 a_coneR1_ = a_coneR_ + innerR;
0287 a_coneR2_ = a_coneR_ + outerR;
0288
0289
0290
0291
0292
0293 std::string labelBS = iConfig.getParameter<std::string>("labelBeamSpot");
0294 edm::InputTag algTag = iConfig.getParameter<edm::InputTag>("algInputTag");
0295 edm::InputTag extTag = iConfig.getParameter<edm::InputTag>("extInputTag");
0296 std::string labelMuon = iConfig.getParameter<std::string>("labelMuon");
0297
0298 for (unsigned int k = 0; k < oldID_.size(); ++k) {
0299 oldDet_.emplace_back((oldID_[k] / 10000) % 10);
0300 oldEta_.emplace_back((oldID_[k] / 100) % 100);
0301 oldDepth_.emplace_back(oldID_[k] % 100);
0302 }
0303 l1GtUtils_ = new l1t::L1TGlobalUtil(iConfig, consumesCollector(), *this, algTag, extTag, l1t::UseEventSetupIn::Event);
0304
0305 tok_trigEvt_ = consumes<trigger::TriggerEvent>(triggerEvent_);
0306 tok_trigRes_ = consumes<edm::TriggerResults>(theTriggerResultsLabel_);
0307 tok_bs_ = consumes<reco::BeamSpot>(labelBS);
0308 tok_genTrack_ = consumes<reco::TrackCollection>(labelGenTrack_);
0309 tok_ew_ = consumes<GenEventInfoProduct>(edm::InputTag("generator"));
0310 tok_parts_ = consumes<reco::GenParticleCollection>(edm::InputTag("genParticles"));
0311 tok_cala_ = consumes<CaloTowerCollection>(labelTower_);
0312 tok_alg_ = consumes<BXVector<GlobalAlgBlk>>(algTag);
0313 tok_Muon_ = consumes<reco::MuonCollection>(labelMuon);
0314 tok_recVtx_ = consumes<reco::VertexCollection>(labelRecVtx_);
0315 tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", labelEB_));
0316 tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", labelEE_));
0317 tok_hbhe_ = consumes<HBHERecHitCollection>(labelHBHE_);
0318 edm::LogVerbatim("HcalIsoTrack") << "Labels used " << triggerEvent_ << " " << theTriggerResultsLabel_ << " "
0319 << labelBS << " " << labelRecVtx_ << " " << labelGenTrack_ << " "
0320 << edm::InputTag("ecalRecHit", labelEB_) << " "
0321 << edm::InputTag("ecalRecHit", labelEE_) << " " << labelHBHE_ << " " << labelTower_
0322 << " " << labelMuon;
0323
0324 tok_ddrec_ = esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>();
0325 tok_bFieldH_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0326 tok_ecalChStatus_ = esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0327 tok_sevlv_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0328 tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0329 tok_caloTopology_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0330 tok_htopo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0331 tok_resp_ = esConsumes<HcalRespCorrs, HcalRespCorrsRcd>();
0332 tok_ecalPFRecHitThresholds_ = esConsumes<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd>();
0333
0334 edm::LogVerbatim("HcalIsoTrack")
0335 << "Parameters read from config file \n"
0336 << "\t minPt " << selectionParameter_.minPt << "\t theTrackQuality " << theTrackQuality_ << "\t minQuality "
0337 << selectionParameter_.minQuality << "\t maxDxyPV " << selectionParameter_.maxDxyPV << "\t maxDzPV "
0338 << selectionParameter_.maxDzPV << "\t maxChi2 " << selectionParameter_.maxChi2 << "\t maxDpOverP "
0339 << selectionParameter_.maxDpOverP << "\t minOuterHit " << selectionParameter_.minOuterHit << "\t minLayerCrossed "
0340 << selectionParameter_.minLayerCrossed << "\t maxInMiss " << selectionParameter_.maxInMiss << "\t maxOutMiss "
0341 << selectionParameter_.maxOutMiss << "\t a_coneR " << a_coneR_ << ":" << a_coneR1_ << ":" << a_coneR2_
0342 << "\t a_charIsoR " << a_charIsoR_ << "\t a_mipR " << a_mipR_ << "\t a_mipR2 " << a_mipR2_ << "\t a_mipR3 "
0343 << a_mipR3_ << "\t a_mipR4 " << a_mipR4_ << "\t a_mipR5 " << a_mipR5_ << "\n pTrackMin_ " << pTrackMin_
0344 << "\t eEcalMax_ " << eEcalMax_ << "\t maxRestrictionP_ " << maxRestrictionP_ << "\t slopeRestrictionP_ "
0345 << slopeRestrictionP_ << "\t eIsolateStrong_ " << eIsolate1_ << "\t eIsolateSoft_ " << eIsolate2_
0346 << "\t hcalScale_ " << hcalScale_ << "\n\t momentumLow_ " << pTrackLow_ << "\t momentumHigh_ " << pTrackHigh_
0347 << "\n\t ignoreTrigger_ " << ignoreTrigger_ << "\n\t useL1Trigger_ " << useL1Trigger_ << "\t unCorrect_ "
0348 << unCorrect_ << "\t collapseDepth_ " << collapseDepth_ << "\t L1TrigName_ " << l1TrigName_
0349 << "\nThreshold flag used " << usePFThresh_ << " value for EB " << hitEthrEB_ << " EE " << hitEthrEE0_ << ":"
0350 << hitEthrEE1_ << ":" << hitEthrEE2_ << ":" << hitEthrEE3_ << ":" << hitEthrEELo_ << ":" << hitEthrEEHi_;
0351 edm::LogVerbatim("HcalIsoTrack") << "Process " << processName_ << " L1Filter:" << l1Filter_
0352 << " L2Filter:" << l2Filter_ << " L3Filter:" << l3Filter_ << " and "
0353 << debEvents_.size() << " events to be debugged";
0354 for (unsigned int k = 0; k < trigNames_.size(); ++k) {
0355 edm::LogVerbatim("HcalIsoTrack") << "Trigger[" << k << "] " << trigNames_[k];
0356 }
0357 edm::LogVerbatim("HcalIsoTrack") << oldID_.size() << " DetIDs to be corrected with HEP17 flag:" << hep17_;
0358 for (unsigned int k = 0; k < oldID_.size(); ++k)
0359 edm::LogVerbatim("HcalIsoTrack") << "[" << k << "] Det " << oldDet_[k] << " EtaAbs " << oldEta_[k] << " Depth "
0360 << oldDepth_[k] << ":" << newDepth_[k];
0361
0362 for (int i = 0; i < 10; i++)
0363 phibins_.push_back(-M_PI + 0.1 * (2 * i + 1) * M_PI);
0364 for (int i = 0; i < 8; ++i)
0365 etabins_.push_back(-2.1 + 0.6 * i);
0366 etadist_ = etabins_[1] - etabins_[0];
0367 phidist_ = phibins_[1] - phibins_[0];
0368 etahalfdist_ = 0.5 * etadist_;
0369 phihalfdist_ = 0.5 * phidist_;
0370 edm::LogVerbatim("HcalIsoTrack") << "EtaDist " << etadist_ << " " << etahalfdist_ << " PhiDist " << phidist_ << " "
0371 << phihalfdist_;
0372
0373 produces<HcalIsoTrkCalibVariablesCollection>(labelIsoTkVar_);
0374 produces<HcalIsoTrkEventVariablesCollection>(labelIsoTkEvtVar_);
0375 edm::LogVerbatim("HcalIsoTrack") << " Expected to produce the collections:\n"
0376 << "HcalIsoTrkCalibVariablesCollection with label " << labelIsoTkVar_
0377 << "\nand HcalIsoTrkEventVariablesCollection with label " << labelIsoTkEvtVar_;
0378 }
0379
0380 void AlCaHcalIsotrkProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0381 edm::ParameterSetDescription desc;
0382 desc.add<std::string>("processName", "HLT");
0383 std::vector<std::string> trig;
0384 desc.add<std::vector<std::string>>("triggers", trig);
0385 desc.add<std::string>("l1Filter", "");
0386 desc.add<std::string>("l2Filter", "L2Filter");
0387 desc.add<std::string>("l3Filter", "Filter");
0388
0389 desc.add<std::string>("trackQuality", "highPurity");
0390 desc.add<double>("minTrackPt", 1.0);
0391 desc.add<double>("maxDxyPV", 0.02);
0392 desc.add<double>("maxDzPV", 0.02);
0393 desc.add<double>("maxChi2", 5.0);
0394 desc.add<double>("maxDpOverP", 0.1);
0395 desc.add<int>("minOuterHit", 4);
0396 desc.add<int>("minLayerCrossed", 8);
0397 desc.add<int>("maxInMiss", 0);
0398 desc.add<int>("maxOutMiss", 0);
0399
0400 desc.add<double>("minimumTrackP", 10.0);
0401 desc.add<double>("coneRadius", 34.98);
0402
0403 desc.add<double>("coneRadiusMIP", 14.0);
0404 desc.add<double>("coneRadiusMIP2", 18.0);
0405 desc.add<double>("coneRadiusMIP3", 20.0);
0406 desc.add<double>("coneRadiusMIP4", 22.0);
0407 desc.add<double>("coneRadiusMIP5", 24.0);
0408 desc.add<double>("maximumEcalEnergy", 2.0);
0409
0410 desc.add<double>("maxTrackP", 8.0);
0411 desc.add<double>("slopeTrackP", 0.05090504066);
0412 desc.add<double>("isolationEnergyTight", 2.0);
0413 desc.add<double>("isolationEnergyLoose", 10.0);
0414
0415 desc.add<double>("EBHitEnergyThreshold", 0.08);
0416 desc.add<double>("EEHitEnergyThreshold0", 0.30);
0417 desc.add<double>("EEHitEnergyThreshold1", 0.00);
0418 desc.add<double>("EEHitEnergyThreshold2", 0.00);
0419 desc.add<double>("EEHitEnergyThreshold3", 0.00);
0420 desc.add<double>("EEHitEnergyThresholdLow", 0.30);
0421 desc.add<double>("EEHitEnergyThresholdHigh", 0.30);
0422
0423 desc.add<double>("momentumLow", 40.0);
0424 desc.add<double>("momentumHigh", 60.0);
0425 desc.add<int>("prescaleLow", 1);
0426 desc.add<int>("prescaleHigh", 1);
0427
0428 desc.add<edm::InputTag>("labelTriggerEvent", edm::InputTag("hltTriggerSummaryAOD", "", "HLT"));
0429 desc.add<edm::InputTag>("labelTriggerResult", edm::InputTag("TriggerResults", "", "HLT"));
0430 desc.add<std::string>("labelTrack", "generalTracks");
0431 desc.add<std::string>("labelVertex", "offlinePrimaryVertices");
0432 desc.add<std::string>("labelEBRecHit", "EcalRecHitsEB");
0433 desc.add<std::string>("labelEERecHit", "EcalRecHitsEE");
0434 desc.add<std::string>("labelHBHERecHit", "hbhereco");
0435 desc.add<std::string>("labelBeamSpot", "offlineBeamSpot");
0436 desc.add<std::string>("labelCaloTower", "towerMaker");
0437 desc.add<std::string>("labelMuon", "muons");
0438 desc.add<edm::InputTag>("algInputTag", edm::InputTag("gtStage2Digis"));
0439 desc.add<edm::InputTag>("extInputTag", edm::InputTag("gtStage2Digis"));
0440 desc.add<std::string>("isoTrackLabel", "HcalIsoTrack");
0441 desc.add<std::string>("isoTrackEventLabel", "HcalIsoTrackEvent");
0442
0443
0444 desc.addUntracked<bool>("ignoreTriggers", false);
0445 desc.addUntracked<bool>("useL1Trigger", false);
0446 desc.addUntracked<double>("hcalScale", 1.0);
0447 desc.addUntracked<bool>("unCorrect", false);
0448 desc.addUntracked<bool>("collapseDepth", false);
0449 desc.addUntracked<std::string>("l1TrigName", "L1_SingleJet60");
0450 std::vector<int> dummy;
0451 desc.addUntracked<std::vector<int>>("oldID", dummy);
0452 desc.addUntracked<std::vector<int>>("newDepth", dummy);
0453 desc.addUntracked<bool>("hep17", false);
0454 std::vector<int> events;
0455 desc.add<std::vector<int>>("debugEvents", events);
0456 desc.add<bool>("usePFThreshold", true);
0457 descriptions.add("alcaHcalIsotrkProducer", desc);
0458 }
0459
0460 void AlCaHcalIsotrkProducer::produce(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0461 nAll_++;
0462 debug_ = (debEvents_.empty())
0463 ? true
0464 : (std::find(debEvents_.begin(), debEvents_.end(), iEvent.id().event()) != debEvents_.end());
0465 #ifdef EDM_ML_DEBUG
0466 if (debug_)
0467 edm::LogVerbatim("HcalIsoTrack") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event()
0468 << " Luminosity " << iEvent.luminosityBlock() << " Bunch "
0469 << iEvent.bunchCrossing();
0470 #endif
0471
0472 HcalIsoTrkEventVariables isoTrkEvent;
0473
0474 const MagneticField* bField = &iSetup.getData(tok_bFieldH_);
0475 const EcalChannelStatus* theEcalChStatus = &iSetup.getData(tok_ecalChStatus_);
0476 const EcalSeverityLevelAlgo* theEcalSevlv = &iSetup.getData(tok_sevlv_);
0477 eThresholds_ = &iSetup.getData(tok_ecalPFRecHitThresholds_);
0478
0479
0480 const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0481 const CaloTopology* caloTopology = &iSetup.getData(tok_caloTopology_);
0482 const HcalTopology* theHBHETopology = &iSetup.getData(tok_htopo_);
0483
0484
0485 const HcalRespCorrs* respCorrs = &iSetup.getData(tok_resp_);
0486
0487 bool okC(true);
0488
0489 auto trkCollection = iEvent.getHandle(tok_genTrack_);
0490 if (!trkCollection.isValid()) {
0491 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
0492 okC = false;
0493 }
0494
0495
0496 auto const& muonh = iEvent.getHandle(tok_Muon_);
0497
0498
0499 auto const& recVtxs = iEvent.getHandle(tok_recVtx_);
0500 auto const& beamSpotH = iEvent.getHandle(tok_bs_);
0501 math::XYZPoint leadPV(0, 0, 0);
0502 int goodPV(0);
0503 if (recVtxs.isValid() && !(recVtxs->empty())) {
0504 isoTrkEvent.allvertex_ = recVtxs->size();
0505 for (unsigned int k = 0; k < recVtxs->size(); ++k) {
0506 if (!((*recVtxs)[k].isFake()) && ((*recVtxs)[k].ndof() > 4)) {
0507 if (goodPV == 0)
0508 leadPV = math::XYZPoint((*recVtxs)[k].x(), (*recVtxs)[k].y(), (*recVtxs)[k].z());
0509 goodPV++;
0510 }
0511 }
0512 }
0513 if (goodPV == 0 && beamSpotH.isValid()) {
0514 leadPV = beamSpotH->position();
0515 }
0516 #ifdef EDM_ML_DEBUG
0517 if (debug_) {
0518 edm::LogVerbatim("HcalIsoTrack") << "Primary Vertex " << leadPV << " out of " << goodPV << " vertex";
0519 if (beamSpotH.isValid())
0520 edm::LogVerbatim("HcalIsoTrack") << " Beam Spot " << beamSpotH->position();
0521 }
0522 #endif
0523
0524
0525 auto barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0526 if (!barrelRecHitsHandle.isValid()) {
0527 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEB_;
0528 okC = false;
0529 }
0530 auto endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0531 if (!endcapRecHitsHandle.isValid()) {
0532 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEE_;
0533 okC = false;
0534 }
0535 auto hbhe = iEvent.getHandle(tok_hbhe_);
0536 if (!hbhe.isValid()) {
0537 edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelHBHE_;
0538 okC = false;
0539 }
0540 auto caloTower = iEvent.getHandle(tok_cala_);
0541
0542
0543 auto const& genParticles = iEvent.getHandle(tok_parts_);
0544 auto const& genEventInfo = iEvent.getHandle(tok_ew_);
0545 double eventWeight = (genEventInfo.isValid()) ? genEventInfo->weight() : 1.0;
0546
0547
0548 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
0549 spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, false);
0550 std::vector<spr::propagatedTrackID> trkCaloDets;
0551 spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDets, false);
0552 std::vector<math::XYZTLorentzVector> vecL1, vecL3;
0553 isoTrkEvent.tracks_ = trkCollection->size();
0554 isoTrkEvent.tracksProp_ = trkCaloDirections.size();
0555 isoTrkEvent.hltbits_.assign(trigNames_.size(), false);
0556
0557 if (!ignoreTrigger_) {
0558
0559 l1GtUtils_->retrieveL1(iEvent, iSetup, tok_alg_);
0560 const std::vector<std::pair<std::string, bool>>& finalDecisions = l1GtUtils_->decisionsFinal();
0561 for (const auto& decision : finalDecisions) {
0562 if (decision.first.find(l1TrigName_) != std::string::npos) {
0563 isoTrkEvent.l1Bit_ = decision.second;
0564 break;
0565 }
0566 }
0567 #ifdef EDM_ML_DEBUG
0568 if (debug_)
0569 edm::LogVerbatim("HcalIsoTrack") << "Trigger Information for " << l1TrigName_ << " is " << isoTrkEvent.l1Bit_
0570 << " from a list of " << finalDecisions.size() << " decisions";
0571 #endif
0572
0573
0574 auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
0575 if (triggerResults.isValid()) {
0576 const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0577 const std::vector<std::string>& names = triggerNames.triggerNames();
0578 if (!trigNames_.empty()) {
0579 for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0580 int hlt = triggerResults->accept(iHLT);
0581 for (unsigned int i = 0; i < trigNames_.size(); ++i) {
0582 if (names[iHLT].find(trigNames_[i]) != std::string::npos) {
0583 isoTrkEvent.hltbits_[i] = (hlt > 0);
0584 if (hlt > 0)
0585 isoTrkEvent.trigPass_ = true;
0586 #ifdef EDM_ML_DEBUG
0587 if (debug_)
0588 edm::LogVerbatim("HcalIsoTrack")
0589 << "This trigger " << names[iHLT] << " Flag " << hlt << ":" << isoTrkEvent.hltbits_[i];
0590 #endif
0591 }
0592 }
0593 }
0594 }
0595 }
0596 #ifdef EDM_ML_DEBUG
0597 if (debug_)
0598 edm::LogVerbatim("HcalIsoTrack") << "HLT Information shows " << isoTrkEvent.trigPass_ << ":" << trigNames_.empty()
0599 << ":" << okC;
0600 #endif
0601 }
0602
0603 auto outputHcalIsoTrackColl = std::make_unique<HcalIsoTrkCalibVariablesCollection>();
0604 std::array<int, 3> ntksave{{0, 0, 0}};
0605 if (ignoreTrigger_ || useL1Trigger_) {
0606 if (ignoreTrigger_ || isoTrkEvent.l1Bit_)
0607 ntksave = getProducts(goodPV,
0608 eventWeight,
0609 vecL1,
0610 vecL3,
0611 leadPV,
0612 trkCaloDirections,
0613 trkCaloDets,
0614 geo,
0615 caloTopology,
0616 theHBHETopology,
0617 theEcalChStatus,
0618 theEcalSevlv,
0619 barrelRecHitsHandle,
0620 endcapRecHitsHandle,
0621 hbhe,
0622 caloTower,
0623 genParticles,
0624 respCorrs,
0625 muonh,
0626 *outputHcalIsoTrackColl,
0627 isoTrkEvent,
0628 iEvent.id());
0629 isoTrkEvent.tracksSaved_ = ntksave[0];
0630 isoTrkEvent.tracksLoose_ = ntksave[1];
0631 isoTrkEvent.tracksTight_ = ntksave[2];
0632 } else {
0633 trigger::TriggerEvent triggerEvent;
0634 auto const& triggerEventHandle = iEvent.getHandle(tok_trigEvt_);
0635 if (!triggerEventHandle.isValid()) {
0636 edm::LogWarning("HcalIsoTrack") << "Error! Can't get the product " << triggerEvent_.label();
0637 } else if (okC) {
0638 triggerEvent = *(triggerEventHandle.product());
0639 const trigger::TriggerObjectCollection& TOC(triggerEvent.getObjects());
0640 bool done(false);
0641 auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
0642 if (triggerResults.isValid()) {
0643 std::vector<std::string> modules;
0644 const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0645 const std::vector<std::string>& names = triggerNames.triggerNames();
0646 for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0647 bool ok = (isoTrkEvent.trigPass_) || (trigNames_.empty());
0648 if (ok) {
0649 unsigned int triggerindx = hltConfig_.triggerIndex(names[iHLT]);
0650 const std::vector<std::string>& moduleLabels(hltConfig_.moduleLabels(triggerindx));
0651 std::vector<math::XYZTLorentzVector> vecL2;
0652 vecL1.clear();
0653 vecL3.clear();
0654
0655 for (unsigned int ifilter = 0; ifilter < triggerEvent.sizeFilters(); ++ifilter) {
0656 std::vector<int> Keys;
0657 std::string label = triggerEvent.filterTag(ifilter).label();
0658
0659 for (unsigned int imodule = 0; imodule < moduleLabels.size(); imodule++) {
0660 if (label.find(moduleLabels[imodule]) != std::string::npos) {
0661 #ifdef EDM_ML_DEBUG
0662 if (debug_)
0663 edm::LogVerbatim("HcalIsoTrack") << "FilterName " << label;
0664 #endif
0665 for (unsigned int ifiltrKey = 0; ifiltrKey < triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
0666 Keys.push_back(triggerEvent.filterKeys(ifilter)[ifiltrKey]);
0667 const trigger::TriggerObject& TO(TOC[Keys[ifiltrKey]]);
0668 math::XYZTLorentzVector v4(TO.px(), TO.py(), TO.pz(), TO.energy());
0669 if (label.find(l2Filter_) != std::string::npos) {
0670 vecL2.push_back(v4);
0671 } else if (label.find(l3Filter_) != std::string::npos) {
0672 vecL3.push_back(v4);
0673 } else if ((label.find(l1Filter_) != std::string::npos) || (l1Filter_.empty())) {
0674 vecL1.push_back(v4);
0675 }
0676 #ifdef EDM_ML_DEBUG
0677 if (debug_)
0678 edm::LogVerbatim("HcalIsoTrack")
0679 << "key " << ifiltrKey << " : pt " << TO.pt() << " eta " << TO.eta() << " phi " << TO.phi()
0680 << " mass " << TO.mass() << " Id " << TO.id();
0681 #endif
0682 }
0683 #ifdef EDM_ML_DEBUG
0684 if (debug_)
0685 edm::LogVerbatim("HcalIsoTrack")
0686 << "sizes " << vecL1.size() << ":" << vecL2.size() << ":" << vecL3.size();
0687 #endif
0688 }
0689 }
0690 }
0691
0692
0693 if (!done) {
0694 ntksave = getProducts(goodPV,
0695 eventWeight,
0696 vecL1,
0697 vecL3,
0698 leadPV,
0699 trkCaloDirections,
0700 trkCaloDets,
0701 geo,
0702 caloTopology,
0703 theHBHETopology,
0704 theEcalChStatus,
0705 theEcalSevlv,
0706 barrelRecHitsHandle,
0707 endcapRecHitsHandle,
0708 hbhe,
0709 caloTower,
0710 genParticles,
0711 respCorrs,
0712 muonh,
0713 *outputHcalIsoTrackColl,
0714 isoTrkEvent,
0715 iEvent.id());
0716 isoTrkEvent.tracksSaved_ += ntksave[0];
0717 isoTrkEvent.tracksLoose_ += ntksave[1];
0718 isoTrkEvent.tracksTight_ += ntksave[2];
0719 done = true;
0720 }
0721 }
0722 }
0723 }
0724 }
0725 }
0726 isoTrkEvent.trigPassSel_ = (isoTrkEvent.tracksSaved_ > 0);
0727 if (isoTrkEvent.trigPassSel_) {
0728 ++nGood_;
0729 for (auto itr = outputHcalIsoTrackColl->begin(); itr != outputHcalIsoTrackColl->end(); ++itr) {
0730 if (itr->p_ > pTrackLow_ && itr->p_ < pTrackHigh_)
0731 ++nRange_;
0732 }
0733 }
0734
0735 auto outputEventcol = std::make_unique<HcalIsoTrkEventVariablesCollection>();
0736 outputEventcol->emplace_back(isoTrkEvent);
0737 iEvent.put(std::move(outputEventcol), labelIsoTkEvtVar_);
0738 iEvent.put(std::move(outputHcalIsoTrackColl), labelIsoTkVar_);
0739 }
0740
0741 void AlCaHcalIsotrkProducer::endStream() {
0742 globalCache()->nAll_ += nAll_;
0743 globalCache()->nGood_ += nGood_;
0744 globalCache()->nRange_ += nRange_;
0745 }
0746
0747 void AlCaHcalIsotrkProducer::globalEndJob(const alcaHcalIsoTracks::Counters* count) {
0748 edm::LogVerbatim("HcalIsoTrack") << "Finds " << count->nGood_ << " good tracks in " << count->nAll_ << " events and "
0749 << count->nRange_ << " events in the momentum range";
0750 }
0751
0752 void AlCaHcalIsotrkProducer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0753 hdc_ = &iSetup.getData(tok_ddrec_);
0754 if (!ignoreTrigger_) {
0755 bool changed(false);
0756 edm::LogVerbatim("HcalIsoTrack") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init "
0757 << hltConfig_.init(iRun, iSetup, processName_, changed);
0758 }
0759 }
0760
0761 void AlCaHcalIsotrkProducer::endRun(edm::Run const& iRun, edm::EventSetup const&) {
0762 edm::LogVerbatim("HcalIsoTrack") << "endRun [" << nRun_ << "] " << iRun.run();
0763 ++nRun_;
0764 }
0765
0766 std::array<int, 3> AlCaHcalIsotrkProducer::getProducts(int goodPV,
0767 double eventWeight,
0768 std::vector<math::XYZTLorentzVector>& vecL1,
0769 std::vector<math::XYZTLorentzVector>& vecL3,
0770 math::XYZPoint& leadPV,
0771 std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
0772 std::vector<spr::propagatedTrackID>& trkCaloDets,
0773 const CaloGeometry* geo,
0774 const CaloTopology* caloTopology,
0775 const HcalTopology* theHBHETopology,
0776 const EcalChannelStatus* theEcalChStatus,
0777 const EcalSeverityLevelAlgo* theEcalSevlv,
0778 edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
0779 edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
0780 edm::Handle<HBHERecHitCollection>& hbhe,
0781 edm::Handle<CaloTowerCollection>& tower,
0782 const edm::Handle<reco::GenParticleCollection>& genParticles,
0783 const HcalRespCorrs* respCorrs,
0784 const edm::Handle<reco::MuonCollection>& muonh,
0785 std::vector<HcalIsoTrkCalibVariables>& hocalib,
0786 HcalIsoTrkEventVariables& hocalibEvent,
0787 const edm::EventID& eventId) {
0788 int nSave(0), nLoose(0), nTight(0);
0789 unsigned int nTracks(0), nselTracks(0);
0790 double rhohEV = (tower.isValid()) ? rhoh(tower) : 0;
0791
0792
0793 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
0794 for (trkDetItr = trkCaloDirections.begin(), nTracks = 0; trkDetItr != trkCaloDirections.end();
0795 trkDetItr++, nTracks++) {
0796 const reco::Track* pTrack = &(*(trkDetItr->trkItr));
0797 math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
0798 bool accept(false);
0799 HcalIsoTrkCalibVariables isoTk;
0800 isoTk.eventWeight_ = eventWeight;
0801 isoTk.goodPV_ = goodPV;
0802 isoTk.nVtx_ = hocalibEvent.allvertex_;
0803 isoTk.nTrk_ = trkCaloDirections.size();
0804 isoTk.rhoh_ = rhohEV;
0805 for (const auto& trig : hocalibEvent.hltbits_)
0806 isoTk.trgbits_.emplace_back(trig);
0807 if (!vecL1.empty()) {
0808 isoTk.l1pt_ = vecL1[0].pt();
0809 isoTk.l1eta_ = vecL1[0].eta();
0810 isoTk.l1phi_ = vecL1[0].phi();
0811 }
0812 if (!vecL3.empty()) {
0813 isoTk.l3pt_ = vecL3[0].pt();
0814 isoTk.l3eta_ = vecL3[0].eta();
0815 isoTk.l3phi_ = vecL3[0].phi();
0816 }
0817
0818 isoTk.p_ = pTrack->p();
0819 isoTk.pt_ = pTrack->pt();
0820 isoTk.phi_ = pTrack->phi();
0821 #ifdef EDM_ML_DEBUG
0822 if (debug_)
0823 edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) : " << isoTk.pt_ << "|"
0824 << pTrack->eta() << "|" << isoTk.phi_ << "|" << isoTk.p_;
0825 int flag(0);
0826 #endif
0827 isoTk.mindR2_ = 999;
0828 for (unsigned int k = 0; k < vecL3.size(); ++k) {
0829 double dr = dR(vecL3[k], v4);
0830 if (dr < isoTk.mindR2_) {
0831 isoTk.mindR2_ = dr;
0832 }
0833 }
0834 isoTk.mindR1_ = (!vecL1.empty()) ? dR(vecL1[0], v4) : 999;
0835 #ifdef EDM_ML_DEBUG
0836 if (debug_)
0837 edm::LogVerbatim("HcalIsoTrack") << "Closest L3 object at dr : " << isoTk.mindR2_ << " and from L1 "
0838 << isoTk.mindR1_;
0839 #endif
0840 isoTk.ieta_ = isoTk.iphi_ = 0;
0841 if (trkDetItr->okHCAL) {
0842 HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
0843 isoTk.ieta_ = detId.ieta();
0844 isoTk.iphi_ = detId.iphi();
0845 if (isoTk.p_ > 40.0 && isoTk.p_ <= 60.0)
0846 hocalibEvent.ietaAll_.emplace_back(isoTk.ieta_);
0847 }
0848
0849 isoTk.selectTk_ = spr::goodTrack(pTrack, leadPV, selectionParameter_, false);
0850 spr::trackSelectionParameters oneCutParameters = selectionParameter_;
0851 oneCutParameters.maxDxyPV = 10;
0852 oneCutParameters.maxDzPV = 100;
0853 oneCutParameters.maxInMiss = 2;
0854 oneCutParameters.maxOutMiss = 2;
0855 bool qltyFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, false);
0856 oneCutParameters = selectionParameter_;
0857 oneCutParameters.maxDxyPV = 10;
0858 oneCutParameters.maxDzPV = 100;
0859 isoTk.qltyMissFlag_ = spr::goodTrack(pTrack, leadPV, oneCutParameters, false);
0860 oneCutParameters = selectionParameter_;
0861 oneCutParameters.maxInMiss = 2;
0862 oneCutParameters.maxOutMiss = 2;
0863 isoTk.qltyPVFlag_ = spr::goodTrack(pTrack, leadPV, oneCutParameters, false);
0864 double eIsolation = maxRestrictionP_ * exp(slopeRestrictionP_ * std::abs((double)isoTk.ieta_));
0865 if (eIsolation < eIsolate1_)
0866 eIsolation = eIsolate1_;
0867 if (eIsolation < eIsolate2_)
0868 eIsolation = eIsolate2_;
0869 #ifdef EDM_ML_DEBUG
0870 if (debug_)
0871 edm::LogVerbatim("HcalIsoTrack") << "qltyFlag|okECAL|okHCAL : " << qltyFlag << "|" << trkDetItr->okECAL << "|"
0872 << trkDetItr->okHCAL << " eIsolation " << eIsolation;
0873 if (qltyFlag)
0874 flag += 1;
0875 if (trkDetItr->okECAL)
0876 flag += 2;
0877 if (trkDetItr->okHCAL)
0878 flag += 4;
0879 #endif
0880 isoTk.qltyFlag_ = (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL);
0881 bool notMuon = (muonh.isValid()) ? notaMuon(pTrack, muonh) : true;
0882 #ifdef EDM_ML_DEBUG
0883 if (notMuon)
0884 flag += 8;
0885 #endif
0886 if (isoTk.qltyFlag_ && notMuon) {
0887 nselTracks++;
0888 int nNearTRKs(0);
0889
0890 std::vector<DetId> eIds;
0891 std::vector<double> eHit;
0892 #ifdef EDM_ML_DEBUG
0893 double eMipDR =
0894 #endif
0895 spr::eCone_ecal(geo,
0896 barrelRecHitsHandle,
0897 endcapRecHitsHandle,
0898 trkDetItr->pointHCAL,
0899 trkDetItr->pointECAL,
0900 a_mipR_,
0901 trkDetItr->directionECAL,
0902 eIds,
0903 eHit);
0904 double eEcal(0);
0905 for (unsigned int k = 0; k < eIds.size(); ++k) {
0906 if (eHit[k] > eThreshold(eIds[k], geo))
0907 eEcal += eHit[k];
0908 }
0909 #ifdef EDM_ML_DEBUG
0910 if (debug_)
0911 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR << ":" << eEcal;
0912 #endif
0913 isoTk.eMipDR_.emplace_back(eEcal);
0914
0915
0916 std::vector<DetId> eIds2;
0917 std::vector<double> eHit2;
0918 #ifdef EDM_ML_DEBUG
0919 double eMipDR2 =
0920 #endif
0921 spr::eCone_ecal(geo,
0922 barrelRecHitsHandle,
0923 endcapRecHitsHandle,
0924 trkDetItr->pointHCAL,
0925 trkDetItr->pointECAL,
0926 a_mipR2_,
0927 trkDetItr->directionECAL,
0928 eIds2,
0929 eHit2);
0930 double eEcal2(0);
0931 for (unsigned int k = 0; k < eIds2.size(); ++k) {
0932 if (eHit2[k] > eThreshold(eIds2[k], geo))
0933 eEcal2 += eHit2[k];
0934 }
0935 #ifdef EDM_ML_DEBUG
0936 if (debug_)
0937 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR2 << ":" << eEcal2;
0938 #endif
0939 isoTk.eMipDR_.emplace_back(eEcal2);
0940
0941
0942 std::vector<DetId> eIds3;
0943 std::vector<double> eHit3;
0944 #ifdef EDM_ML_DEBUG
0945 double eMipDR3 =
0946 #endif
0947 spr::eCone_ecal(geo,
0948 barrelRecHitsHandle,
0949 endcapRecHitsHandle,
0950 trkDetItr->pointHCAL,
0951 trkDetItr->pointECAL,
0952 a_mipR3_,
0953 trkDetItr->directionECAL,
0954 eIds3,
0955 eHit3);
0956 double eEcal3(0);
0957 for (unsigned int k = 0; k < eIds3.size(); ++k) {
0958 if (eHit3[k] > eThreshold(eIds3[k], geo))
0959 eEcal3 += eHit3[k];
0960 }
0961 #ifdef EDM_ML_DEBUG
0962 if (debug_)
0963 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR3 << ":" << eEcal3;
0964 #endif
0965 isoTk.eMipDR_.emplace_back(eEcal3);
0966
0967
0968 std::vector<DetId> eIds4;
0969 std::vector<double> eHit4;
0970 #ifdef EDM_ML_DEBUG
0971 double eMipDR4 =
0972 #endif
0973 spr::eCone_ecal(geo,
0974 barrelRecHitsHandle,
0975 endcapRecHitsHandle,
0976 trkDetItr->pointHCAL,
0977 trkDetItr->pointECAL,
0978 a_mipR4_,
0979 trkDetItr->directionECAL,
0980 eIds4,
0981 eHit4);
0982 double eEcal4(0);
0983 for (unsigned int k = 0; k < eIds4.size(); ++k) {
0984 if (eHit4[k] > eThreshold(eIds4[k], geo))
0985 eEcal4 += eHit4[k];
0986 }
0987 #ifdef EDM_ML_DEBUG
0988 if (debug_)
0989 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR4 << ":" << eEcal4;
0990 #endif
0991 isoTk.eMipDR_.emplace_back(eEcal4);
0992
0993
0994 std::vector<DetId> eIds5;
0995 std::vector<double> eHit5;
0996 #ifdef EDM_ML_DEBUG
0997 double eMipDR5 =
0998 #endif
0999 spr::eCone_ecal(geo,
1000 barrelRecHitsHandle,
1001 endcapRecHitsHandle,
1002 trkDetItr->pointHCAL,
1003 trkDetItr->pointECAL,
1004 a_mipR5_,
1005 trkDetItr->directionECAL,
1006 eIds5,
1007 eHit5);
1008 double eEcal5(0);
1009 for (unsigned int k = 0; k < eIds5.size(); ++k) {
1010 if (eHit5[k] > eThreshold(eIds5[k], geo))
1011 eEcal5 += eHit5[k];
1012 }
1013 #ifdef EDM_ML_DEBUG
1014 if (debug_)
1015 edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR5 << ":" << eEcal5;
1016 #endif
1017 isoTk.eMipDR_.emplace_back(eEcal5);
1018
1019
1020 isoTk.emaxNearP_ = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 15, 15);
1021 const DetId cellE(trkDetItr->detIdECAL);
1022 std::pair<double, bool> e11x11P = spr::eECALmatrix(cellE,
1023 barrelRecHitsHandle,
1024 endcapRecHitsHandle,
1025 *theEcalChStatus,
1026 geo,
1027 caloTopology,
1028 theEcalSevlv,
1029 5,
1030 5,
1031 -100.0,
1032 -100.0,
1033 -100.0,
1034 100.0);
1035 std::pair<double, bool> e15x15P = spr::eECALmatrix(cellE,
1036 barrelRecHitsHandle,
1037 endcapRecHitsHandle,
1038 *theEcalChStatus,
1039 geo,
1040 caloTopology,
1041 theEcalSevlv,
1042 7,
1043 7,
1044 -100.0,
1045 -100.0,
1046 -100.0,
1047 100.0);
1048 if (e11x11P.second && e15x15P.second) {
1049 isoTk.eAnnular_ = (e15x15P.first - e11x11P.first);
1050 } else {
1051 isoTk.eAnnular_ = -(e15x15P.first - e11x11P.first);
1052 }
1053 isoTk.hmaxNearP_ = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, false);
1054 const DetId cellH(trkDetItr->detIdHCAL);
1055 double h5x5 = spr::eHCALmatrix(
1056 theHBHETopology, cellH, hbhe, 2, 2, false, true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1057 double h7x7 = spr::eHCALmatrix(
1058 theHBHETopology, cellH, hbhe, 3, 3, false, true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1059 isoTk.hAnnular_ = h7x7 - h5x5;
1060 #ifdef EDM_ML_DEBUG
1061 if (debug_)
1062 edm::LogVerbatim("HcalIsoTrack") << "max p Near (Ecal) " << isoTk.emaxNearP_ << " (Hcal) " << isoTk.hmaxNearP_
1063 << " Annular E (Ecal) " << e11x11P.first << ":" << e15x15P.first << ":"
1064 << isoTk.eAnnular_ << " (Hcal) " << h5x5 << ":" << h7x7 << ":"
1065 << isoTk.hAnnular_;
1066 if (isoTk.eMipDR_[0] < eEcalMax_)
1067 flag += 16;
1068 if (isoTk.hmaxNearP_ < eIsolation)
1069 flag += 32;
1070 #endif
1071 isoTk.gentrackP_ = trackP(pTrack, genParticles);
1072 if (isoTk.eMipDR_[0] < eEcalMax_ && isoTk.hmaxNearP_ < eIsolation) {
1073 int nRecHits(-999), nRecHits1(-999), nRecHits3(-999);
1074 std::vector<DetId> ids, ids1, ids3;
1075 std::vector<double> edet0, edet1, edet3;
1076 isoTk.eHcal_ = spr::eCone_hcal(geo,
1077 hbhe,
1078 trkDetItr->pointHCAL,
1079 trkDetItr->pointECAL,
1080 a_coneR_,
1081 trkDetItr->directionHCAL,
1082 nRecHits,
1083 ids,
1084 edet0,
1085 0);
1086 if (!oldID_.empty()) {
1087 for (unsigned k = 0; k < ids.size(); ++k)
1088 ids[k] = newId(ids[k]);
1089 }
1090 storeEnergy(respCorrs, ids, edet0, isoTk.eHcal_, isoTk.detIds_, isoTk.hitEnergies_);
1091 std::pair<double, double> ehcal0 =
1092 storeEnergy(respCorrs, hbhe, ids, isoTk.hitEnergiesRaw_, isoTk.hitEnergiesAux_);
1093 isoTk.eHcalRaw_ = ehcal0.first;
1094 isoTk.eHcalAux_ = ehcal0.second;
1095
1096
1097 isoTk.eHcal10_ = spr::eCone_hcal(geo,
1098 hbhe,
1099 trkDetItr->pointHCAL,
1100 trkDetItr->pointECAL,
1101 a_coneR1_,
1102 trkDetItr->directionHCAL,
1103 nRecHits1,
1104 ids1,
1105 edet1,
1106 0);
1107 if (!oldID_.empty()) {
1108 for (unsigned k = 0; k < ids1.size(); ++k)
1109 ids1[k] = newId(ids1[k]);
1110 }
1111 storeEnergy(respCorrs, ids1, edet1, isoTk.eHcal10_, isoTk.detIds1_, isoTk.hitEnergies1_);
1112 std::pair<double, double> ehcal1 =
1113 storeEnergy(respCorrs, hbhe, ids1, isoTk.hitEnergies1Raw_, isoTk.hitEnergies1Aux_);
1114 isoTk.eHcal10Raw_ = ehcal1.first;
1115 isoTk.eHcal10Aux_ = ehcal1.second;
1116
1117
1118 isoTk.eHcal30_ = spr::eCone_hcal(geo,
1119 hbhe,
1120 trkDetItr->pointHCAL,
1121 trkDetItr->pointECAL,
1122 a_coneR2_,
1123 trkDetItr->directionHCAL,
1124 nRecHits3,
1125 ids3,
1126 edet3,
1127 0);
1128 if (!oldID_.empty()) {
1129 for (unsigned k = 0; k < ids3.size(); ++k)
1130 ids3[k] = newId(ids3[k]);
1131 }
1132 storeEnergy(respCorrs, ids3, edet3, isoTk.eHcal30_, isoTk.detIds3_, isoTk.hitEnergies3_);
1133 std::pair<double, double> ehcal3 =
1134 storeEnergy(respCorrs, hbhe, ids3, isoTk.hitEnergies3Raw_, isoTk.hitEnergies3Aux_);
1135 isoTk.eHcal30Raw_ = ehcal3.first;
1136 isoTk.eHcal30Aux_ = ehcal3.second;
1137
1138 if (isoTk.p_ > pTrackMin_)
1139 accept = true;
1140 #ifdef EDM_ML_DEBUG
1141 if (accept)
1142 flag += 64;
1143 if (debug_) {
1144 std::string ctype = accept ? " ***** ACCEPT *****" : "";
1145 edm::LogVerbatim("HcalIsoTrack")
1146 << "This track : " << nTracks << " (pt|eta|phi|p) : " << isoTk.pt_ << "|" << pTrack->eta() << "|"
1147 << isoTk.phi_ << "|" << isoTk.p_ << " Generator Level p " << isoTk.gentrackP_;
1148 edm::LogVerbatim("HcalIsoTrack")
1149 << "e_MIP " << isoTk.eMipDR_[0] << " Chg Isolation " << isoTk.hmaxNearP_ << " eHcal" << isoTk.eHcal_
1150 << ":" << isoTk.eHcalRaw_ << ":" << isoTk.eHcalAux_ << " ieta " << isoTk.ieta_ << " Quality "
1151 << isoTk.qltyMissFlag_ << ":" << isoTk.qltyPVFlag_ << ":" << isoTk.selectTk_ << ctype;
1152 for (unsigned int ll = 0; ll < isoTk.detIds_.size(); ll++) {
1153 edm::LogVerbatim("HcalIsoTrack")
1154 << "det id is = " << HcalDetId(isoTk.detIds_[ll]) << " hit enery is = " << isoTk.hitEnergies_[ll]
1155 << " : " << isoTk.hitEnergiesRaw_[ll] << " : " << isoTk.hitEnergiesAux_[ll];
1156 }
1157 for (unsigned int ll = 0; ll < isoTk.detIds1_.size(); ll++) {
1158 edm::LogVerbatim("HcalIsoTrack")
1159 << "det id is = " << HcalDetId(isoTk.detIds1_[ll]) << " hit enery is = " << isoTk.hitEnergies1_[ll]
1160 << " : " << isoTk.hitEnergies1Raw_[ll] << " : " << isoTk.hitEnergies1Aux_[ll];
1161 }
1162 for (unsigned int ll = 0; ll < isoTk.detIds3_.size(); ll++) {
1163 edm::LogVerbatim("HcalIsoTrack")
1164 << "det id is = " << HcalDetId(isoTk.detIds3_[ll]) << " hit enery is = " << isoTk.hitEnergies3_[ll]
1165 << " : " << isoTk.hitEnergies3Raw_[ll] << " : " << isoTk.hitEnergies3Aux_[ll];
1166 }
1167 }
1168 #endif
1169 if (accept) {
1170 edm::LogVerbatim("HcalIsoTrackX")
1171 << "Run " << eventId.run() << " Event " << eventId.event() << " Track " << nTracks << " p " << isoTk.p_;
1172 hocalib.emplace_back(isoTk);
1173 nSave++;
1174 int type(0);
1175 if (isoTk.eMipDR_[0] < 1.0) {
1176 if (isoTk.hmaxNearP_ < eIsolate2_) {
1177 ++nLoose;
1178 type = 1;
1179 }
1180 if (isoTk.hmaxNearP_ < eIsolate1_) {
1181 ++nTight;
1182 type = 2;
1183 }
1184 }
1185 if (isoTk.p_ > 40.0 && isoTk.p_ <= 60.0 && isoTk.selectTk_) {
1186 hocalibEvent.ietaGood_.emplace_back(isoTk.ieta_);
1187 hocalibEvent.trackType_.emplace_back(type);
1188 }
1189 #ifdef EDM_ML_DEBUG
1190 for (unsigned int k = 0; k < isoTk.trgbits_.size(); k++) {
1191 edm::LogVerbatim("HcalIsoTrack") << "trigger bit is = " << isoTk.trgbits_[k];
1192 }
1193 #endif
1194 }
1195 }
1196 }
1197 #ifdef EDM_ML_DEBUG
1198 if (debug_) {
1199 if (isoTk.eMipDR_.empty())
1200 edm::LogVerbatim("HcalIsoTrack") << "Track " << nTracks << " Selection Flag " << std::hex << flag << std::dec
1201 << " Accept " << accept << " Momentum " << isoTk.p_ << ":" << pTrackMin_;
1202 else
1203 edm::LogVerbatim("HcalIsoTrack") << "Track " << nTracks << " Selection Flag " << std::hex << flag << std::dec
1204 << " Accept " << accept << " Momentum " << isoTk.p_ << ":" << pTrackMin_
1205 << " Ecal Energy " << isoTk.eMipDR_[0] << ":" << eEcalMax_
1206 << " Charge Isolation " << isoTk.hmaxNearP_ << ":" << eIsolation;
1207 }
1208 #endif
1209 }
1210 std::array<int, 3> i3{{nSave, nLoose, nTight}};
1211 return i3;
1212 }
1213
1214 double AlCaHcalIsotrkProducer::dR(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
1215 return reco::deltaR(vec1.eta(), vec1.phi(), vec2.eta(), vec2.phi());
1216 }
1217
1218 double AlCaHcalIsotrkProducer::trackP(const reco::Track* pTrack,
1219 const edm::Handle<reco::GenParticleCollection>& genParticles) {
1220 double pmom = -1.0;
1221 if (genParticles.isValid()) {
1222 double mindR(999.9);
1223 for (const auto& p : (*genParticles)) {
1224 double dR = reco::deltaR(pTrack->eta(), pTrack->phi(), p.momentum().Eta(), p.momentum().Phi());
1225 if (dR < mindR) {
1226 mindR = dR;
1227 pmom = p.momentum().R();
1228 }
1229 }
1230 }
1231 return pmom;
1232 }
1233
1234 double AlCaHcalIsotrkProducer::rhoh(const edm::Handle<CaloTowerCollection>& tower) {
1235 std::vector<double> sumPFNallSMDQH2;
1236 sumPFNallSMDQH2.reserve(phibins_.size() * etabins_.size());
1237 for (auto eta : etabins_) {
1238 for (auto phi : phibins_) {
1239 double hadder = 0;
1240 for (const auto& pf_it : (*tower)) {
1241 if (fabs(eta - pf_it.eta()) > etahalfdist_)
1242 continue;
1243 if (fabs(reco::deltaPhi(phi, pf_it.phi())) > phihalfdist_)
1244 continue;
1245 hadder += pf_it.hadEt();
1246 }
1247 sumPFNallSMDQH2.emplace_back(hadder);
1248 }
1249 }
1250
1251 double evt_smdq(0);
1252 std::sort(sumPFNallSMDQH2.begin(), sumPFNallSMDQH2.end());
1253 if (sumPFNallSMDQH2.size() % 2)
1254 evt_smdq = sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 1) / 2];
1255 else
1256 evt_smdq = (sumPFNallSMDQH2[sumPFNallSMDQH2.size() / 2] + sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 2) / 2]) / 2.;
1257 double rhoh = evt_smdq / (etadist_ * phidist_);
1258 #ifdef EDM_ML_DEBUG
1259 if (debug_)
1260 edm::LogVerbatim("HcalIsoTrack") << "Rho " << evt_smdq << ":" << rhoh;
1261 #endif
1262 return rhoh;
1263 }
1264
1265 double AlCaHcalIsotrkProducer::eThreshold(const DetId& id, const CaloGeometry* geo) const {
1266 const GlobalPoint& pos = geo->getPosition(id);
1267 double eta = std::abs(pos.eta());
1268 double eThr(hitEthrEB_);
1269 if (usePFThresh_) {
1270 eThr = static_cast<double>((*eThresholds_)[id]);
1271 } else {
1272 if (id.subdetId() != EcalBarrel) {
1273 eThr = (((eta * hitEthrEE3_ + hitEthrEE2_) * eta + hitEthrEE1_) * eta + hitEthrEE0_);
1274 if (eThr < hitEthrEELo_)
1275 eThr = hitEthrEELo_;
1276 else if (eThr > hitEthrEEHi_)
1277 eThr = hitEthrEEHi_;
1278 }
1279 }
1280 return eThr;
1281 }
1282
1283 DetId AlCaHcalIsotrkProducer::newId(const DetId& id) {
1284 HcalDetId hid(id);
1285 if (hep17_ && ((hid.iphi() < 63) || (hid.iphi() > 66) || (hid.zside() < 0)))
1286 return id;
1287 for (unsigned int k = 0; k < oldID_.size(); ++k) {
1288 if ((hid.subdetId() == oldDet_[k]) && (hid.ietaAbs() == oldEta_[k]) && (hid.depth() == oldDepth_[k])) {
1289 return static_cast<DetId>(HcalDetId(hid.subdet(), hid.ieta(), hid.iphi(), newDepth_[k]));
1290 }
1291 }
1292 return id;
1293 }
1294
1295 void AlCaHcalIsotrkProducer::storeEnergy(const HcalRespCorrs* respCorrs,
1296 const std::vector<DetId>& ids,
1297 std::vector<double>& edet,
1298 double& eHcal,
1299 std::vector<unsigned int>& detIds,
1300 std::vector<double>& hitEnergies) {
1301 double ehcal(0);
1302 if (unCorrect_) {
1303 for (unsigned int k = 0; k < ids.size(); ++k) {
1304 double corr = (respCorrs->getValues(ids[k]))->getValue();
1305 if (corr != 0)
1306 edet[k] /= corr;
1307 ehcal += edet[k];
1308 }
1309 } else {
1310 for (const auto& en : edet)
1311 ehcal += en;
1312 }
1313 if ((std::abs(ehcal - eHcal) > 0.001) && (!unCorrect_))
1314 edm::LogWarning("HcalIsoTrack") << "Check inconsistent energies: " << eHcal << ":" << ehcal << " from "
1315 << ids.size() << " cells";
1316 eHcal = hcalScale_ * ehcal;
1317
1318 if (collapseDepth_) {
1319 std::map<HcalDetId, double> hitMap;
1320 for (unsigned int k = 0; k < ids.size(); ++k) {
1321 HcalDetId id = hdc_->mergedDepthDetId(HcalDetId(ids[k]));
1322 auto itr = hitMap.find(id);
1323 if (itr == hitMap.end()) {
1324 hitMap[id] = edet[k];
1325 } else {
1326 (itr->second) += edet[k];
1327 }
1328 }
1329 detIds.reserve(hitMap.size());
1330 hitEnergies.reserve(hitMap.size());
1331 for (const auto& hit : hitMap) {
1332 detIds.emplace_back(hit.first.rawId());
1333 hitEnergies.emplace_back(hit.second);
1334 }
1335 } else {
1336 detIds.reserve(ids.size());
1337 hitEnergies.reserve(ids.size());
1338 for (unsigned int k = 0; k < ids.size(); ++k) {
1339 detIds.emplace_back(ids[k].rawId());
1340 hitEnergies.emplace_back(edet[k]);
1341 }
1342 }
1343 #ifdef EDM_ML_DEBUG
1344 if (debug_) {
1345 edm::LogVerbatim("HcalIsoTrack") << "StoreEnergy1::Input to storeEnergy with " << ids.size() << " cells";
1346 for (unsigned int k = 0; k < ids.size(); ++k)
1347 edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << HcalDetId(ids[k]) << " E " << edet[k];
1348 edm::LogVerbatim("HcalIsoTrack") << "StoreEnergy1::Output of storeEnergy with " << detIds.size()
1349 << " cells and Etot " << eHcal;
1350 for (unsigned int k = 0; k < detIds.size(); ++k)
1351 edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << HcalDetId(detIds[k]) << " E " << hitEnergies[k];
1352 }
1353 #endif
1354 }
1355
1356 std::pair<double, double> AlCaHcalIsotrkProducer::storeEnergy(const HcalRespCorrs* respCorrs,
1357 edm::Handle<HBHERecHitCollection>& hbhe,
1358 const std::vector<DetId>& ids,
1359 std::vector<double>& hitEnergy1,
1360 std::vector<double>& hitEnergy2) {
1361 double ehcal1(0), ehcal2(0);
1362 std::vector<double> edet1, edet2;
1363 for (unsigned int k = 0; k < ids.size(); ++k) {
1364 double e1(0), e2(0);
1365 for (auto itr = hbhe->begin(); itr != hbhe->end(); ++itr) {
1366 if (itr->id() == ids[k]) {
1367 e1 = itr->eraw();
1368 e2 = itr->eaux();
1369 break;
1370 }
1371 }
1372 if (e1 < 1.e-10)
1373 e1 = 0;
1374 if (e2 < 1.e-10)
1375 e2 = 0;
1376 edet1.emplace_back(e1);
1377 edet2.emplace_back(e2);
1378 }
1379 if (unCorrect_) {
1380 for (unsigned int k = 0; k < ids.size(); ++k) {
1381 double corr = (respCorrs->getValues(ids[k]))->getValue();
1382 if (corr != 0) {
1383 edet1[k] /= corr;
1384 edet2[k] /= corr;
1385 }
1386 }
1387 }
1388 for (unsigned int k = 0; k < ids.size(); ++k) {
1389 ehcal1 += edet1[k];
1390 ehcal2 += edet2[k];
1391 }
1392 ehcal1 *= hcalScale_;
1393 ehcal2 *= hcalScale_;
1394
1395 if (collapseDepth_) {
1396 std::map<HcalDetId, std::pair<double, double>> hitMap;
1397 for (unsigned int k = 0; k < ids.size(); ++k) {
1398 HcalDetId id = hdc_->mergedDepthDetId(HcalDetId(ids[k]));
1399 auto itr = hitMap.find(id);
1400 if (itr == hitMap.end()) {
1401 hitMap[id] = std::make_pair(edet1[k], edet2[k]);
1402 } else {
1403 (itr->second).first += edet1[k];
1404 (itr->second).second += edet2[k];
1405 }
1406 }
1407 hitEnergy1.reserve(hitMap.size());
1408 hitEnergy2.reserve(hitMap.size());
1409 for (const auto& hit : hitMap) {
1410 hitEnergy1.emplace_back(hit.second.first);
1411 hitEnergy2.emplace_back(hit.second.second);
1412 }
1413 } else {
1414 hitEnergy1.reserve(ids.size());
1415 hitEnergy2.reserve(ids.size());
1416 for (unsigned int k = 0; k < ids.size(); ++k) {
1417 hitEnergy1.emplace_back(edet1[k]);
1418 hitEnergy2.emplace_back(edet2[k]);
1419 }
1420 }
1421 #ifdef EDM_ML_DEBUG
1422 if (debug_) {
1423 edm::LogVerbatim("HcalIsoTrack") << "StoreEnergy2::Input to storeEnergy with " << ids.size() << " cells";
1424 edm::LogVerbatim("HcalIsoTrack") << "StoreEnergy2::Output of storeEnergy with " << hitEnergy1.size()
1425 << " cells and Etot " << ehcal1 << ":" << ehcal2;
1426 for (unsigned int k = 0; k < hitEnergy1.size(); ++k)
1427 edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << hitEnergy1[k] << " : " << hitEnergy2[k];
1428 }
1429 #endif
1430 return std::make_pair(ehcal1, ehcal2);
1431 }
1432
1433 bool AlCaHcalIsotrkProducer::notaMuon(const reco::Track* pTrack0, const edm::Handle<reco::MuonCollection>& muonh) {
1434 bool flag(true);
1435 for (reco::MuonCollection::const_iterator recMuon = muonh->begin(); recMuon != muonh->end(); ++recMuon) {
1436 if (recMuon->innerTrack().isNonnull()) {
1437 const reco::Track* pTrack = (recMuon->innerTrack()).get();
1438 bool mediumMuon = (((recMuon->isPFMuon()) && (recMuon->isGlobalMuon() || recMuon->isTrackerMuon())) &&
1439 (recMuon->innerTrack()->validFraction() > 0.49));
1440 if (mediumMuon) {
1441 double chiGlobal = ((recMuon->globalTrack().isNonnull()) ? recMuon->globalTrack()->normalizedChi2() : 999);
1442 bool goodGlob = (recMuon->isGlobalMuon() && chiGlobal < 3 &&
1443 recMuon->combinedQuality().chi2LocalPosition < 12 && recMuon->combinedQuality().trkKink < 20);
1444 mediumMuon = muon::segmentCompatibility(*recMuon) > (goodGlob ? 0.303 : 0.451);
1445 }
1446 if (mediumMuon) {
1447 double dR = reco::deltaR(pTrack->eta(), pTrack->phi(), pTrack0->eta(), pTrack0->phi());
1448 if (dR < 0.1) {
1449 flag = false;
1450 break;
1451 }
1452 }
1453 }
1454 }
1455 return flag;
1456 }
1457
1458 #include "FWCore/Framework/interface/MakerMacros.h"
1459
1460 DEFINE_FWK_MODULE(AlCaHcalIsotrkProducer);