File indexing completed on 2023-01-25 23:30:06
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/one/EDAnalyzer.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 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0019
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 class HcalHBHEMuonAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0073 public:
0074 explicit HcalHBHEMuonAnalyzer(const edm::ParameterSet&);
0075
0076 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0077
0078 private:
0079 void beginJob() override;
0080 void analyze(edm::Event const&, edm::EventSetup const&) override;
0081 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0082 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0083 void clearVectors();
0084 int matchId(const HcalDetId&, const HcalDetId&);
0085 double activeLength(const DetId&);
0086 bool isGoodVertex(const reco::Vertex& vtx);
0087 double respCorr(const DetId& id);
0088 double gainFactor(const HcalDbService* dbserv, const HcalDetId& id);
0089 int depth16HE(int ieta, int iphi);
0090 bool goodCell(const HcalDetId& hcid, const reco::Track* pTrack, const CaloGeometry* geo, const MagneticField* bField);
0091
0092
0093 HLTConfigProvider hltConfig_;
0094 const edm::InputTag hlTriggerResults_;
0095 const edm::InputTag labelEBRecHit_, labelEERecHit_, labelHBHERecHit_;
0096 const std::string labelVtx_, labelMuon_, fileInCorr_;
0097 const std::vector<std::string> triggers_;
0098 const double pMinMuon_;
0099 const int verbosity_, useRaw_;
0100 const bool unCorrect_, collapseDepth_, isItPlan1_;
0101 const bool ignoreHECorr_, isItPreRecHit_;
0102 const bool getCharge_, writeRespCorr_;
0103 const int maxDepth_;
0104 const std::string modnam_, procnm_;
0105
0106 const edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0107 const edm::EDGetTokenT<reco::VertexCollection> tok_Vtx_;
0108 const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0109 const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0110 const edm::EDGetTokenT<HBHERecHitCollection> tok_HBHE_;
0111 const edm::EDGetTokenT<reco::MuonCollection> tok_Muon_;
0112
0113 const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec_;
0114 const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
0115 const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_respcorr_;
0116 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0117 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0118 const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_chan_;
0119 const edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0120 const edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_topo_;
0121 const edm::ESGetToken<HcalDbService, HcalDbRecord> tok_dbservice_;
0122
0123 const HcalDDDRecConstants* hdc_;
0124 const HcalTopology* theHBHETopology_;
0125 const CaloGeometry* geo_;
0126 HcalRespCorrs* respCorrs_;
0127
0128 bool mergedDepth_, useMyCorr_;
0129 int kount_;
0130
0131
0132 static const int depthMax_ = 7;
0133 TTree* tree_;
0134 unsigned int runNumber_, eventNumber_, lumiNumber_, bxNumber_;
0135 unsigned int goodVertex_;
0136 bool muon_is_good_, muon_global_, muon_tracker_;
0137 bool muon_is_tight_, muon_is_medium_;
0138 double ptGlob_, etaGlob_, phiGlob_, energyMuon_, pMuon_;
0139 float muon_trkKink_, muon_chi2LocalPosition_, muon_segComp_;
0140 int trackerLayer_, numPixelLayers_, tight_PixelHits_;
0141 bool innerTrack_, outerTrack_, globalTrack_;
0142 double chiTracker_, dxyTracker_, dzTracker_;
0143 double innerTrackpt_, innerTracketa_, innerTrackphi_;
0144 double tight_validFraction_, outerTrackChi_;
0145 double outerTrackPt_, outerTrackEta_, outerTrackPhi_;
0146 int outerTrackHits_, outerTrackRHits_;
0147 double globalTrckPt_, globalTrckEta_, globalTrckPhi_;
0148 int globalMuonHits_, matchedStat_;
0149 double chiGlobal_, tight_LongPara_, tight_TransImpara_;
0150 double isolationR04_, isolationR03_;
0151 double ecalEnergy_, hcalEnergy_, hoEnergy_;
0152 bool matchedId_, hcalHot_;
0153 double ecal3x3Energy_, hcal1x1Energy_;
0154 unsigned int ecalDetId_, hcalDetId_, ehcalDetId_;
0155 int hcal_ieta_, hcal_iphi_;
0156 double hcalDepthEnergy_[depthMax_];
0157 double hcalDepthActiveLength_[depthMax_];
0158 double hcalDepthEnergyHot_[depthMax_];
0159 double hcalDepthActiveLengthHot_[depthMax_];
0160 double hcalDepthChargeHot_[depthMax_];
0161 double hcalDepthChargeHotBG_[depthMax_];
0162 double hcalDepthEnergyCorr_[depthMax_];
0163 double hcalDepthEnergyHotCorr_[depthMax_];
0164 bool hcalDepthMatch_[depthMax_];
0165 bool hcalDepthMatchHot_[depthMax_];
0166 double hcalActiveLength_, hcalActiveLengthHot_;
0167 std::vector<std::string> all_triggers_;
0168 std::vector<int> hltresults_;
0169
0170 std::vector<HcalDDDRecConstants::HcalActiveLength> actHB, actHE;
0171 std::map<DetId, double> corrValue_;
0172
0173 };
0174
0175 HcalHBHEMuonAnalyzer::HcalHBHEMuonAnalyzer(const edm::ParameterSet& iConfig)
0176 : hlTriggerResults_(iConfig.getParameter<edm::InputTag>("hlTriggerResults")),
0177 labelEBRecHit_(iConfig.getParameter<edm::InputTag>("labelEBRecHit")),
0178 labelEERecHit_(iConfig.getParameter<edm::InputTag>("labelEERecHit")),
0179 labelHBHERecHit_(iConfig.getParameter<edm::InputTag>("labelHBHERecHit")),
0180 labelVtx_(iConfig.getParameter<std::string>("labelVertex")),
0181 labelMuon_(iConfig.getParameter<std::string>("labelMuon")),
0182 fileInCorr_(iConfig.getUntrackedParameter<std::string>("fileInCorr", "")),
0183 triggers_(iConfig.getParameter<std::vector<std::string>>("triggers")),
0184 pMinMuon_(iConfig.getParameter<double>("pMinMuon")),
0185 verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0186 useRaw_(iConfig.getParameter<int>("useRaw")),
0187 unCorrect_(iConfig.getParameter<bool>("unCorrect")),
0188 collapseDepth_(iConfig.getParameter<bool>("collapseDepth")),
0189 isItPlan1_(iConfig.getParameter<bool>("isItPlan1")),
0190 ignoreHECorr_(iConfig.getUntrackedParameter<bool>("ignoreHECorr", false)),
0191 isItPreRecHit_(iConfig.getUntrackedParameter<bool>("isItPreRecHit", false)),
0192 getCharge_(iConfig.getParameter<bool>("getCharge")),
0193 writeRespCorr_(iConfig.getUntrackedParameter<bool>("writeRespCorr", false)),
0194 maxDepth_(iConfig.getUntrackedParameter<int>("maxDepth", 7)),
0195 modnam_(iConfig.getUntrackedParameter<std::string>("moduleName", "")),
0196 procnm_(iConfig.getUntrackedParameter<std::string>("processName", "")),
0197 tok_trigRes_(consumes<edm::TriggerResults>(hlTriggerResults_)),
0198 tok_Vtx_((modnam_.empty()) ? consumes<reco::VertexCollection>(labelVtx_)
0199 : consumes<reco::VertexCollection>(edm::InputTag(modnam_, labelVtx_, procnm_))),
0200 tok_EB_(consumes<EcalRecHitCollection>(labelEBRecHit_)),
0201 tok_EE_(consumes<EcalRecHitCollection>(labelEERecHit_)),
0202 tok_HBHE_(consumes<HBHERecHitCollection>(labelHBHERecHit_)),
0203 tok_Muon_((modnam_.empty()) ? consumes<reco::MuonCollection>(labelMuon_)
0204 : consumes<reco::MuonCollection>(edm::InputTag(modnam_, labelMuon_, procnm_))),
0205 tok_ddrec_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0206 tok_htopo_(esConsumes<HcalTopology, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0207 tok_respcorr_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd, edm::Transition::BeginRun>()),
0208 tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0209 tok_magField_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0210 tok_chan_(esConsumes<EcalChannelStatus, EcalChannelStatusRcd>()),
0211 tok_sevlv_(esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>()),
0212 tok_topo_(esConsumes<CaloTopology, CaloTopologyRecord>()),
0213 tok_dbservice_(esConsumes<HcalDbService, HcalDbRecord>()),
0214 hdc_(nullptr),
0215 theHBHETopology_(nullptr),
0216 respCorrs_(nullptr) {
0217 usesResource(TFileService::kSharedResource);
0218
0219 kount_ = 0;
0220 mergedDepth_ = (!isItPreRecHit_) || (collapseDepth_);
0221
0222 if (modnam_.empty()) {
0223 edm::LogVerbatim("HBHEMuon") << "Labels used: Trig " << hlTriggerResults_ << " Vtx " << labelVtx_ << " EB "
0224 << labelEBRecHit_ << " EE " << labelEERecHit_ << " HBHE " << labelHBHERecHit_ << " MU "
0225 << labelMuon_;
0226 } else {
0227 edm::LogVerbatim("HBHEMuon") << "Labels used Trig " << hlTriggerResults_ << "\n Vtx "
0228 << edm::InputTag(modnam_, labelVtx_, procnm_) << "\n EB " << labelEBRecHit_
0229 << "\n EE " << labelEERecHit_ << "\n HBHE " << labelHBHERecHit_ << "\n MU "
0230 << edm::InputTag(modnam_, labelMuon_, procnm_);
0231 }
0232
0233 if (!fileInCorr_.empty()) {
0234 std::ifstream infile(fileInCorr_.c_str());
0235 if (infile.is_open()) {
0236 while (true) {
0237 unsigned int id;
0238 double cfac;
0239 infile >> id >> cfac;
0240 if (!infile.good())
0241 break;
0242 corrValue_[DetId(id)] = cfac;
0243 }
0244 infile.close();
0245 }
0246 }
0247 useMyCorr_ = (!corrValue_.empty());
0248 edm::LogVerbatim("HBHEMuon") << "Flags used: UseRaw " << useRaw_ << " GetCharge " << getCharge_ << " UnCorrect "
0249 << unCorrect_ << " IgnoreHECorr " << ignoreHECorr_ << " CollapseDepth " << collapseDepth_
0250 << ":" << mergedDepth_ << " IsItPlan1 " << isItPlan1_ << " IsItPreRecHit "
0251 << isItPreRecHit_ << " UseMyCorr " << useMyCorr_ << " pMinMuon " << pMinMuon_;
0252 }
0253
0254
0255
0256
0257
0258
0259 void HcalHBHEMuonAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0260 ++kount_;
0261 clearVectors();
0262 std::vector<bool> muon_is_good, muon_global, muon_tracker;
0263 std::vector<bool> muon_is_tight, muon_is_medium;
0264 std::vector<double> ptGlob, etaGlob, phiGlob, energyMuon, pMuon;
0265 std::vector<float> muon_trkKink, muon_chi2LocalPosition, muon_segComp;
0266 std::vector<int> trackerLayer, numPixelLayers, tight_PixelHits;
0267 std::vector<bool> innerTrack, outerTrack, globalTrack;
0268 std::vector<double> chiTracker, dxyTracker, dzTracker;
0269 std::vector<double> innerTrackpt, innerTracketa, innerTrackphi;
0270 std::vector<double> tight_validFraction, outerTrackChi;
0271 std::vector<double> outerTrackPt, outerTrackEta, outerTrackPhi;
0272 std::vector<int> outerTrackHits, outerTrackRHits;
0273 std::vector<double> globalTrckPt, globalTrckEta, globalTrckPhi;
0274 std::vector<int> globalMuonHits, matchedStat;
0275 std::vector<double> chiGlobal, tight_LongPara, tight_TransImpara;
0276 std::vector<double> isolationR04, isolationR03;
0277 std::vector<double> ecalEnergy, hcalEnergy, hoEnergy;
0278 std::vector<bool> matchedId, hcalHot;
0279 std::vector<double> ecal3x3Energy, hcal1x1Energy;
0280 std::vector<unsigned int> ecalDetId, hcalDetId, ehcalDetId;
0281 std::vector<int> hcal_ieta, hcal_iphi;
0282 std::vector<double> hcalDepthEnergy[depthMax_];
0283 std::vector<double> hcalDepthActiveLength[depthMax_];
0284 std::vector<double> hcalDepthEnergyHot[depthMax_];
0285 std::vector<double> hcalDepthActiveLengthHot[depthMax_];
0286 std::vector<double> hcalDepthChargeHot[depthMax_];
0287 std::vector<double> hcalDepthChargeHotBG[depthMax_];
0288 std::vector<double> hcalDepthEnergyCorr[depthMax_];
0289 std::vector<double> hcalDepthEnergyHotCorr[depthMax_];
0290 std::vector<bool> hcalDepthMatch[depthMax_];
0291 std::vector<bool> hcalDepthMatchHot[depthMax_];
0292 std::vector<double> hcalActiveLength, hcalActiveLengthHot;
0293 runNumber_ = iEvent.id().run();
0294 eventNumber_ = iEvent.id().event();
0295 lumiNumber_ = iEvent.id().luminosityBlock();
0296 bxNumber_ = iEvent.bunchCrossing();
0297 #ifdef EDM_ML_DEBUG
0298 edm::LogVerbatim("HBHEMuon") << "Run " << runNumber_ << " Event " << eventNumber_ << " Lumi " << lumiNumber_ << " BX "
0299 << bxNumber_ << std::endl;
0300 #endif
0301 const edm::Handle<edm::TriggerResults>& _Triggers = iEvent.getHandle(tok_trigRes_);
0302 #ifdef EDM_ML_DEBUG
0303 if ((verbosity_ / 10000) % 10 > 0)
0304 edm::LogVerbatim("HBHEMuon") << "Size of all triggers " << all_triggers_.size();
0305 #endif
0306 int Ntriggers = static_cast<int>(all_triggers_.size());
0307 #ifdef EDM_ML_DEBUG
0308 if ((verbosity_ / 10000) % 10 > 0)
0309 edm::LogVerbatim("HBHEMuon") << "Size of HLT MENU: " << _Triggers->size();
0310 #endif
0311 if (_Triggers.isValid()) {
0312 const edm::TriggerNames& triggerNames_ = iEvent.triggerNames(*_Triggers);
0313 std::vector<int> index;
0314 for (int i = 0; i < Ntriggers; i++) {
0315 index.push_back(triggerNames_.triggerIndex(all_triggers_[i]));
0316 int triggerSize = static_cast<int>(_Triggers->size());
0317 #ifdef EDM_ML_DEBUG
0318 if ((verbosity_ / 10000) % 10 > 0)
0319 edm::LogVerbatim("HBHEMuon") << "outside loop " << index[i] << "\ntriggerSize " << triggerSize;
0320 #endif
0321 if (index[i] < triggerSize) {
0322 hltresults_.push_back(_Triggers->accept(index[i]));
0323 #ifdef EDM_ML_DEBUG
0324 if ((verbosity_ / 10000) % 10 > 0)
0325 edm::LogVerbatim("HBHEMuon") << "Trigger_info " << triggerSize << " triggerSize " << index[i]
0326 << " trigger_index " << hltresults_.at(i) << " hltresult";
0327 #endif
0328 } else {
0329 if ((verbosity_ / 10000) % 10 > 0)
0330 edm::LogVerbatim("HBHEMuon") << "Requested HLT path \""
0331 << "\" does not exist";
0332 }
0333 }
0334 }
0335
0336
0337 const MagneticField* bField = &iSetup.getData(tok_magField_);
0338 const EcalChannelStatus* theEcalChStatus = &iSetup.getData(tok_chan_);
0339 const EcalSeverityLevelAlgo* sevlv = &iSetup.getData(tok_sevlv_);
0340 const CaloTopology* caloTopology = &iSetup.getData(tok_topo_);
0341 const HcalDbService* conditions = &iSetup.getData(tok_dbservice_);
0342
0343
0344 const edm::Handle<reco::VertexCollection>& vtx = iEvent.getHandle(tok_Vtx_);
0345
0346 edm::Handle<EcalRecHitCollection> barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0347 edm::Handle<EcalRecHitCollection> endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0348
0349 edm::Handle<HBHERecHitCollection> hbhe = iEvent.getHandle(tok_HBHE_);
0350
0351 const edm::Handle<reco::MuonCollection>& _Muon = iEvent.getHandle(tok_Muon_);
0352
0353
0354 math::XYZPoint pvx;
0355 goodVertex_ = 0;
0356 if (!vtx.isValid()) {
0357 #ifdef EDM_ML_DEBUG
0358 edm::LogVerbatim("HBHEMuon") << "No Good Vertex found == Reject";
0359 #endif
0360 return;
0361 }
0362 reco::VertexCollection::const_iterator firstGoodVertex = vtx->end();
0363 for (reco::VertexCollection::const_iterator it = vtx->begin(); it != vtx->end(); it++) {
0364 if (isGoodVertex(*it)) {
0365 if (firstGoodVertex == vtx->end())
0366 firstGoodVertex = it;
0367 ++goodVertex_;
0368 }
0369 }
0370 if (firstGoodVertex != vtx->end())
0371 pvx = firstGoodVertex->position();
0372
0373 bool accept(false);
0374 if (_Muon.isValid() && barrelRecHitsHandle.isValid() && endcapRecHitsHandle.isValid() && hbhe.isValid()) {
0375 for (const auto& RecMuon : (*(_Muon.product()))) {
0376 muon_is_good.push_back(RecMuon.isPFMuon());
0377 muon_global.push_back(RecMuon.isGlobalMuon());
0378 muon_tracker.push_back(RecMuon.isTrackerMuon());
0379 ptGlob.push_back(RecMuon.pt());
0380 etaGlob.push_back(RecMuon.eta());
0381 phiGlob.push_back(RecMuon.phi());
0382 energyMuon.push_back(RecMuon.energy());
0383 pMuon.push_back(RecMuon.p());
0384 #ifdef EDM_ML_DEBUG
0385 edm::LogVerbatim("HBHEMuon") << "Energy:" << RecMuon.energy() << " P:" << RecMuon.p();
0386 #endif
0387 muon_is_tight.push_back(muon::isTightMuon(RecMuon, *firstGoodVertex));
0388 muon_is_medium.push_back(muon::isMediumMuon(RecMuon));
0389 muon_trkKink.push_back(RecMuon.combinedQuality().trkKink);
0390 muon_chi2LocalPosition.push_back(RecMuon.combinedQuality().chi2LocalPosition);
0391 muon_segComp.push_back(muon::segmentCompatibility(RecMuon));
0392
0393 if (RecMuon.track().isNonnull()) {
0394 trackerLayer.push_back(RecMuon.track()->hitPattern().trackerLayersWithMeasurement());
0395 } else {
0396 trackerLayer.push_back(-1);
0397 }
0398 if (RecMuon.innerTrack().isNonnull()) {
0399 innerTrack.push_back(true);
0400 numPixelLayers.push_back(RecMuon.innerTrack()->hitPattern().pixelLayersWithMeasurement());
0401 chiTracker.push_back(RecMuon.innerTrack()->normalizedChi2());
0402 dxyTracker.push_back(fabs(RecMuon.innerTrack()->dxy(pvx)));
0403 dzTracker.push_back(fabs(RecMuon.innerTrack()->dz(pvx)));
0404 innerTrackpt.push_back(RecMuon.innerTrack()->pt());
0405 innerTracketa.push_back(RecMuon.innerTrack()->eta());
0406 innerTrackphi.push_back(RecMuon.innerTrack()->phi());
0407 tight_PixelHits.push_back(RecMuon.innerTrack()->hitPattern().numberOfValidPixelHits());
0408 tight_validFraction.push_back(RecMuon.innerTrack()->validFraction());
0409 } else {
0410 innerTrack.push_back(false);
0411 numPixelLayers.push_back(0);
0412 chiTracker.push_back(0);
0413 dxyTracker.push_back(0);
0414 dzTracker.push_back(0);
0415 innerTrackpt.push_back(0);
0416 innerTracketa.push_back(0);
0417 innerTrackphi.push_back(0);
0418 tight_PixelHits.push_back(0);
0419 tight_validFraction.push_back(-99);
0420 }
0421
0422 if (RecMuon.outerTrack().isNonnull()) {
0423 outerTrack.push_back(true);
0424 outerTrackPt.push_back(RecMuon.outerTrack()->pt());
0425 outerTrackEta.push_back(RecMuon.outerTrack()->eta());
0426 outerTrackPhi.push_back(RecMuon.outerTrack()->phi());
0427 outerTrackChi.push_back(RecMuon.outerTrack()->normalizedChi2());
0428 outerTrackHits.push_back(RecMuon.outerTrack()->numberOfValidHits());
0429 outerTrackRHits.push_back(RecMuon.outerTrack()->recHitsSize());
0430 } else {
0431 outerTrack.push_back(false);
0432 outerTrackPt.push_back(0);
0433 outerTrackEta.push_back(0);
0434 outerTrackPhi.push_back(0);
0435 outerTrackChi.push_back(0);
0436 outerTrackHits.push_back(0);
0437 outerTrackRHits.push_back(0);
0438 }
0439
0440 if (RecMuon.globalTrack().isNonnull()) {
0441 globalTrack.push_back(true);
0442 chiGlobal.push_back(RecMuon.globalTrack()->normalizedChi2());
0443 globalMuonHits.push_back(RecMuon.globalTrack()->hitPattern().numberOfValidMuonHits());
0444 matchedStat.push_back(RecMuon.numberOfMatchedStations());
0445 globalTrckPt.push_back(RecMuon.globalTrack()->pt());
0446 globalTrckEta.push_back(RecMuon.globalTrack()->eta());
0447 globalTrckPhi.push_back(RecMuon.globalTrack()->phi());
0448 tight_TransImpara.push_back(fabs(RecMuon.muonBestTrack()->dxy(pvx)));
0449 tight_LongPara.push_back(fabs(RecMuon.muonBestTrack()->dz(pvx)));
0450 } else {
0451 globalTrack.push_back(false);
0452 chiGlobal.push_back(0);
0453 globalMuonHits.push_back(0);
0454 matchedStat.push_back(0);
0455 globalTrckPt.push_back(0);
0456 globalTrckEta.push_back(0);
0457 globalTrckPhi.push_back(0);
0458 tight_TransImpara.push_back(0);
0459 tight_LongPara.push_back(0);
0460 }
0461
0462 isolationR04.push_back(
0463 ((RecMuon.pfIsolationR04().sumChargedHadronPt +
0464 std::max(0.,
0465 RecMuon.pfIsolationR04().sumNeutralHadronEt + RecMuon.pfIsolationR04().sumPhotonEt -
0466 (0.5 * RecMuon.pfIsolationR04().sumPUPt))) /
0467 RecMuon.pt()));
0468
0469 isolationR03.push_back(
0470 ((RecMuon.pfIsolationR03().sumChargedHadronPt +
0471 std::max(0.,
0472 RecMuon.pfIsolationR03().sumNeutralHadronEt + RecMuon.pfIsolationR03().sumPhotonEt -
0473 (0.5 * RecMuon.pfIsolationR03().sumPUPt))) /
0474 RecMuon.pt()));
0475
0476 ecalEnergy.push_back(RecMuon.calEnergy().emS9);
0477 hcalEnergy.push_back(RecMuon.calEnergy().hadS9);
0478 hoEnergy.push_back(RecMuon.calEnergy().hoS9);
0479
0480 double eEcal(0), eHcal(0), activeLengthTot(0), activeLengthHotTot(0);
0481 double eHcalDepth[depthMax_], eHcalDepthHot[depthMax_];
0482 double eHcalDepthC[depthMax_], eHcalDepthHotC[depthMax_];
0483 double cHcalDepthHot[depthMax_], cHcalDepthHotBG[depthMax_];
0484 double activeL[depthMax_], activeHotL[depthMax_];
0485 bool matchDepth[depthMax_], matchDepthHot[depthMax_];
0486 HcalDetId eHcalDetId[depthMax_];
0487 unsigned int isHot(0);
0488 bool tmpmatch(false);
0489 int ieta(-1000), iphi(-1000);
0490 for (int i = 0; i < depthMax_; ++i) {
0491 eHcalDepth[i] = eHcalDepthHot[i] = 0;
0492 eHcalDepthC[i] = eHcalDepthHotC[i] = 0;
0493 cHcalDepthHot[i] = cHcalDepthHotBG[i] = 0;
0494 activeL[i] = activeHotL[i] = 0;
0495 matchDepth[i] = matchDepthHot[i] = true;
0496 }
0497 if (RecMuon.innerTrack().isNonnull()) {
0498 const reco::Track* pTrack = (RecMuon.innerTrack()).get();
0499 spr::propagatedTrackID trackID = spr::propagateCALO(pTrack, geo_, bField, (((verbosity_ / 100) % 10 > 0)));
0500 if ((RecMuon.p() > pMinMuon_) && (trackID.okHCAL))
0501 accept = true;
0502
0503 ecalDetId.push_back((trackID.detIdECAL)());
0504 hcalDetId.push_back((trackID.detIdHCAL)());
0505 ehcalDetId.push_back((trackID.detIdEHCAL)());
0506
0507 HcalDetId check;
0508 std::pair<bool, HcalDetId> info = spr::propagateHCALBack(pTrack, geo_, bField, (((verbosity_ / 100) % 10 > 0)));
0509 if (info.first) {
0510 check = info.second;
0511 }
0512
0513 bool okE = trackID.okECAL;
0514 if (okE) {
0515 const DetId isoCell(trackID.detIdECAL);
0516 std::pair<double, bool> e3x3 = spr::eECALmatrix(isoCell,
0517 barrelRecHitsHandle,
0518 endcapRecHitsHandle,
0519 *theEcalChStatus,
0520 geo_,
0521 caloTopology,
0522 sevlv,
0523 1,
0524 1,
0525 -100.0,
0526 -100.0,
0527 -500.0,
0528 500.0,
0529 false);
0530 eEcal = e3x3.first;
0531 okE = e3x3.second;
0532 }
0533 #ifdef EDM_ML_DEBUG
0534 edm::LogVerbatim("HBHEMuon") << "Propagate Track to ECAL: " << okE << ":" << trackID.okECAL << " E " << eEcal;
0535 #endif
0536
0537 if (trackID.okHCAL) {
0538 DetId closestCell(trackID.detIdHCAL);
0539 HcalDetId hcidt(closestCell.rawId());
0540 if ((hcidt.ieta() == check.ieta()) && (hcidt.iphi() == check.iphi()))
0541 tmpmatch = true;
0542 #ifdef EDM_ML_DEBUG
0543 edm::LogVerbatim("HBHEMuon") << "Front " << hcidt << " Back " << info.first << ":" << check << " Match "
0544 << tmpmatch;
0545 #endif
0546
0547 HcalSubdetector subdet = hcidt.subdet();
0548 ieta = hcidt.ieta();
0549 iphi = hcidt.iphi();
0550 bool hborhe = (std::abs(ieta) == 16);
0551
0552 eHcal = spr::eHCALmatrix(theHBHETopology_,
0553 closestCell,
0554 hbhe,
0555 0,
0556 0,
0557 false,
0558 true,
0559 -100.0,
0560 -100.0,
0561 -100.0,
0562 -100.0,
0563 -500.,
0564 500.,
0565 useRaw_);
0566 std::vector<std::pair<double, int>> ehdepth;
0567 spr::energyHCALCell((HcalDetId)closestCell,
0568 hbhe,
0569 ehdepth,
0570 maxDepth_,
0571 -100.0,
0572 -100.0,
0573 -100.0,
0574 -100.0,
0575 -500.0,
0576 500.0,
0577 useRaw_,
0578 depth16HE(ieta, iphi),
0579 (((verbosity_ / 1000) % 10) > 0));
0580 for (int i = 0; i < depthMax_; ++i)
0581 eHcalDetId[i] = HcalDetId();
0582 for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0583 HcalSubdetector subdet0 =
0584 (hborhe) ? ((ehdepth[i].second >= depth16HE(ieta, iphi)) ? HcalEndcap : HcalBarrel) : subdet;
0585 HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
0586 double actL = activeLength(DetId(hcid0));
0587 double ene = ehdepth[i].first;
0588 bool tmpC(false);
0589 if (ene > 0.0) {
0590 if (!(theHBHETopology_->validHcal(hcid0))) {
0591 edm::LogWarning("HBHEMuon") << "(1) Invalid ID " << hcid0 << " with E = " << ene;
0592 edm::LogWarning("HBHEMuon") << HcalDetId(closestCell) << " with " << ehdepth.size() << " depths:";
0593 for (const auto& ehd : ehdepth)
0594 edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0595 } else {
0596 tmpC = goodCell(hcid0, pTrack, geo_, bField);
0597 double enec(ene);
0598 if (unCorrect_) {
0599 double corr = (ignoreHECorr_ && (subdet0 == HcalEndcap)) ? 1.0 : respCorr(DetId(hcid0));
0600 if (corr != 0)
0601 ene /= corr;
0602 #ifdef EDM_ML_DEBUG
0603 HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc_->mergedDepthDetId(hcid0) : hcid0;
0604 edm::LogVerbatim("HBHEMuon") << hcid0 << ":" << id << " Corr " << corr;
0605 #endif
0606 }
0607 int depth = ehdepth[i].second - 1;
0608 if (collapseDepth_) {
0609 HcalDetId id = hdc_->mergedDepthDetId(hcid0);
0610 depth = id.depth() - 1;
0611 }
0612 eHcalDepth[depth] += ene;
0613 eHcalDepthC[depth] += enec;
0614 activeL[depth] += actL;
0615 activeLengthTot += actL;
0616 matchDepth[depth] = (matchDepth[depth] && tmpC);
0617 #ifdef EDM_ML_DEBUG
0618 if ((verbosity_ % 10) > 0)
0619 edm::LogVerbatim("HBHEMuon")
0620 << hcid0 << " E " << ene << ":" << enec << " L " << actL << " Match " << tmpC;
0621 #endif
0622 }
0623 }
0624 }
0625 #ifdef EDM_ML_DEBUG
0626 if ((verbosity_ % 10) > 0) {
0627 edm::LogVerbatim("HBHEMuon") << hcidt << " Match " << tmpmatch << " Depths " << ehdepth.size();
0628 for (unsigned int k = 0; k < ehdepth.size(); ++k)
0629 edm::LogVerbatim("HBHEMuon") << " [" << k << ":" << ehdepth[k].second << "] " << matchDepth[k];
0630 }
0631 #endif
0632 HcalDetId hotCell;
0633 spr::eHCALmatrix(geo_, theHBHETopology_, closestCell, hbhe, 1, 1, hotCell, false, useRaw_, false);
0634 isHot = matchId(closestCell, hotCell);
0635 if (hotCell != HcalDetId()) {
0636 subdet = HcalDetId(hotCell).subdet();
0637 ieta = HcalDetId(hotCell).ieta();
0638 iphi = HcalDetId(hotCell).iphi();
0639 hborhe = (std::abs(ieta) == 16);
0640 std::vector<std::pair<double, int>> ehdepth;
0641 spr::energyHCALCell(hotCell,
0642 hbhe,
0643 ehdepth,
0644 maxDepth_,
0645 -100.0,
0646 -100.0,
0647 -100.0,
0648 -100.0,
0649 -500.0,
0650 500.0,
0651 useRaw_,
0652 depth16HE(ieta, iphi),
0653 false);
0654 for (int i = 0; i < depthMax_; ++i)
0655 eHcalDetId[i] = HcalDetId();
0656 for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0657 HcalSubdetector subdet0 =
0658 (hborhe) ? ((ehdepth[i].second >= depth16HE(ieta, iphi)) ? HcalEndcap : HcalBarrel) : subdet;
0659 HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
0660 double actL = activeLength(DetId(hcid0));
0661 double ene = ehdepth[i].first;
0662 bool tmpC(false);
0663 if (ene > 0.0) {
0664 if (!(theHBHETopology_->validHcal(hcid0))) {
0665 edm::LogWarning("HBHEMuon") << "(2) Invalid ID " << hcid0 << " with E = " << ene;
0666 edm::LogWarning("HBHEMuon") << HcalDetId(hotCell) << " with " << ehdepth.size() << " depths:";
0667 for (const auto& ehd : ehdepth)
0668 edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0669 } else {
0670 tmpC = goodCell(hcid0, pTrack, geo_, bField);
0671 double chg(ene), enec(ene);
0672 if (unCorrect_) {
0673 double corr = (ignoreHECorr_ && (subdet0 == HcalEndcap)) ? 1.0 : respCorr(DetId(hcid0));
0674 if (corr != 0)
0675 ene /= corr;
0676 #ifdef EDM_ML_DEBUG
0677 HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc_->mergedDepthDetId(hcid0) : hcid0;
0678 edm::LogVerbatim("HBHEMuon")
0679 << hcid0 << ":" << id << " Corr " << corr << " E " << ene << ":" << enec;
0680 #endif
0681 }
0682 if (getCharge_) {
0683 double gain = gainFactor(conditions, hcid0);
0684 if (gain != 0)
0685 chg /= gain;
0686 #ifdef EDM_ML_DEBUG
0687 edm::LogVerbatim("HBHEMuon") << hcid0 << " Gain " << gain << " C " << chg;
0688 #endif
0689 }
0690 int depth = ehdepth[i].second - 1;
0691 if (collapseDepth_) {
0692 HcalDetId id = hdc_->mergedDepthDetId(hcid0);
0693 depth = id.depth() - 1;
0694 }
0695 eHcalDepthHot[depth] += ene;
0696 eHcalDepthHotC[depth] += enec;
0697 cHcalDepthHot[depth] += chg;
0698 activeHotL[depth] += actL;
0699 activeLengthHotTot += actL;
0700 matchDepthHot[depth] = (matchDepthHot[depth] && tmpC);
0701 #ifdef EDM_ML_DEBUG
0702 if ((verbosity_ % 10) > 0)
0703 edm::LogVerbatim("HBHEMuon") << hcid0 << " depth " << depth << " E " << ene << ":" << enec << " C "
0704 << chg << " L " << actL << " Match " << tmpC;
0705 #endif
0706 }
0707 }
0708 }
0709
0710 HcalDetId oppCell(subdet, -ieta, iphi, HcalDetId(hotCell).depth());
0711 std::vector<std::pair<double, int>> ehdeptho;
0712 spr::energyHCALCell(oppCell,
0713 hbhe,
0714 ehdeptho,
0715 maxDepth_,
0716 -100.0,
0717 -100.0,
0718 -100.0,
0719 -100.0,
0720 -500.0,
0721 500.0,
0722 useRaw_,
0723 depth16HE(-ieta, iphi),
0724 false);
0725 for (unsigned int i = 0; i < ehdeptho.size(); ++i) {
0726 HcalSubdetector subdet0 =
0727 (hborhe) ? ((ehdeptho[i].second >= depth16HE(-ieta, iphi)) ? HcalEndcap : HcalBarrel) : subdet;
0728 HcalDetId hcid0(subdet0, -ieta, iphi, ehdeptho[i].second);
0729 double ene = ehdeptho[i].first;
0730 if (ene > 0.0) {
0731 if (!(theHBHETopology_->validHcal(hcid0))) {
0732 edm::LogWarning("HBHEMuon") << "(3) Invalid ID " << hcid0 << " with E = " << ene;
0733 edm::LogWarning("HBHEMuon") << oppCell << " with " << ehdeptho.size() << " depths:";
0734 for (const auto& ehd : ehdeptho)
0735 edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0736 } else {
0737 double chg(ene);
0738 if (unCorrect_) {
0739 double corr = (ignoreHECorr_ && (subdet0 == HcalEndcap)) ? 1.0 : respCorr(DetId(hcid0));
0740 if (corr != 0)
0741 ene /= corr;
0742 #ifdef EDM_ML_DEBUG
0743 HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc_->mergedDepthDetId(hcid0) : hcid0;
0744 edm::LogVerbatim("HBHEMuon")
0745 << hcid0 << ":" << id << " Corr " << corr << " E " << ene << ":" << ehdeptho[i].first;
0746 #endif
0747 }
0748 if (getCharge_) {
0749 double gain = gainFactor(conditions, hcid0);
0750 if (gain != 0)
0751 chg /= gain;
0752 #ifdef EDM_ML_DEBUG
0753 edm::LogVerbatim("HBHEMuon") << hcid0 << " Gain " << gain << " C " << chg;
0754 #endif
0755 }
0756 int depth = ehdeptho[i].second - 1;
0757 if (collapseDepth_) {
0758 HcalDetId id = hdc_->mergedDepthDetId(hcid0);
0759 depth = id.depth() - 1;
0760 }
0761 cHcalDepthHotBG[depth] += chg;
0762 #ifdef EDM_ML_DEBUG
0763 if ((verbosity_ % 10) > 0)
0764 edm::LogVerbatim("HBHEMuon") << hcid0 << " Depth " << depth << " E " << ene << " C " << chg;
0765 #endif
0766 }
0767 }
0768 }
0769 }
0770 }
0771 #ifdef EDM_ML_DEBUG
0772 edm::LogVerbatim("HBHEMuon") << "Propagate Track to HCAL: " << trackID.okHCAL << " Match " << tmpmatch
0773 << " Hot " << isHot << " Energy " << eHcal << std::endl;
0774 #endif
0775
0776 } else {
0777 ecalDetId.push_back(0);
0778 hcalDetId.push_back(0);
0779 ehcalDetId.push_back(0);
0780 }
0781
0782 matchedId.push_back(tmpmatch);
0783 ecal3x3Energy.push_back(eEcal);
0784 hcal1x1Energy.push_back(eHcal);
0785 hcal_ieta.push_back(ieta);
0786 hcal_iphi.push_back(iphi);
0787 for (int i = 0; i < depthMax_; ++i) {
0788 hcalDepthEnergy[i].push_back(eHcalDepth[i]);
0789 hcalDepthActiveLength[i].push_back(activeL[i]);
0790 hcalDepthEnergyHot[i].push_back(eHcalDepthHot[i]);
0791 hcalDepthActiveLengthHot[i].push_back(activeHotL[i]);
0792 hcalDepthEnergyCorr[i].push_back(eHcalDepthC[i]);
0793 hcalDepthEnergyHotCorr[i].push_back(eHcalDepthHotC[i]);
0794 hcalDepthChargeHot[i].push_back(cHcalDepthHot[i]);
0795 hcalDepthChargeHotBG[i].push_back(cHcalDepthHotBG[i]);
0796 hcalDepthMatch[i].push_back(matchDepth[i]);
0797 hcalDepthMatchHot[i].push_back(matchDepthHot[i]);
0798 }
0799 hcalActiveLength.push_back(activeLengthTot);
0800 hcalHot.push_back(isHot);
0801 hcalActiveLengthHot.push_back(activeLengthHotTot);
0802 }
0803 }
0804 if (accept) {
0805 #ifdef EDM_ML_DEBUG
0806 for (unsigned int i = 0; i < hcal_ieta.size(); ++i)
0807 edm::LogVerbatim("HBHEMuon") << "[" << i << "] ieta/iphi for entry to "
0808 << "HCAL has value of " << hcal_ieta[i] << ":" << hcal_iphi[i];
0809 #endif
0810 for (unsigned int k = 0; k < muon_is_good.size(); ++k) {
0811 muon_is_good_ = muon_is_good[k];
0812 muon_global_ = muon_global[k];
0813 muon_tracker_ = muon_tracker[k];
0814 muon_is_tight_ = muon_is_tight[k];
0815 muon_is_medium_ = muon_is_medium[k];
0816 ptGlob_ = ptGlob[k];
0817 etaGlob_ = etaGlob[k];
0818 phiGlob_ = phiGlob[k];
0819 energyMuon_ = energyMuon[k];
0820 pMuon_ = pMuon[k];
0821 muon_trkKink_ = muon_trkKink[k];
0822 muon_chi2LocalPosition_ = muon_chi2LocalPosition[k];
0823 muon_segComp_ = muon_segComp[k];
0824 trackerLayer_ = trackerLayer[k];
0825 numPixelLayers_ = numPixelLayers[k];
0826 tight_PixelHits_ = tight_PixelHits[k];
0827 innerTrack_ = innerTrack[k];
0828 outerTrack_ = outerTrack[k];
0829 globalTrack_ = globalTrack[k];
0830 chiTracker_ = chiTracker[k];
0831 dxyTracker_ = dxyTracker[k];
0832 dzTracker_ = dzTracker[k];
0833 innerTrackpt_ = innerTrackpt[k];
0834 innerTracketa_ = innerTracketa[k];
0835 innerTrackphi_ = innerTrackphi[k];
0836 tight_validFraction_ = tight_validFraction[k];
0837 outerTrackChi_ = outerTrackChi[k];
0838 outerTrackPt_ = outerTrackPt[k];
0839 outerTrackEta_ = outerTrackEta[k];
0840 outerTrackPhi_ = outerTrackPhi[k];
0841 outerTrackHits_ = outerTrackHits[k];
0842 outerTrackRHits_ = outerTrackRHits[k];
0843 globalTrckPt_ = globalTrckPt[k];
0844 globalTrckEta_ = globalTrckEta[k];
0845 globalTrckPhi_ = globalTrckPhi[k];
0846 globalMuonHits_ = globalMuonHits[k];
0847 matchedStat_ = matchedStat[k];
0848 chiGlobal_ = chiGlobal[k];
0849 tight_LongPara_ = tight_LongPara[k];
0850 tight_TransImpara_ = tight_TransImpara[k];
0851 isolationR04_ = isolationR04[k];
0852 isolationR03_ = isolationR03[k];
0853 ecalEnergy_ = ecalEnergy[k];
0854 hcalEnergy_ = hcalEnergy[k];
0855 hoEnergy_ = hoEnergy[k];
0856 matchedId_ = matchedId[k];
0857 hcalHot_ = hcalHot[k];
0858 ecal3x3Energy_ = ecal3x3Energy[k];
0859 hcal1x1Energy_ = hcal1x1Energy[k];
0860 ecalDetId_ = ecalDetId[k];
0861 hcalDetId_ = hcalDetId[k];
0862 ehcalDetId_ = ehcalDetId[k];
0863 hcal_ieta_ = hcal_ieta[k];
0864 hcal_iphi_ = hcal_iphi[k];
0865 for (int i = 0; i < depthMax_; ++i) {
0866 hcalDepthEnergy_[i] = hcalDepthEnergy[i][k];
0867 hcalDepthActiveLength_[i] = hcalDepthActiveLength[i][k];
0868 hcalDepthEnergyHot_[i] = hcalDepthEnergyHot[i][k];
0869 hcalDepthActiveLengthHot_[i] = hcalDepthActiveLengthHot[i][k];
0870 hcalDepthChargeHot_[i] = hcalDepthChargeHot[i][k];
0871 hcalDepthChargeHotBG_[i] = hcalDepthChargeHotBG[i][k];
0872 hcalDepthEnergyCorr_[i] = hcalDepthEnergyCorr[i][k];
0873 hcalDepthEnergyHotCorr_[i] = hcalDepthEnergyHotCorr[i][k];
0874 hcalDepthMatch_[i] = hcalDepthMatch[i][k];
0875 hcalDepthMatchHot_[i] = hcalDepthMatchHot[i][k];
0876 }
0877 hcalActiveLength_ = hcalActiveLength[k];
0878 hcalActiveLengthHot_ = hcalActiveLengthHot[k];
0879 tree_->Fill();
0880 }
0881 }
0882 }
0883
0884
0885 void HcalHBHEMuonAnalyzer::beginJob() {
0886 edm::Service<TFileService> fs;
0887 tree_ = fs->make<TTree>("TREE", "TREE");
0888 tree_->Branch("Event_No", &eventNumber_);
0889 tree_->Branch("Run_No", &runNumber_);
0890 tree_->Branch("LumiNumber", &lumiNumber_);
0891 tree_->Branch("BXNumber", &bxNumber_);
0892 tree_->Branch("GoodVertex", &goodVertex_);
0893 tree_->Branch("PF_Muon", &muon_is_good_);
0894 tree_->Branch("Global_Muon", &muon_global_);
0895 tree_->Branch("Tracker_muon", &muon_tracker_);
0896 tree_->Branch("MuonIsTight", &muon_is_tight_);
0897 tree_->Branch("MuonIsMedium", &muon_is_medium_);
0898 tree_->Branch("pt_of_muon", &ptGlob_);
0899 tree_->Branch("eta_of_muon", &etaGlob_);
0900 tree_->Branch("phi_of_muon", &phiGlob_);
0901 tree_->Branch("energy_of_muon", &energyMuon_);
0902 tree_->Branch("p_of_muon", &pMuon_);
0903 tree_->Branch("muon_trkKink", &muon_trkKink_);
0904 tree_->Branch("muon_chi2LocalPosition", &muon_chi2LocalPosition_);
0905 tree_->Branch("muon_segComp", &muon_segComp_);
0906
0907 tree_->Branch("TrackerLayer", &trackerLayer_);
0908 tree_->Branch("NumPixelLayers", &numPixelLayers_);
0909 tree_->Branch("InnerTrackPixelHits", &tight_PixelHits_);
0910 tree_->Branch("innerTrack", &innerTrack_);
0911 tree_->Branch("chiTracker", &chiTracker_);
0912 tree_->Branch("DxyTracker", &dxyTracker_);
0913 tree_->Branch("DzTracker", &dzTracker_);
0914 tree_->Branch("innerTrackpt", &innerTrackpt_);
0915 tree_->Branch("innerTracketa", &innerTracketa_);
0916 tree_->Branch("innerTrackphi", &innerTrackphi_);
0917 tree_->Branch("tight_validFraction", &tight_validFraction_);
0918
0919 tree_->Branch("OuterTrack", &outerTrack_);
0920 tree_->Branch("OuterTrackChi", &outerTrackChi_);
0921 tree_->Branch("OuterTrackPt", &outerTrackPt_);
0922 tree_->Branch("OuterTrackEta", &outerTrackEta_);
0923 tree_->Branch("OuterTrackPhi", &outerTrackPhi_);
0924 tree_->Branch("OuterTrackHits", &outerTrackHits_);
0925 tree_->Branch("OuterTrackRHits", &outerTrackRHits_);
0926
0927 tree_->Branch("GlobalTrack", &globalTrack_);
0928 tree_->Branch("GlobalTrckPt", &globalTrckPt_);
0929 tree_->Branch("GlobalTrckEta", &globalTrckEta_);
0930 tree_->Branch("GlobalTrckPhi", &globalTrckPhi_);
0931 tree_->Branch("Global_Muon_Hits", &globalMuonHits_);
0932 tree_->Branch("MatchedStations", &matchedStat_);
0933 tree_->Branch("GlobTrack_Chi", &chiGlobal_);
0934 tree_->Branch("Tight_LongitudinalImpactparameter", &tight_LongPara_);
0935 tree_->Branch("Tight_TransImpactparameter", &tight_TransImpara_);
0936
0937 tree_->Branch("IsolationR04", &isolationR04_);
0938 tree_->Branch("IsolationR03", &isolationR03_);
0939 tree_->Branch("ecal_3into3", &ecalEnergy_);
0940 tree_->Branch("hcal_3into3", &hcalEnergy_);
0941 tree_->Branch("tracker_3into3", &hoEnergy_);
0942
0943 tree_->Branch("matchedId", &matchedId_);
0944 tree_->Branch("hcal_cellHot", &hcalHot_);
0945
0946 tree_->Branch("ecal_3x3", &ecal3x3Energy_);
0947 tree_->Branch("hcal_1x1", &hcal1x1Energy_);
0948 tree_->Branch("ecal_detID", &ecalDetId_);
0949 tree_->Branch("hcal_detID", &hcalDetId_);
0950 tree_->Branch("ehcal_detID", &ehcalDetId_);
0951 tree_->Branch("hcal_ieta", &hcal_ieta_);
0952 tree_->Branch("hcal_iphi", &hcal_iphi_);
0953
0954 char name[100];
0955 for (int k = 0; k < maxDepth_; ++k) {
0956 sprintf(name, "hcal_edepth%d", (k + 1));
0957 tree_->Branch(name, &hcalDepthEnergy_[k]);
0958 sprintf(name, "hcal_activeL%d", (k + 1));
0959 tree_->Branch(name, &hcalDepthActiveLength_[k]);
0960 sprintf(name, "hcal_edepthHot%d", (k + 1));
0961 tree_->Branch(name, &hcalDepthEnergyHot_[k]);
0962 sprintf(name, "hcal_activeHotL%d", (k + 1));
0963 tree_->Branch(name, &hcalDepthActiveLengthHot_[k]);
0964 sprintf(name, "hcal_cdepthHot%d", (k + 1));
0965 tree_->Branch(name, &hcalDepthChargeHot_[k]);
0966 sprintf(name, "hcal_cdepthHotBG%d", (k + 1));
0967 tree_->Branch(name, &hcalDepthChargeHotBG_[k]);
0968 sprintf(name, "hcal_edepthCorrect%d", (k + 1));
0969 tree_->Branch(name, &hcalDepthEnergyCorr_[k]);
0970 sprintf(name, "hcal_edepthHotCorrect%d", (k + 1));
0971 tree_->Branch(name, &hcalDepthEnergyHotCorr_[k]);
0972 sprintf(name, "hcal_depthMatch%d", (k + 1));
0973 tree_->Branch(name, &hcalDepthMatch_[k]);
0974 sprintf(name, "hcal_depthMatchHot%d", (k + 1));
0975 tree_->Branch(name, &hcalDepthMatchHot_[k]);
0976 }
0977
0978 tree_->Branch("activeLength", &hcalActiveLength_);
0979 tree_->Branch("activeLengthHot", &hcalActiveLengthHot_);
0980
0981 tree_->Branch("hltresults", &hltresults_);
0982 tree_->Branch("all_triggers", &all_triggers_);
0983 }
0984
0985
0986 void HcalHBHEMuonAnalyzer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0987 hdc_ = &iSetup.getData(tok_ddrec_);
0988 actHB.clear();
0989 actHE.clear();
0990 actHB = hdc_->getThickActive(0);
0991 actHE = hdc_->getThickActive(1);
0992 #ifdef EDM_ML_DEBUG
0993 unsigned int k1(0), k2(0);
0994 edm::LogVerbatim("HBHEMuon") << actHB.size() << " Active Length for HB";
0995 for (const auto& act : actHB) {
0996 edm::LogVerbatim("HBHEMuon") << "[" << k1 << "] ieta " << act.ieta << " depth " << act.depth << " zside "
0997 << act.zside << " type " << act.stype << " phi " << act.iphis.size() << ":"
0998 << act.iphis[0] << " L " << act.thick;
0999 HcalDetId hcid1(HcalBarrel, (act.ieta) * (act.zside), act.iphis[0], act.depth);
1000 HcalDetId hcid2 = mergedDepth_ ? hdc_->mergedDepthDetId(hcid1) : hcid1;
1001 edm::LogVerbatim("HBHEMuon") << hcid1 << " | " << hcid2 << " L " << activeLength(DetId(hcid2));
1002 ++k1;
1003 }
1004 edm::LogVerbatim("HBHEMuon") << actHE.size() << " Active Length for HE";
1005 for (const auto& act : actHE) {
1006 edm::LogVerbatim("HBHEMuon") << "[" << k2 << "] ieta " << act.ieta << " depth " << act.depth << " zside "
1007 << act.zside << " type " << act.stype << " phi " << act.iphis.size() << ":"
1008 << act.iphis[0] << " L " << act.thick;
1009 HcalDetId hcid1(HcalEndcap, (act.ieta) * (act.zside), act.iphis[0], act.depth);
1010 HcalDetId hcid2 = mergedDepth_ ? hdc_->mergedDepthDetId(hcid1) : hcid1;
1011 edm::LogVerbatim("HBHEMuon") << hcid1 << " | " << hcid2 << " L " << activeLength(DetId(hcid2));
1012 ++k2;
1013 }
1014 #endif
1015
1016 bool changed = true;
1017 all_triggers_.clear();
1018 if (hltConfig_.init(iRun, iSetup, "HLT", changed)) {
1019
1020 #ifdef EDM_ML_DEBUG
1021 edm::LogVerbatim("HBHEMuon") << "HLT config with process name HLT successfully extracted";
1022 #endif
1023 unsigned int ntriggers = hltConfig_.size();
1024 for (unsigned int t = 0; t < ntriggers; ++t) {
1025 std::string hltname(hltConfig_.triggerName(t));
1026 for (unsigned int ik = 0; ik < triggers_.size(); ++ik) {
1027 if (hltname.find(triggers_[ik]) != std::string::npos) {
1028 all_triggers_.push_back(hltname);
1029 break;
1030 }
1031 }
1032 }
1033 edm::LogVerbatim("HBHEMuon") << "All triggers size in begin run " << all_triggers_.size();
1034 } else {
1035 edm::LogError("HBHEMuon") << "Error! HLT config extraction with process name HLT failed";
1036 }
1037
1038 theHBHETopology_ = &iSetup.getData(tok_htopo_);
1039 const HcalRespCorrs* resp = &iSetup.getData(tok_respcorr_);
1040 respCorrs_ = new HcalRespCorrs(*resp);
1041 respCorrs_->setTopo(theHBHETopology_);
1042 geo_ = &iSetup.getData(tok_geom_);
1043
1044
1045 if (writeRespCorr_) {
1046 const HcalGeometry* gHcal = static_cast<const HcalGeometry*>(geo_->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
1047 const std::vector<DetId>& ids = gHcal->getValidDetIds(DetId::Hcal, 0);
1048 edm::LogVerbatim("HBHEMuon") << "\nTable of Correction Factors for Run " << iRun.run() << "\n";
1049 for (auto const& id : ids) {
1050 if ((id.det() == DetId::Hcal) && ((id.subdetId() == HcalBarrel) || (id.subdetId() == HcalEndcap))) {
1051 edm::LogVerbatim("HBHEMuon") << HcalDetId(id) << " " << id.rawId() << " "
1052 << (respCorrs_->getValues(id))->getValue();
1053 }
1054 }
1055 }
1056 }
1057
1058
1059 void HcalHBHEMuonAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1060 edm::ParameterSetDescription desc;
1061 desc.add<edm::InputTag>("hlTriggerResults", edm::InputTag("TriggerResults", "", "HLT"));
1062 desc.add<edm::InputTag>("labelEBRecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
1063 desc.add<edm::InputTag>("labelEERecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
1064 desc.add<edm::InputTag>("labelHBHERecHit", edm::InputTag("hbhereco"));
1065 desc.add<std::string>("labelVertex", "offlinePrimaryVertices");
1066 desc.add<std::string>("labelMuon", "muons");
1067 std::vector<std::string> trig = {};
1068 desc.add<std::vector<std::string>>("triggers", trig);
1069 desc.add<double>("pMinMuon", 5.0);
1070 desc.addUntracked<int>("verbosity", 0);
1071 desc.add<int>("useRaw", 0);
1072 desc.add<bool>("unCorrect", true);
1073 desc.add<bool>("getCharge", true);
1074 desc.add<bool>("collapseDepth", false);
1075 desc.add<bool>("isItPlan1", false);
1076 desc.addUntracked<bool>("ignoreHECorr", false);
1077 desc.addUntracked<bool>("isItPreRecHit", false);
1078 desc.addUntracked<std::string>("moduleName", "");
1079 desc.addUntracked<std::string>("processName", "");
1080 desc.addUntracked<int>("maxDepth", 7);
1081 desc.addUntracked<std::string>("fileInCorr", "");
1082 desc.addUntracked<bool>("writeRespCorr", false);
1083 descriptions.add("hcalHBHEMuon", desc);
1084 }
1085
1086 void HcalHBHEMuonAnalyzer::clearVectors() {
1087
1088 eventNumber_ = -99999;
1089 runNumber_ = -99999;
1090 lumiNumber_ = -99999;
1091 bxNumber_ = -99999;
1092 goodVertex_ = -99999;
1093
1094 muon_is_good_ = false;
1095 muon_global_ = false;
1096 muon_tracker_ = false;
1097 ptGlob_ = 0;
1098 etaGlob_ = 0;
1099 phiGlob_ = 0;
1100 energyMuon_ = 0;
1101 pMuon_ = 0;
1102 muon_trkKink_ = 0;
1103 muon_chi2LocalPosition_ = 0;
1104 muon_segComp_ = 0;
1105 muon_is_tight_ = false;
1106 muon_is_medium_ = false;
1107
1108 trackerLayer_ = 0;
1109 numPixelLayers_ = 0;
1110 tight_PixelHits_ = 0;
1111 innerTrack_ = false;
1112 chiTracker_ = 0;
1113 dxyTracker_ = 0;
1114 dzTracker_ = 0;
1115 innerTrackpt_ = 0;
1116 innerTracketa_ = 0;
1117 innerTrackphi_ = 0;
1118 tight_validFraction_ = 0;
1119
1120 outerTrack_ = false;
1121 outerTrackPt_ = 0;
1122 outerTrackEta_ = 0;
1123 outerTrackPhi_ = 0;
1124 outerTrackHits_ = 0;
1125 outerTrackRHits_ = 0;
1126 outerTrackChi_ = 0;
1127
1128 globalTrack_ = false;
1129 globalTrckPt_ = 0;
1130 globalTrckEta_ = 0;
1131 globalTrckPhi_ = 0;
1132 globalMuonHits_ = 0;
1133 matchedStat_ = 0;
1134 chiGlobal_ = 0;
1135 tight_LongPara_ = 0;
1136 tight_TransImpara_ = 0;
1137
1138 isolationR04_ = 0;
1139 isolationR03_ = 0;
1140 ecalEnergy_ = 0;
1141 hcalEnergy_ = 0;
1142 hoEnergy_ = 0;
1143 matchedId_ = false;
1144 hcalHot_ = false;
1145 ecal3x3Energy_ = 0;
1146 hcal1x1Energy_ = 0;
1147 ecalDetId_ = 0;
1148 hcalDetId_ = 0;
1149 ehcalDetId_ = 0;
1150 hcal_ieta_ = 0;
1151 hcal_iphi_ = 0;
1152 for (int i = 0; i < maxDepth_; ++i) {
1153 hcalDepthEnergy_[i] = 0;
1154 hcalDepthActiveLength_[i] = 0;
1155 hcalDepthEnergyHot_[i] = 0;
1156 hcalDepthActiveLengthHot_[i] = 0;
1157 hcalDepthChargeHot_[i] = 0;
1158 hcalDepthChargeHotBG_[i] = 0;
1159 hcalDepthEnergyCorr_[i] = 0;
1160 hcalDepthEnergyHotCorr_[i] = 0;
1161 hcalDepthMatch_[i] = false;
1162 hcalDepthMatchHot_[i] = false;
1163 }
1164 hcalActiveLength_ = 0;
1165 hcalActiveLengthHot_ = 0;
1166 hltresults_.clear();
1167 }
1168
1169 int HcalHBHEMuonAnalyzer::matchId(const HcalDetId& id1, const HcalDetId& id2) {
1170 HcalDetId kd1(id1.subdet(), id1.ieta(), id1.iphi(), 1);
1171 HcalDetId kd2(id1.subdet(), id2.ieta(), id2.iphi(), 1);
1172 int match = ((kd1 == kd2) ? 1 : 0);
1173 return match;
1174 }
1175
1176 double HcalHBHEMuonAnalyzer::activeLength(const DetId& id_) {
1177 HcalDetId id(id_);
1178 int ieta = id.ietaAbs();
1179 int zside = id.zside();
1180 int iphi = id.iphi();
1181 std::vector<int> dpths;
1182 if (mergedDepth_) {
1183 std::vector<HcalDetId> ids;
1184 hdc_->unmergeDepthDetId(id, ids);
1185 for (auto idh : ids)
1186 dpths.emplace_back(idh.depth());
1187 } else {
1188 dpths.emplace_back(id.depth());
1189 }
1190 double lx(0);
1191 if (id.subdet() == HcalBarrel) {
1192 for (unsigned int i = 0; i < actHB.size(); ++i) {
1193 if ((ieta == actHB[i].ieta) && (zside == actHB[i].zside) &&
1194 (std::find(dpths.begin(), dpths.end(), actHB[i].depth) != dpths.end()) &&
1195 (std::find(actHB[i].iphis.begin(), actHB[i].iphis.end(), iphi) != actHB[i].iphis.end())) {
1196 lx += actHB[i].thick;
1197 }
1198 }
1199 } else {
1200 for (unsigned int i = 0; i < actHE.size(); ++i) {
1201 if ((ieta == actHE[i].ieta) && (zside == actHE[i].zside) &&
1202 (std::find(dpths.begin(), dpths.end(), actHE[i].depth) != dpths.end()) &&
1203 (std::find(actHE[i].iphis.begin(), actHE[i].iphis.end(), iphi) != actHE[i].iphis.end())) {
1204 lx += actHE[i].thick;
1205 }
1206 }
1207 }
1208 return lx;
1209 }
1210
1211 bool HcalHBHEMuonAnalyzer::isGoodVertex(const reco::Vertex& vtx) {
1212 if (vtx.isFake())
1213 return false;
1214 if (vtx.ndof() < 4)
1215 return false;
1216 if (vtx.position().Rho() > 2.)
1217 return false;
1218 if (fabs(vtx.position().Z()) > 24.)
1219 return false;
1220 return true;
1221 }
1222
1223 double HcalHBHEMuonAnalyzer::respCorr(const DetId& id) {
1224 double cfac(1.0);
1225 if (useMyCorr_) {
1226 auto itr = corrValue_.find(id);
1227 if (itr != corrValue_.end())
1228 cfac = itr->second;
1229 } else if (respCorrs_ != nullptr) {
1230 cfac = (respCorrs_->getValues(id))->getValue();
1231 }
1232 return cfac;
1233 }
1234
1235 double HcalHBHEMuonAnalyzer::gainFactor(const HcalDbService* conditions, const HcalDetId& id) {
1236 double gain(0.0);
1237 const HcalCalibrations& calibs = conditions->getHcalCalibrations(id);
1238 for (int capid = 0; capid < 4; ++capid)
1239 gain += (0.25 * calibs.respcorrgain(capid));
1240 return gain;
1241 }
1242
1243 int HcalHBHEMuonAnalyzer::depth16HE(int ieta, int iphi) {
1244
1245
1246
1247 int zside = (ieta > 0) ? 1 : -1;
1248 int depth = theHBHETopology_->dddConstants()->getMinDepth(1, 16, iphi, zside);
1249 if (isItPlan1_ && (!isItPreRecHit_))
1250 depth = 3;
1251 #ifdef EDM_ML_DEBUG
1252 edm::LogVerbatim("HBHEMuon") << "Plan1 " << isItPlan1_ << " PreRecHit " << isItPreRecHit_ << " phi " << iphi
1253 << " depth " << depth;
1254 #endif
1255 return depth;
1256 }
1257
1258 bool HcalHBHEMuonAnalyzer::goodCell(const HcalDetId& hcid,
1259 const reco::Track* pTrack,
1260 const CaloGeometry* geo,
1261 const MagneticField* bField) {
1262 std::pair<double, double> rz = hdc_->getRZ(hcid);
1263 bool typeRZ = (hcid.subdet() == HcalEndcap) ? false : true;
1264 bool match = spr::propagateHCAL(pTrack, geo, bField, typeRZ, rz, (((verbosity_ / 10000) % 10) > 0));
1265 return match;
1266 }
1267
1268
1269 #include "FWCore/Framework/interface/MakerMacros.h"
1270
1271 DEFINE_FWK_MODULE(HcalHBHEMuonAnalyzer);
