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