File indexing completed on 2024-04-06 11:58:44
0001 #include <memory>
0002 #include <iostream>
0003 #include <fstream>
0004 #include <vector>
0005 #include <TFile.h>
0006 #include <TTree.h>
0007 #include "TPRegexp.h"
0008
0009
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/stream/EDProducer.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0016 #include "FWCore/ServiceRegistry/interface/Service.h"
0017 #include "FWCore/Common/interface/TriggerNames.h"
0018
0019 #include "DataFormats/HcalCalibObjects/interface/HcalHBHEMuonVariables.h"
0020 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0021 #include "DataFormats/MuonReco/interface/Muon.h"
0022 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0023 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0024 #include "DataFormats/VertexReco/interface/Vertex.h"
0025 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027
0028
0029
0030 #include "DataFormats/Common/interface/TriggerResults.h"
0031 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0032 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0033 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0034 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0035 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0036 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0037 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0038 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0039
0040 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0041 #include "HLTrigger/HLTcore/interface/HLTConfigData.h"
0042
0043 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0044 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0045
0046 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0047 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0048 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0049 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0050 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0051
0052 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0053 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0054 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0055
0056 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0057 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0058 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0059
0060 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0061 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0062 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0063 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0064 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0065 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0066 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0067 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0068 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0069 #include "MagneticField/Engine/interface/MagneticField.h"
0070 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0071
0072
0073
0074 namespace alcaHcalHBHEMuonProducer {
0075 struct Counters {
0076 Counters() : nAll_(0), nGood_(0) {}
0077 mutable std::atomic<unsigned int> nAll_, nGood_;
0078 };
0079 }
0080
0081 class AlCaHcalHBHEMuonProducer : public edm::stream::EDProducer<edm::GlobalCache<alcaHcalHBHEMuonProducer::Counters>> {
0082 public:
0083 explicit AlCaHcalHBHEMuonProducer(const edm::ParameterSet&, const alcaHcalHBHEMuonProducer::Counters*);
0084 ~AlCaHcalHBHEMuonProducer() override = default;
0085
0086 static std::unique_ptr<alcaHcalHBHEMuonProducer::Counters> initializeGlobalCache(edm::ParameterSet const&) {
0087 return std::make_unique<alcaHcalHBHEMuonProducer::Counters>();
0088 }
0089
0090 void produce(edm::Event&, const edm::EventSetup&) override;
0091
0092 void endStream() override;
0093
0094 static void globalEndJob(const alcaHcalHBHEMuonProducer::Counters* counters);
0095
0096 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0097
0098 private:
0099 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0100 void endRun(edm::Run const&, edm::EventSetup const&) override;
0101 int matchId(const HcalDetId&, const HcalDetId&);
0102 double activeLength(const DetId&, const HcalDDDRecConstants* hdc);
0103 bool isGoodVertex(const reco::Vertex& vtx);
0104 double respCorr(const DetId& id, const HcalRespCorrs* respCorrs);
0105 double gainFactor(const HcalDbService* dbserv, const HcalDetId& id);
0106 int depth16HE(int ieta, int iphi, const HcalTopology* theHBHETopology);
0107 bool goodCell(const HcalDetId& hcid,
0108 const reco::Track* pTrack,
0109 const CaloGeometry* geo,
0110 const MagneticField* bField,
0111 const HcalDDDRecConstants* hdc);
0112
0113
0114 HLTConfigProvider hltConfig_;
0115 const std::vector<std::string> trigNames_;
0116 const std::string processName_;
0117 const edm::InputTag triggerResults_;
0118 const edm::InputTag labelEBRecHit_, labelEERecHit_, labelHBHERecHit_;
0119 const std::string labelVtx_, labelMuon_, labelHBHEMuon_;
0120 const bool collapseDepth_, isItPlan1_;
0121 const int verbosity_;
0122 const bool isItPreRecHit_, writeRespCorr_;
0123 const std::string fileInCorr_;
0124 const int maxDepth_;
0125 const bool mergedDepth_, usePFThresh_;
0126
0127 bool useMyCorr_;
0128 int nRun_, nAll_, nGood_;
0129
0130 const edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0131 const edm::EDGetTokenT<reco::VertexCollection> tok_Vtx_;
0132 const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0133 const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0134 const edm::EDGetTokenT<HBHERecHitCollection> tok_HBHE_;
0135 const edm::EDGetTokenT<reco::MuonCollection> tok_Muon_;
0136
0137 const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec0_;
0138 const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec1_;
0139 const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_respcorr0_;
0140 const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_respcorr1_;
0141 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom0_;
0142 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom1_;
0143 const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo0_;
0144 const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo1_;
0145 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0146 const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_chan_;
0147 const edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0148 const edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_topo_;
0149 const edm::ESGetToken<HcalDbService, HcalDbRecord> tok_dbservice_;
0150 const edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> tok_ecalPFRecHitThresholds_;
0151
0152
0153 static const int depthMax_ = 7;
0154 std::vector<std::string> all_triggers_;
0155 std::vector<int> hltresults_;
0156
0157 std::vector<HcalDDDRecConstants::HcalActiveLength> actHB, actHE;
0158 std::map<DetId, double> corrValue_;
0159
0160 };
0161
0162 AlCaHcalHBHEMuonProducer::AlCaHcalHBHEMuonProducer(const edm::ParameterSet& iConfig,
0163 const alcaHcalHBHEMuonProducer::Counters*)
0164 : trigNames_(iConfig.getParameter<std::vector<std::string>>("triggers")),
0165 processName_(iConfig.getParameter<std::string>("processName")),
0166 triggerResults_(iConfig.getParameter<edm::InputTag>("triggerResults")),
0167 labelEBRecHit_(iConfig.getParameter<edm::InputTag>("labelEBRecHit")),
0168 labelEERecHit_(iConfig.getParameter<edm::InputTag>("labelEERecHit")),
0169 labelHBHERecHit_(iConfig.getParameter<edm::InputTag>("labelHBHERecHit")),
0170 labelVtx_(iConfig.getParameter<std::string>("labelVertex")),
0171 labelMuon_(iConfig.getParameter<std::string>("labelMuon")),
0172 labelHBHEMuon_(iConfig.getParameter<std::string>("labelHBHEMuon")),
0173 collapseDepth_(iConfig.getParameter<bool>("collapseDepth")),
0174 isItPlan1_(iConfig.getParameter<bool>("isItPlan1")),
0175 verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0176 isItPreRecHit_(iConfig.getUntrackedParameter<bool>("isItPreRecHit", false)),
0177 writeRespCorr_(iConfig.getUntrackedParameter<bool>("writeRespCorr", false)),
0178 fileInCorr_(iConfig.getUntrackedParameter<std::string>("fileInCorr", "")),
0179 maxDepth_(iConfig.getUntrackedParameter<int>("maxDepth", 7)),
0180 mergedDepth_((!isItPreRecHit_) || (collapseDepth_)),
0181 usePFThresh_(iConfig.getParameter<bool>("usePFThreshold")),
0182 nRun_(0),
0183 nAll_(0),
0184 nGood_(0),
0185 tok_trigRes_(consumes<edm::TriggerResults>(triggerResults_)),
0186 tok_Vtx_(consumes<reco::VertexCollection>(labelVtx_)),
0187 tok_EB_(consumes<EcalRecHitCollection>(labelEBRecHit_)),
0188 tok_EE_(consumes<EcalRecHitCollection>(labelEERecHit_)),
0189 tok_HBHE_(consumes<HBHERecHitCollection>(labelHBHERecHit_)),
0190 tok_Muon_(consumes<reco::MuonCollection>(labelMuon_)),
0191 tok_ddrec0_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0192 tok_ddrec1_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>()),
0193 tok_respcorr0_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd, edm::Transition::BeginRun>()),
0194 tok_respcorr1_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd>()),
0195 tok_geom0_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0196 tok_geom1_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0197 tok_htopo0_(esConsumes<HcalTopology, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0198 tok_htopo1_(esConsumes<HcalTopology, HcalRecNumberingRecord>()),
0199 tok_magField_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0200 tok_chan_(esConsumes<EcalChannelStatus, EcalChannelStatusRcd>()),
0201 tok_sevlv_(esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>()),
0202 tok_topo_(esConsumes<CaloTopology, CaloTopologyRecord>()),
0203 tok_dbservice_(esConsumes<HcalDbService, HcalDbRecord>()),
0204 tok_ecalPFRecHitThresholds_(esConsumes<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd>()) {
0205
0206 edm::LogVerbatim("HBHEMuon") << "Labels used: Trig " << triggerResults_ << " Vtx " << labelVtx_ << " EB "
0207 << labelEBRecHit_ << " EE " << labelEERecHit_ << " HBHE " << labelHBHERecHit_ << " MU "
0208 << labelMuon_;
0209
0210 if (!fileInCorr_.empty()) {
0211 std::ifstream infile(fileInCorr_.c_str());
0212 if (infile.is_open()) {
0213 while (true) {
0214 unsigned int id;
0215 double cfac;
0216 infile >> id >> cfac;
0217 if (!infile.good())
0218 break;
0219 corrValue_[DetId(id)] = cfac;
0220 }
0221 infile.close();
0222 }
0223 }
0224 useMyCorr_ = (!corrValue_.empty());
0225 edm::LogVerbatim("HBHEMuon") << "Flags used: ollapseDepth " << collapseDepth_ << ":" << mergedDepth_ << " IsItPlan1 "
0226 << isItPlan1_ << " IsItPreRecHit " << isItPreRecHit_ << " UseMyCorr " << useMyCorr_;
0227
0228
0229 produces<HcalHBHEMuonVariablesCollection>(labelHBHEMuon_);
0230 edm::LogVerbatim("HcalIsoTrack") << " Expected to produce the collections:\n"
0231 << "HcalHBHEMuonVariablesCollection with label " << labelHBHEMuon_;
0232 }
0233
0234
0235
0236
0237
0238
0239 void AlCaHcalHBHEMuonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0240 ++nAll_;
0241 auto outputHcalHBHEMuonColl = std::make_unique<HcalHBHEMuonVariablesCollection>();
0242
0243 unsigned int runNumber = iEvent.id().run();
0244 unsigned int eventNumber = iEvent.id().event();
0245 unsigned int lumiNumber = iEvent.id().luminosityBlock();
0246 unsigned int bxNumber = iEvent.bunchCrossing();
0247 #ifdef EDM_ML_DEBUG
0248 edm::LogVerbatim("HBHEMuon") << "Run " << runNumber << " Event " << eventNumber << " Lumi " << lumiNumber << " BX "
0249 << bxNumber << std::endl;
0250 #endif
0251
0252
0253 bool ok(false);
0254
0255 if (trigNames_.empty()) {
0256 ok = true;
0257 } else {
0258 auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
0259 if (triggerResults.isValid()) {
0260 std::vector<std::string> modules;
0261 const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0262 const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
0263 for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0264 int hlt = triggerResults->accept(iHLT);
0265 for (unsigned int i = 0; i < trigNames_.size(); ++i) {
0266 if (triggerNames_[iHLT].find(trigNames_[i]) != std::string::npos) {
0267 if (hlt > 0) {
0268 ok = true;
0269 }
0270 #ifdef EDM_ML_DEBUG
0271 edm::LogVerbatim("HBHEMuon") << "AlCaHcalHBHEMuonFilter::Trigger " << triggerNames_[iHLT] << " Flag " << hlt
0272 << ":" << ok << std::endl;
0273 #endif
0274 }
0275 }
0276 }
0277 }
0278 }
0279
0280
0281 const HcalDDDRecConstants* hdc = &iSetup.getData(tok_ddrec1_);
0282 const HcalRespCorrs* resp = &iSetup.getData(tok_respcorr1_);
0283 const CaloGeometry* geo = &iSetup.getData(tok_geom1_);
0284 const HcalTopology* theHBHETopology = &iSetup.getData(tok_htopo1_);
0285 const MagneticField* bField = &iSetup.getData(tok_magField_);
0286 const EcalChannelStatus* theEcalChStatus = &iSetup.getData(tok_chan_);
0287 const EcalSeverityLevelAlgo* sevlv = &iSetup.getData(tok_sevlv_);
0288 const CaloTopology* caloTopology = &iSetup.getData(tok_topo_);
0289 const HcalDbService* conditions = &iSetup.getData(tok_dbservice_);
0290 HcalRespCorrs respCorrsObj(*resp);
0291 HcalRespCorrs* respCorrs = &respCorrsObj;
0292 respCorrs->setTopo(theHBHETopology);
0293 const EcalPFRecHitThresholds* eThresholds = &iSetup.getData(tok_ecalPFRecHitThresholds_);
0294 ;
0295
0296
0297 auto const& vtx = iEvent.getHandle(tok_Vtx_);
0298 auto barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0299 auto endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0300 auto hbhe = iEvent.getHandle(tok_HBHE_);
0301 auto const& muons = iEvent.getHandle(tok_Muon_);
0302
0303
0304 math::XYZPoint pvx;
0305 unsigned int goodVertex = 0;
0306 reco::VertexCollection::const_iterator firstGoodVertex;
0307 if (!vtx.isValid()) {
0308 #ifdef EDM_ML_DEBUG
0309 edm::LogVerbatim("HBHEMuon") << "No Good Vertex found == Reject";
0310 #endif
0311 } else {
0312 firstGoodVertex = vtx->end();
0313 for (reco::VertexCollection::const_iterator it = vtx->begin(); it != vtx->end(); it++) {
0314 if (isGoodVertex(*it)) {
0315 if (firstGoodVertex == vtx->end())
0316 firstGoodVertex = it;
0317 ++goodVertex;
0318 }
0319 }
0320 if (firstGoodVertex != vtx->end())
0321 pvx = firstGoodVertex->position();
0322 }
0323
0324 if (ok && (goodVertex > 0) && muons.isValid() && barrelRecHitsHandle.isValid() && endcapRecHitsHandle.isValid() &&
0325 hbhe.isValid()) {
0326 for (reco::MuonCollection::const_iterator recMuon = muons->begin(); recMuon != muons->end(); ++recMuon) {
0327 HcalHBHEMuonVariables hbheMuon;
0328 hbheMuon.runNumber_ = runNumber;
0329 hbheMuon.eventNumber_ = eventNumber;
0330 hbheMuon.lumiNumber_ = lumiNumber;
0331 hbheMuon.bxNumber_ = bxNumber;
0332 hbheMuon.goodVertex_ = goodVertex;
0333 hbheMuon.muonGood_ = (recMuon->isPFMuon());
0334 hbheMuon.muonGlobal_ = (recMuon->isGlobalMuon());
0335 hbheMuon.muonTracker_ = (recMuon->isTrackerMuon());
0336 hbheMuon.ptGlob_ = ((recMuon)->pt());
0337 hbheMuon.etaGlob_ = (recMuon->eta());
0338 hbheMuon.phiGlob_ = (recMuon->phi());
0339 hbheMuon.energyMuon_ = (recMuon->energy());
0340 hbheMuon.pMuon_ = (recMuon->p());
0341 #ifdef EDM_ML_DEBUG
0342 edm::LogVerbatim("HBHEMuon") << "Energy:" << recMuon->energy() << " P:" << recMuon->p();
0343 #endif
0344 hbheMuon.muonTight_ = (muon::isTightMuon(*recMuon, *firstGoodVertex));
0345 hbheMuon.muonMedium_ = (muon::isMediumMuon(*recMuon));
0346 hbheMuon.muonTrkKink_ = (recMuon->combinedQuality().trkKink);
0347 hbheMuon.muonChi2LocalPosition_ = (recMuon->combinedQuality().chi2LocalPosition);
0348 hbheMuon.muonSegComp_ = (muon::segmentCompatibility(*recMuon));
0349
0350 if (recMuon->track().isNonnull()) {
0351 hbheMuon.trackerLayer_ = (recMuon->track()->hitPattern().trackerLayersWithMeasurement());
0352 } else {
0353 hbheMuon.trackerLayer_ = -1;
0354 }
0355 if (recMuon->innerTrack().isNonnull()) {
0356 hbheMuon.innerTrack_ = true;
0357 hbheMuon.numPixelLayers_ = (recMuon->innerTrack()->hitPattern().pixelLayersWithMeasurement());
0358 hbheMuon.chiTracker_ = (recMuon->innerTrack()->normalizedChi2());
0359 hbheMuon.dxyTracker_ = (fabs(recMuon->innerTrack()->dxy(pvx)));
0360 hbheMuon.dzTracker_ = (fabs(recMuon->innerTrack()->dz(pvx)));
0361 hbheMuon.innerTrackPt_ = (recMuon->innerTrack()->pt());
0362 hbheMuon.innerTrackEta_ = (recMuon->innerTrack()->eta());
0363 hbheMuon.innerTrackPhi_ = (recMuon->innerTrack()->phi());
0364 hbheMuon.tightPixelHits_ = (recMuon->innerTrack()->hitPattern().numberOfValidPixelHits());
0365 hbheMuon.tightValidFraction_ = (recMuon->innerTrack()->validFraction());
0366 }
0367
0368 if (recMuon->outerTrack().isNonnull()) {
0369 hbheMuon.outerTrack_ = true;
0370 hbheMuon.outerTrackPt_ = (recMuon->outerTrack()->pt());
0371 hbheMuon.outerTrackEta_ = (recMuon->outerTrack()->eta());
0372 hbheMuon.outerTrackPhi_ = (recMuon->outerTrack()->phi());
0373 hbheMuon.outerTrackChi_ = (recMuon->outerTrack()->normalizedChi2());
0374 hbheMuon.outerTrackHits_ = (recMuon->outerTrack()->numberOfValidHits());
0375 hbheMuon.outerTrackRHits_ = (recMuon->outerTrack()->recHitsSize());
0376 }
0377
0378 if (recMuon->globalTrack().isNonnull()) {
0379 hbheMuon.globalTrack_ = true;
0380 hbheMuon.chiGlobal_ = (recMuon->globalTrack()->normalizedChi2());
0381 hbheMuon.globalMuonHits_ = (recMuon->globalTrack()->hitPattern().numberOfValidMuonHits());
0382 hbheMuon.matchedStat_ = (recMuon->numberOfMatchedStations());
0383 hbheMuon.globalTrackPt_ = (recMuon->globalTrack()->pt());
0384 hbheMuon.globalTrackEta_ = (recMuon->globalTrack()->eta());
0385 hbheMuon.globalTrackPhi_ = (recMuon->globalTrack()->phi());
0386 hbheMuon.tightTransImpara_ = (fabs(recMuon->muonBestTrack()->dxy(pvx)));
0387 hbheMuon.tightLongPara_ = (fabs(recMuon->muonBestTrack()->dz(pvx)));
0388 }
0389
0390 hbheMuon.isolationR04_ =
0391 ((recMuon->pfIsolationR04().sumChargedHadronPt +
0392 std::max(0.,
0393 recMuon->pfIsolationR04().sumNeutralHadronEt + recMuon->pfIsolationR04().sumPhotonEt -
0394 (0.5 * recMuon->pfIsolationR04().sumPUPt))) /
0395 recMuon->pt());
0396
0397 hbheMuon.isolationR03_ =
0398 ((recMuon->pfIsolationR03().sumChargedHadronPt +
0399 std::max(0.,
0400 recMuon->pfIsolationR03().sumNeutralHadronEt + recMuon->pfIsolationR03().sumPhotonEt -
0401 (0.5 * recMuon->pfIsolationR03().sumPUPt))) /
0402 recMuon->pt());
0403
0404 hbheMuon.ecalEnergy_ = (recMuon->calEnergy().emS9);
0405 hbheMuon.hcalEnergy_ = (recMuon->calEnergy().hadS9);
0406 hbheMuon.hoEnergy_ = (recMuon->calEnergy().hoS9);
0407
0408 if (recMuon->innerTrack().isNonnull()) {
0409 const reco::Track* pTrack = (recMuon->innerTrack()).get();
0410 spr::propagatedTrackID trackID = spr::propagateCALO(pTrack, geo, bField, (((verbosity_ / 100) % 10 > 0)));
0411
0412 double activeLengthTot(0), activeLengthHotTot(0);
0413 double eHcalDepth[depthMax_], eHcalDepthHot[depthMax_];
0414 double eHcalDepthC[depthMax_], eHcalDepthHotC[depthMax_];
0415 double cHcalDepthHot[depthMax_], cHcalDepthHotBG[depthMax_];
0416 double eHcalDepthRaw[depthMax_], eHcalDepthHotRaw[depthMax_];
0417 double eHcalDepthCRaw[depthMax_], eHcalDepthHotCRaw[depthMax_];
0418 double cHcalDepthHotRaw[depthMax_], cHcalDepthHotBGRaw[depthMax_];
0419 double eHcalDepthAux[depthMax_], eHcalDepthHotAux[depthMax_];
0420 double eHcalDepthCAux[depthMax_], eHcalDepthHotCAux[depthMax_];
0421 double cHcalDepthHotAux[depthMax_], cHcalDepthHotBGAux[depthMax_];
0422 double activeL[depthMax_], activeHotL[depthMax_];
0423 bool matchDepth[depthMax_], matchDepthHot[depthMax_];
0424 HcalDetId eHcalDetId[depthMax_];
0425 unsigned int isHot(0);
0426 int ieta(-1000), iphi(-1000);
0427 for (int i = 0; i < depthMax_; ++i) {
0428 eHcalDepth[i] = eHcalDepthHot[i] = 0;
0429 eHcalDepthC[i] = eHcalDepthHotC[i] = 0;
0430 cHcalDepthHot[i] = cHcalDepthHotBG[i] = 0;
0431 eHcalDepthRaw[i] = eHcalDepthHotRaw[i] = 0;
0432 eHcalDepthCRaw[i] = eHcalDepthHotCRaw[i] = 0;
0433 cHcalDepthHotRaw[i] = cHcalDepthHotBGRaw[i] = 0;
0434 eHcalDepthAux[i] = eHcalDepthHotAux[i] = 0;
0435 eHcalDepthCAux[i] = eHcalDepthHotCAux[i] = 0;
0436 cHcalDepthHotAux[i] = cHcalDepthHotBGAux[i] = 0;
0437 activeL[i] = activeHotL[i] = 0;
0438 matchDepth[i] = matchDepthHot[i] = true;
0439 }
0440 #ifdef EDM_ML_DEBUG
0441 double eHcal(0);
0442 #endif
0443
0444 hbheMuon.ecalDetId_ = ((trackID.detIdECAL)());
0445 hbheMuon.hcalDetId_ = ((trackID.detIdHCAL)());
0446 hbheMuon.ehcalDetId_ = ((trackID.detIdEHCAL)());
0447
0448 HcalDetId check(false);
0449 std::pair<bool, HcalDetId> info = spr::propagateHCALBack(pTrack, geo, bField, (((verbosity_ / 100) % 10 > 0)));
0450 if (info.first) {
0451 check = info.second;
0452 }
0453
0454 bool okE = trackID.okECAL;
0455 if (okE) {
0456 const DetId isoCell(trackID.detIdECAL);
0457 std::pair<double, bool> e3x3 = (usePFThresh_ ? spr::eECALmatrix(isoCell,
0458 barrelRecHitsHandle,
0459 endcapRecHitsHandle,
0460 *theEcalChStatus,
0461 geo,
0462 caloTopology,
0463 sevlv,
0464 eThresholds,
0465 1,
0466 1,
0467 false)
0468 : spr::eECALmatrix(isoCell,
0469 barrelRecHitsHandle,
0470 endcapRecHitsHandle,
0471 *theEcalChStatus,
0472 geo,
0473 caloTopology,
0474 sevlv,
0475 1,
0476 1,
0477 -100.0,
0478 -100.0,
0479 -500.0,
0480 500.0,
0481 false));
0482 hbheMuon.ecal3x3Energy_ = e3x3.first;
0483 okE = e3x3.second;
0484 }
0485 #ifdef EDM_ML_DEBUG
0486 edm::LogVerbatim("HBHEMuon") << "Propagate Track to ECAL: " << okE << ":" << trackID.okECAL << " E "
0487 << hbheMuon.ecal3x3Energy_;
0488 #endif
0489
0490 if (trackID.okHCAL) {
0491 DetId closestCell(trackID.detIdHCAL);
0492 HcalDetId hcidt(closestCell.rawId());
0493 if ((hcidt.ieta() == check.ieta()) && (hcidt.iphi() == check.iphi()))
0494 hbheMuon.matchedId_ = true;
0495 #ifdef EDM_ML_DEBUG
0496 edm::LogVerbatim("HBHEMuon") << "Front " << hcidt << " Back " << info.first << ":" << check << " Match "
0497 << hbheMuon.matchedId_;
0498 #endif
0499 HcalSubdetector subdet = hcidt.subdet();
0500 ieta = hcidt.ieta();
0501 iphi = hcidt.iphi();
0502 bool hborhe = (std::abs(ieta) == 16);
0503
0504 hbheMuon.hcal1x1Energy_ = spr::eHCALmatrix(
0505 theHBHETopology, closestCell, hbhe, 0, 0, false, true, -100.0, -100.0, -100.0, -100.0, -500., 500., 0);
0506 hbheMuon.hcal1x1EnergyAux_ = spr::eHCALmatrix(
0507 theHBHETopology, closestCell, hbhe, 0, 0, false, true, -100.0, -100.0, -100.0, -100.0, -500., 500., 1);
0508 hbheMuon.hcal1x1EnergyRaw_ = spr::eHCALmatrix(
0509 theHBHETopology, closestCell, hbhe, 0, 0, false, true, -100.0, -100.0, -100.0, -100.0, -500., 500., 2);
0510 std::vector<std::pair<double, int>> ehdepth, ehdepthAux, ehdepthRaw;
0511 spr::energyHCALCell((HcalDetId)closestCell,
0512 hbhe,
0513 ehdepth,
0514 maxDepth_,
0515 -100.0,
0516 -100.0,
0517 -100.0,
0518 -100.0,
0519 -500.0,
0520 500.0,
0521 0,
0522 depth16HE(ieta, iphi, theHBHETopology),
0523 (((verbosity_ / 1000) % 10) > 0));
0524 for (int i = 0; i < depthMax_; ++i)
0525 eHcalDetId[i] = HcalDetId();
0526 spr::energyHCALCell((HcalDetId)closestCell,
0527 hbhe,
0528 ehdepthAux,
0529 maxDepth_,
0530 -100.0,
0531 -100.0,
0532 -100.0,
0533 -100.0,
0534 -500.0,
0535 500.0,
0536 1,
0537 depth16HE(ieta, iphi, theHBHETopology),
0538 (((verbosity_ / 1000) % 10) > 0));
0539 spr::energyHCALCell((HcalDetId)closestCell,
0540 hbhe,
0541 ehdepthRaw,
0542 maxDepth_,
0543 -100.0,
0544 -100.0,
0545 -100.0,
0546 -100.0,
0547 -500.0,
0548 500.0,
0549 2,
0550 depth16HE(ieta, iphi, theHBHETopology),
0551 (((verbosity_ / 1000) % 10) > 0));
0552 for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0553 HcalSubdetector subdet0 =
0554 (hborhe) ? ((ehdepth[i].second >= depth16HE(ieta, iphi, theHBHETopology)) ? HcalEndcap : HcalBarrel)
0555 : subdet;
0556 HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
0557 double actL = activeLength(DetId(hcid0), hdc);
0558 double ene = ehdepth[i].first;
0559 double eneAux = ehdepthAux[i].first;
0560 double eneRaw = ehdepthRaw[i].first;
0561 bool tmpC(false);
0562 if (ene > 0.0) {
0563 if (!(theHBHETopology->validHcal(hcid0))) {
0564 edm::LogWarning("HBHEMuon") << "(1) Invalid ID " << hcid0 << " with E = " << ene;
0565 edm::LogWarning("HBHEMuon") << HcalDetId(closestCell) << " with " << ehdepth.size() << " depths:";
0566 for (const auto& ehd : ehdepth)
0567 edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0568 } else {
0569 tmpC = goodCell(hcid0, pTrack, geo, bField, hdc);
0570 double enec(ene);
0571 double corr = respCorr(DetId(hcid0), respCorrs);
0572 if (corr != 0)
0573 ene /= corr;
0574 #ifdef EDM_ML_DEBUG
0575 HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc->mergedDepthDetId(hcid0) : hcid0;
0576 edm::LogVerbatim("HBHEMuon") << hcid0 << ":" << id << " Corr " << corr;
0577 #endif
0578 int depth = ehdepth[i].second - 1;
0579 if (collapseDepth_) {
0580 HcalDetId id = hdc->mergedDepthDetId(hcid0);
0581 depth = id.depth() - 1;
0582 }
0583 eHcalDepth[depth] += ene;
0584 eHcalDepthC[depth] += enec;
0585 activeL[depth] += actL;
0586 activeLengthTot += actL;
0587 matchDepth[depth] = (matchDepth[depth] && tmpC);
0588 #ifdef EDM_ML_DEBUG
0589 if ((verbosity_ % 10) > 0)
0590 edm::LogVerbatim("HBHEMuon")
0591 << hcid0 << " E " << ene << ":" << enec << " L " << actL << " Match " << tmpC;
0592 #endif
0593 }
0594 }
0595 if (eneAux > 0.0) {
0596 if (theHBHETopology->validHcal(hcid0)) {
0597 double enecAux(eneAux);
0598 double corr = respCorr(DetId(hcid0), respCorrs);
0599 if (corr != 0)
0600 eneAux /= corr;
0601 int depth = ehdepthAux[i].second - 1;
0602 if (collapseDepth_) {
0603 HcalDetId id = hdc->mergedDepthDetId(hcid0);
0604 depth = id.depth() - 1;
0605 }
0606 eHcalDepthAux[depth] += eneAux;
0607 eHcalDepthCAux[depth] += enecAux;
0608 #ifdef EDM_ML_DEBUG
0609 if ((verbosity_ % 10) > 0)
0610 edm::LogVerbatim("HBHEMuon")
0611 << hcid0 << " E " << eneAux << ":" << enecAux << " L " << actL << " Match " << tmpC;
0612 #endif
0613 }
0614 }
0615 if (eneRaw > 0.0) {
0616 if (theHBHETopology->validHcal(hcid0)) {
0617 double enecRaw(eneRaw);
0618 double corr = respCorr(DetId(hcid0), respCorrs);
0619 if (corr != 0)
0620 eneRaw /= corr;
0621 int depth = ehdepthRaw[i].second - 1;
0622 if (collapseDepth_) {
0623 HcalDetId id = hdc->mergedDepthDetId(hcid0);
0624 depth = id.depth() - 1;
0625 }
0626 eHcalDepthRaw[depth] += eneRaw;
0627 eHcalDepthCRaw[depth] += enecRaw;
0628 #ifdef EDM_ML_DEBUG
0629 if ((verbosity_ % 10) > 0)
0630 edm::LogVerbatim("HBHEMuon")
0631 << hcid0 << " E " << eneRaw << ":" << enecRaw << " L " << actL << " Match " << tmpC;
0632 #endif
0633 }
0634 }
0635 }
0636 #ifdef EDM_ML_DEBUG
0637 if ((verbosity_ % 10) > 0) {
0638 edm::LogVerbatim("HBHEMuon") << hcidt << " Match " << hbheMuon.matchedId_ << " Depths " << ehdepth.size();
0639 for (unsigned int k = 0; k < ehdepth.size(); ++k)
0640 edm::LogVerbatim("HBHEMuon") << " [" << k << ":" << ehdepth[k].second << "] " << matchDepth[k];
0641 }
0642 #endif
0643 HcalDetId hotCell;
0644 spr::eHCALmatrix(geo, theHBHETopology, closestCell, hbhe, 1, 1, hotCell, false, 0, false);
0645 isHot = matchId(closestCell, hotCell);
0646 if (hotCell != HcalDetId()) {
0647 subdet = HcalDetId(hotCell).subdet();
0648 ieta = HcalDetId(hotCell).ieta();
0649 iphi = HcalDetId(hotCell).iphi();
0650 hborhe = (std::abs(ieta) == 16);
0651 std::vector<std::pair<double, int>> ehdepth, ehdepthAux, ehdepthRaw;
0652 spr::energyHCALCell(hotCell,
0653 hbhe,
0654 ehdepth,
0655 maxDepth_,
0656 -100.0,
0657 -100.0,
0658 -100.0,
0659 -100.0,
0660 -500.0,
0661 500.0,
0662 0,
0663 depth16HE(ieta, iphi, theHBHETopology),
0664 false);
0665 spr::energyHCALCell(hotCell,
0666 hbhe,
0667 ehdepthAux,
0668 maxDepth_,
0669 -100.0,
0670 -100.0,
0671 -100.0,
0672 -100.0,
0673 -500.0,
0674 500.0,
0675 1,
0676 depth16HE(ieta, iphi, theHBHETopology),
0677 false);
0678 spr::energyHCALCell(hotCell,
0679 hbhe,
0680 ehdepthRaw,
0681 maxDepth_,
0682 -100.0,
0683 -100.0,
0684 -100.0,
0685 -100.0,
0686 -500.0,
0687 500.0,
0688 2,
0689 depth16HE(ieta, iphi, theHBHETopology),
0690 false);
0691 for (int i = 0; i < depthMax_; ++i)
0692 eHcalDetId[i] = HcalDetId();
0693 for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0694 HcalSubdetector subdet0 =
0695 (hborhe) ? ((ehdepth[i].second >= depth16HE(ieta, iphi, theHBHETopology)) ? HcalEndcap : HcalBarrel)
0696 : subdet;
0697 HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
0698 double actL = activeLength(DetId(hcid0), hdc);
0699 double ene = ehdepth[i].first;
0700 bool tmpC(false);
0701 if (ene > 0.0) {
0702 if (!(theHBHETopology->validHcal(hcid0))) {
0703 edm::LogWarning("HBHEMuon") << "(2) Invalid ID " << hcid0 << " with E = " << ene;
0704 edm::LogWarning("HBHEMuon") << HcalDetId(hotCell) << " with " << ehdepth.size() << " depths:";
0705 for (const auto& ehd : ehdepth)
0706 edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0707 } else {
0708 tmpC = goodCell(hcid0, pTrack, geo, bField, hdc);
0709 double chg(ene), enec(ene);
0710 double corr = respCorr(DetId(hcid0), respCorrs);
0711 if (corr != 0)
0712 ene /= corr;
0713 #ifdef EDM_ML_DEBUG
0714 HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc->mergedDepthDetId(hcid0) : hcid0;
0715 edm::LogVerbatim("HBHEMuon") << hcid0 << ":" << id << " Corr " << corr << " E " << ene << ":" << enec;
0716 #endif
0717 double gain = gainFactor(conditions, hcid0);
0718 if (gain != 0)
0719 chg /= gain;
0720 #ifdef EDM_ML_DEBUG
0721 edm::LogVerbatim("HBHEMuon") << hcid0 << " Gain " << gain << " C " << chg;
0722 #endif
0723 int depth = ehdepth[i].second - 1;
0724 if (collapseDepth_) {
0725 HcalDetId id = hdc->mergedDepthDetId(hcid0);
0726 depth = id.depth() - 1;
0727 }
0728 eHcalDepthHot[depth] += ene;
0729 eHcalDepthHotC[depth] += enec;
0730 cHcalDepthHot[depth] += chg;
0731 activeHotL[depth] += actL;
0732 activeLengthHotTot += actL;
0733 matchDepthHot[depth] = (matchDepthHot[depth] && tmpC);
0734 #ifdef EDM_ML_DEBUG
0735 eHcal += ene;
0736 if ((verbosity_ % 10) > 0)
0737 edm::LogVerbatim("HBHEMuon") << hcid0 << " depth " << depth << " E " << ene << ":" << enec << " C "
0738 << chg << " L " << actL << " Match " << tmpC;
0739 #endif
0740 }
0741 }
0742 double eneAux = ehdepthAux[i].first;
0743 if (eneAux > 0.0) {
0744 if (theHBHETopology->validHcal(hcid0)) {
0745 double chgAux(eneAux), enecAux(eneAux);
0746 double corr = respCorr(DetId(hcid0), respCorrs);
0747 if (corr != 0)
0748 eneAux /= corr;
0749 double gain = gainFactor(conditions, hcid0);
0750 if (gain != 0)
0751 chgAux /= gain;
0752 int depth = ehdepthAux[i].second - 1;
0753 if (collapseDepth_) {
0754 HcalDetId id = hdc->mergedDepthDetId(hcid0);
0755 depth = id.depth() - 1;
0756 }
0757 eHcalDepthHotAux[depth] += eneAux;
0758 eHcalDepthHotCAux[depth] += enecAux;
0759 cHcalDepthHotAux[depth] += chgAux;
0760 #ifdef EDM_ML_DEBUG
0761 if ((verbosity_ % 10) > 0)
0762 edm::LogVerbatim("HBHEMuon") << hcid0 << " depth " << depth << " E " << eneAux << ":" << enecAux
0763 << " C " << chgAux << " L " << actL << " Match " << tmpC;
0764 #endif
0765 }
0766 }
0767 double eneRaw = ehdepthRaw[i].first;
0768 if (eneRaw > 0.0) {
0769 if (theHBHETopology->validHcal(hcid0)) {
0770 double chgRaw(eneRaw), enecRaw(eneRaw);
0771 double corr = respCorr(DetId(hcid0), respCorrs);
0772 if (corr != 0)
0773 eneRaw /= corr;
0774 double gain = gainFactor(conditions, hcid0);
0775 if (gain != 0)
0776 chgRaw /= gain;
0777 int depth = ehdepthRaw[i].second - 1;
0778 if (collapseDepth_) {
0779 HcalDetId id = hdc->mergedDepthDetId(hcid0);
0780 depth = id.depth() - 1;
0781 }
0782 eHcalDepthHotRaw[depth] += eneRaw;
0783 eHcalDepthHotCRaw[depth] += enecRaw;
0784 cHcalDepthHotRaw[depth] += chgRaw;
0785 #ifdef EDM_ML_DEBUG
0786 if ((verbosity_ % 10) > 0)
0787 edm::LogVerbatim("HBHEMuon") << hcid0 << " depth " << depth << " E " << eneRaw << ":" << enecRaw
0788 << " C " << chgRaw << " L " << actL << " Match " << tmpC;
0789 #endif
0790 }
0791 }
0792 }
0793
0794 HcalDetId oppCell(subdet, -ieta, iphi, HcalDetId(hotCell).depth());
0795 std::vector<std::pair<double, int>> ehdeptho, ehdepthoAux, ehdepthoRaw;
0796 spr::energyHCALCell(oppCell,
0797 hbhe,
0798 ehdeptho,
0799 maxDepth_,
0800 -100.0,
0801 -100.0,
0802 -100.0,
0803 -100.0,
0804 -500.0,
0805 500.0,
0806 0,
0807 depth16HE(-ieta, iphi, theHBHETopology),
0808 false);
0809 spr::energyHCALCell(oppCell,
0810 hbhe,
0811 ehdepthoAux,
0812 maxDepth_,
0813 -100.0,
0814 -100.0,
0815 -100.0,
0816 -100.0,
0817 -500.0,
0818 500.0,
0819 1,
0820 depth16HE(-ieta, iphi, theHBHETopology),
0821 false);
0822 spr::energyHCALCell(oppCell,
0823 hbhe,
0824 ehdepthoRaw,
0825 maxDepth_,
0826 -100.0,
0827 -100.0,
0828 -100.0,
0829 -100.0,
0830 -500.0,
0831 500.0,
0832 2,
0833 depth16HE(-ieta, iphi, theHBHETopology),
0834 false);
0835 for (unsigned int i = 0; i < ehdeptho.size(); ++i) {
0836 HcalSubdetector subdet0 =
0837 (hborhe) ? ((ehdeptho[i].second >= depth16HE(-ieta, iphi, theHBHETopology)) ? HcalEndcap : HcalBarrel)
0838 : subdet;
0839 HcalDetId hcid0(subdet0, -ieta, iphi, ehdeptho[i].second);
0840 double ene = ehdeptho[i].first;
0841 if (ene > 0.0) {
0842 if (!(theHBHETopology->validHcal(hcid0))) {
0843 edm::LogWarning("HBHEMuon") << "(3) Invalid ID " << hcid0 << " with E = " << ene;
0844 edm::LogWarning("HBHEMuon") << oppCell << " with " << ehdeptho.size() << " depths:";
0845 for (const auto& ehd : ehdeptho)
0846 edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0847 } else {
0848 double chg(ene);
0849 double corr = respCorr(DetId(hcid0), respCorrs);
0850 if (corr != 0)
0851 ene /= corr;
0852 #ifdef EDM_ML_DEBUG
0853 HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc->mergedDepthDetId(hcid0) : hcid0;
0854 edm::LogVerbatim("HBHEMuon")
0855 << hcid0 << ":" << id << " Corr " << corr << " E " << ene << ":" << ehdeptho[i].first;
0856 #endif
0857 double gain = gainFactor(conditions, hcid0);
0858 if (gain != 0)
0859 chg /= gain;
0860 #ifdef EDM_ML_DEBUG
0861 edm::LogVerbatim("HBHEMuon") << hcid0 << " Gain " << gain << " C " << chg;
0862 #endif
0863 int depth = ehdeptho[i].second - 1;
0864 if (collapseDepth_) {
0865 HcalDetId id = hdc->mergedDepthDetId(hcid0);
0866 depth = id.depth() - 1;
0867 }
0868 cHcalDepthHotBG[depth] += chg;
0869 #ifdef EDM_ML_DEBUG
0870 if ((verbosity_ % 10) > 0)
0871 edm::LogVerbatim("HBHEMuon") << hcid0 << " Depth " << depth << " E " << ene << " C " << chg;
0872 #endif
0873 }
0874 }
0875 double eneAux = ehdepthoAux[i].first;
0876 if (eneAux > 0.0) {
0877 if (theHBHETopology->validHcal(hcid0)) {
0878 double chgAux(eneAux);
0879 double corr = respCorr(DetId(hcid0), respCorrs);
0880 if (corr != 0)
0881 eneAux /= corr;
0882 double gain = gainFactor(conditions, hcid0);
0883 if (gain != 0)
0884 chgAux /= gain;
0885 int depth = ehdepthoAux[i].second - 1;
0886 if (collapseDepth_) {
0887 HcalDetId id = hdc->mergedDepthDetId(hcid0);
0888 depth = id.depth() - 1;
0889 }
0890 cHcalDepthHotBGAux[depth] += chgAux;
0891 #ifdef EDM_ML_DEBUG
0892 if ((verbosity_ % 10) > 0)
0893 edm::LogVerbatim("HBHEMuon") << hcid0 << " Depth " << depth << " E " << eneAux << " C " << chgAux;
0894 #endif
0895 }
0896 }
0897 double eneRaw = ehdepthoRaw[i].first;
0898 if (eneRaw > 0.0) {
0899 if (theHBHETopology->validHcal(hcid0)) {
0900 double chgRaw(eneRaw);
0901 double corr = respCorr(DetId(hcid0), respCorrs);
0902 if (corr != 0)
0903 eneRaw /= corr;
0904 double gain = gainFactor(conditions, hcid0);
0905 if (gain != 0)
0906 chgRaw /= gain;
0907 int depth = ehdepthoRaw[i].second - 1;
0908 if (collapseDepth_) {
0909 HcalDetId id = hdc->mergedDepthDetId(hcid0);
0910 depth = id.depth() - 1;
0911 }
0912 cHcalDepthHotBGRaw[depth] += chgRaw;
0913 #ifdef EDM_ML_DEBUG
0914 if ((verbosity_ % 10) > 0)
0915 edm::LogVerbatim("HBHEMuon") << hcid0 << " Depth " << depth << " E " << eneRaw << " C " << chgRaw;
0916 #endif
0917 }
0918 }
0919 }
0920 }
0921 }
0922 #ifdef EDM_ML_DEBUG
0923 edm::LogVerbatim("HBHEMuon") << "Propagate Track to HCAL: " << trackID.okHCAL << " Match "
0924 << hbheMuon.matchedId_ << " Hot " << isHot << " Energy " << eHcal;
0925 #endif
0926 hbheMuon.hcalIeta_ = ieta;
0927 hbheMuon.hcalIphi_ = iphi;
0928 for (int i = 0; i < depthMax_; ++i) {
0929 hbheMuon.hcalDepthEnergy_.push_back(eHcalDepth[i]);
0930 hbheMuon.hcalDepthActiveLength_.push_back(activeL[i]);
0931 hbheMuon.hcalDepthEnergyHot_.push_back(eHcalDepthHot[i]);
0932 hbheMuon.hcalDepthActiveLengthHot_.push_back(activeHotL[i]);
0933 hbheMuon.hcalDepthEnergyCorr_.push_back(eHcalDepthC[i]);
0934 hbheMuon.hcalDepthEnergyHotCorr_.push_back(eHcalDepthHotC[i]);
0935 hbheMuon.hcalDepthChargeHot_.push_back(cHcalDepthHot[i]);
0936 hbheMuon.hcalDepthChargeHotBG_.push_back(cHcalDepthHotBG[i]);
0937 hbheMuon.hcalDepthMatch_.push_back(matchDepth[i]);
0938 hbheMuon.hcalDepthMatchHot_.push_back(matchDepthHot[i]);
0939 hbheMuon.hcalDepthEnergyAux_.push_back(eHcalDepthAux[i]);
0940 hbheMuon.hcalDepthEnergyHotAux_.push_back(eHcalDepthHotAux[i]);
0941 hbheMuon.hcalDepthEnergyCorrAux_.push_back(eHcalDepthCAux[i]);
0942 hbheMuon.hcalDepthEnergyHotCorrAux_.push_back(eHcalDepthHotCAux[i]);
0943 hbheMuon.hcalDepthChargeHotAux_.push_back(cHcalDepthHotAux[i]);
0944 hbheMuon.hcalDepthChargeHotBGAux_.push_back(cHcalDepthHotBGAux[i]);
0945 hbheMuon.hcalDepthEnergyRaw_.push_back(eHcalDepthRaw[i]);
0946 hbheMuon.hcalDepthEnergyHotRaw_.push_back(eHcalDepthHotRaw[i]);
0947 hbheMuon.hcalDepthEnergyCorrRaw_.push_back(eHcalDepthCRaw[i]);
0948 hbheMuon.hcalDepthEnergyHotCorrRaw_.push_back(eHcalDepthHotCRaw[i]);
0949 hbheMuon.hcalDepthChargeHotRaw_.push_back(cHcalDepthHotRaw[i]);
0950 hbheMuon.hcalDepthChargeHotBGRaw_.push_back(cHcalDepthHotBGRaw[i]);
0951 }
0952 hbheMuon.hcalActiveLength_ = activeLengthTot;
0953 hbheMuon.hcalHot_ = isHot;
0954 hbheMuon.hcalActiveLengthHot_ = activeLengthHotTot;
0955
0956 if ((recMuon->p() > 10.0) && (trackID.okHCAL))
0957 outputHcalHBHEMuonColl->emplace_back(hbheMuon);
0958 }
0959 }
0960 }
0961 if (!outputHcalHBHEMuonColl->empty())
0962 ++nGood_;
0963 iEvent.put(std::move(outputHcalHBHEMuonColl), labelHBHEMuon_);
0964 }
0965
0966
0967 void AlCaHcalHBHEMuonProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0968 edm::ParameterSetDescription desc;
0969 std::vector<std::string> triggers = {"HLT_IsoMu", "HLT_Mu"};
0970 desc.add<std::vector<std::string>>("triggers", triggers);
0971 desc.add<std::string>("processName", "HLT");
0972 desc.add<edm::InputTag>("triggerResults", edm::InputTag("TriggerResults", "", "HLT"));
0973 desc.add<edm::InputTag>("labelEBRecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0974 desc.add<edm::InputTag>("labelEERecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0975 desc.add<edm::InputTag>("labelHBHERecHit", edm::InputTag("hbhereco"));
0976 desc.add<std::string>("labelVertex", "offlinePrimaryVertices");
0977 desc.add<std::string>("labelMuon", "muons");
0978 desc.add<std::string>("labelHBHEMuon", "hbheMuon");
0979 desc.add<bool>("collapseDepth", false);
0980 desc.add<bool>("isItPlan1", false);
0981 desc.addUntracked<int>("verbosity", 0);
0982 desc.addUntracked<bool>("isItPreRecHit", false);
0983 desc.addUntracked<bool>("writeRespCorr", false);
0984 desc.addUntracked<std::string>("fileInCorr", "");
0985 desc.addUntracked<int>("maxDepth", 7);
0986 desc.add<bool>("usePFThreshold", true);
0987 descriptions.add("alcaHcalHBHEMuonProducer", desc);
0988 }
0989
0990
0991 void AlCaHcalHBHEMuonProducer::endStream() {
0992 globalCache()->nAll_ += nAll_;
0993 globalCache()->nGood_ += nGood_;
0994 }
0995
0996 void AlCaHcalHBHEMuonProducer::globalEndJob(const alcaHcalHBHEMuonProducer::Counters* count) {
0997 edm::LogVerbatim("HBHEMuon") << "Selects " << count->nGood_ << " out of " << count->nAll_ << " total # of events\n";
0998 }
0999
1000
1001 void AlCaHcalHBHEMuonProducer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
1002 const HcalDDDRecConstants* hdc = &iSetup.getData(tok_ddrec0_);
1003 actHB.clear();
1004 actHE.clear();
1005 actHB = hdc->getThickActive(0);
1006 actHE = hdc->getThickActive(1);
1007 #ifdef EDM_ML_DEBUG
1008 unsigned int k1(0), k2(0);
1009 edm::LogVerbatim("HBHEMuon") << actHB.size() << " Active Length for HB";
1010 for (const auto& act : actHB) {
1011 edm::LogVerbatim("HBHEMuon") << "[" << k1 << "] ieta " << act.ieta << " depth " << act.depth << " zside "
1012 << act.zside << " type " << act.stype << " phi " << act.iphis.size() << ":"
1013 << act.iphis[0] << " L " << act.thick;
1014 HcalDetId hcid1(HcalBarrel, (act.ieta) * (act.zside), act.iphis[0], act.depth);
1015 HcalDetId hcid2 = mergedDepth_ ? hdc->mergedDepthDetId(hcid1) : hcid1;
1016 edm::LogVerbatim("HBHEMuon") << hcid1 << " | " << hcid2 << " L " << activeLength(DetId(hcid2), hdc);
1017 ++k1;
1018 }
1019 edm::LogVerbatim("HBHEMuon") << actHE.size() << " Active Length for HE";
1020 for (const auto& act : actHE) {
1021 edm::LogVerbatim("HBHEMuon") << "[" << k2 << "] ieta " << act.ieta << " depth " << act.depth << " zside "
1022 << act.zside << " type " << act.stype << " phi " << act.iphis.size() << ":"
1023 << act.iphis[0] << " L " << act.thick;
1024 HcalDetId hcid1(HcalEndcap, (act.ieta) * (act.zside), act.iphis[0], act.depth);
1025 HcalDetId hcid2 = mergedDepth_ ? hdc->mergedDepthDetId(hcid1) : hcid1;
1026 edm::LogVerbatim("HBHEMuon") << hcid1 << " | " << hcid2 << " L " << activeLength(DetId(hcid2), hdc);
1027 ++k2;
1028 }
1029 #endif
1030
1031 bool changed = true;
1032 bool flag = hltConfig_.init(iRun, iSetup, processName_, changed);
1033 edm::LogVerbatim("HBHEMuon") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init " << flag << std::endl;
1034
1035 const HcalRespCorrs* resp = &iSetup.getData(tok_respcorr0_);
1036 const HcalTopology* theHBHETopology = &iSetup.getData(tok_htopo0_);
1037 const CaloGeometry* geo = &iSetup.getData(tok_geom0_);
1038 HcalRespCorrs respCorrsObj(*resp);
1039 HcalRespCorrs* respCorrs = &respCorrsObj;
1040 respCorrs->setTopo(theHBHETopology);
1041
1042
1043 if (writeRespCorr_) {
1044 const HcalGeometry* gHcal = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
1045 const std::vector<DetId>& ids = gHcal->getValidDetIds(DetId::Hcal, 0);
1046 edm::LogVerbatim("HBHEMuon") << "\nTable of Correction Factors for Run " << iRun.run() << "\n";
1047 for (auto const& id : ids) {
1048 if ((id.det() == DetId::Hcal) && ((id.subdetId() == HcalBarrel) || (id.subdetId() == HcalEndcap))) {
1049 edm::LogVerbatim("HBHEMuon") << HcalDetId(id) << " " << id.rawId() << " "
1050 << (respCorrs->getValues(id))->getValue();
1051 }
1052 }
1053 }
1054 }
1055
1056 void AlCaHcalHBHEMuonProducer::endRun(edm::Run const& iRun, edm::EventSetup const&) {
1057 edm::LogVerbatim("HBHEMuon") << "endRun[" << nRun_ << "] " << iRun.run() << "\n";
1058 ++nRun_;
1059 }
1060
1061
1062
1063 int AlCaHcalHBHEMuonProducer::matchId(const HcalDetId& id1, const HcalDetId& id2) {
1064 HcalDetId kd1(id1.subdet(), id1.ieta(), id1.iphi(), 1);
1065 HcalDetId kd2(id1.subdet(), id2.ieta(), id2.iphi(), 1);
1066 int match = ((kd1 == kd2) ? 1 : 0);
1067 return match;
1068 }
1069
1070 double AlCaHcalHBHEMuonProducer::activeLength(const DetId& id0, const HcalDDDRecConstants* hdc) {
1071 HcalDetId id(id0);
1072 int ieta = id.ietaAbs();
1073 int zside = id.zside();
1074 int iphi = id.iphi();
1075 std::vector<int> dpths;
1076 if (mergedDepth_) {
1077 std::vector<HcalDetId> ids;
1078 hdc->unmergeDepthDetId(id, ids);
1079 for (auto idh : ids)
1080 dpths.emplace_back(idh.depth());
1081 } else {
1082 dpths.emplace_back(id.depth());
1083 }
1084 double lx(0);
1085 if (id.subdet() == HcalBarrel) {
1086 for (unsigned int i = 0; i < actHB.size(); ++i) {
1087 if ((ieta == actHB[i].ieta) && (zside == actHB[i].zside) &&
1088 (std::find(dpths.begin(), dpths.end(), actHB[i].depth) != dpths.end()) &&
1089 (std::find(actHB[i].iphis.begin(), actHB[i].iphis.end(), iphi) != actHB[i].iphis.end())) {
1090 lx += actHB[i].thick;
1091 }
1092 }
1093 } else {
1094 for (unsigned int i = 0; i < actHE.size(); ++i) {
1095 if ((ieta == actHE[i].ieta) && (zside == actHE[i].zside) &&
1096 (std::find(dpths.begin(), dpths.end(), actHE[i].depth) != dpths.end()) &&
1097 (std::find(actHE[i].iphis.begin(), actHE[i].iphis.end(), iphi) != actHE[i].iphis.end())) {
1098 lx += actHE[i].thick;
1099 }
1100 }
1101 }
1102 return lx;
1103 }
1104
1105 bool AlCaHcalHBHEMuonProducer::isGoodVertex(const reco::Vertex& vtx) {
1106 if (vtx.isFake())
1107 return false;
1108 if (vtx.ndof() < 4)
1109 return false;
1110 if (vtx.position().Rho() > 2.)
1111 return false;
1112 if (fabs(vtx.position().Z()) > 24.)
1113 return false;
1114 return true;
1115 }
1116
1117 double AlCaHcalHBHEMuonProducer::respCorr(const DetId& id, const HcalRespCorrs* respCorrs) {
1118 double cfac(1.0);
1119 if (useMyCorr_) {
1120 auto itr = corrValue_.find(id);
1121 if (itr != corrValue_.end())
1122 cfac = itr->second;
1123 } else if (respCorrs != nullptr) {
1124 cfac = (respCorrs->getValues(id))->getValue();
1125 }
1126 return cfac;
1127 }
1128
1129 double AlCaHcalHBHEMuonProducer::gainFactor(const HcalDbService* conditions, const HcalDetId& id) {
1130 double gain(0.0);
1131 const HcalCalibrations& calibs = conditions->getHcalCalibrations(id);
1132 for (int capid = 0; capid < 4; ++capid)
1133 gain += (0.25 * calibs.respcorrgain(capid));
1134 return gain;
1135 }
1136
1137 int AlCaHcalHBHEMuonProducer::depth16HE(int ieta, int iphi, const HcalTopology* theHBHETopology) {
1138
1139
1140
1141 int zside = (ieta > 0) ? 1 : -1;
1142 int depth = theHBHETopology->dddConstants()->getMinDepth(1, 16, iphi, zside);
1143 if (isItPlan1_ && (!isItPreRecHit_))
1144 depth = 3;
1145 #ifdef EDM_ML_DEBUG
1146 edm::LogVerbatim("HBHEMuon") << "Plan1 " << isItPlan1_ << " PreRecHit " << isItPreRecHit_ << " phi " << iphi
1147 << " depth " << depth;
1148 #endif
1149 return depth;
1150 }
1151
1152 bool AlCaHcalHBHEMuonProducer::goodCell(const HcalDetId& hcid,
1153 const reco::Track* pTrack,
1154 const CaloGeometry* geo,
1155 const MagneticField* bField,
1156 const HcalDDDRecConstants* hdc) {
1157 std::pair<double, double> rz = hdc->getRZ(hcid);
1158 bool typeRZ = (hcid.subdet() == HcalEndcap) ? false : true;
1159 bool match = spr::propagateHCAL(pTrack, geo, bField, typeRZ, rz, (((verbosity_ / 10000) % 10) > 0));
1160 return match;
1161 }
1162
1163
1164 #include "FWCore/Framework/interface/MakerMacros.h"
1165
1166 DEFINE_FWK_MODULE(AlCaHcalHBHEMuonProducer);
