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