File indexing completed on 2022-10-14 01:43:21
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 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 #include "DataFormats/TrackReco/interface/HitPattern.h"
0029 #include "DataFormats/TrackReco/interface/TrackBase.h"
0030
0031
0032
0033 #include "DataFormats/Common/interface/TriggerResults.h"
0034 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0035 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0036 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0037 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0038 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0039 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0040 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0041 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0042
0043 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0044 #include "HLTrigger/HLTcore/interface/HLTConfigData.h"
0045
0046 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0047 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0048
0049 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0050 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0051 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0052
0053 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0054 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0055 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0056
0057 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0058 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0059 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0060 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0061 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0062 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0063 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0064 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0065 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0066 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0067 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0068 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0069 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0070 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0071 #include "MagneticField/Engine/interface/MagneticField.h"
0072 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0073
0074
0075
0076 class HcalHBHEMuonHighEtaAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0077 public:
0078 explicit HcalHBHEMuonHighEtaAnalyzer(const edm::ParameterSet&);
0079
0080 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0081
0082 private:
0083 void beginJob() override;
0084 void analyze(edm::Event const&, edm::EventSetup const&) override;
0085 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0086 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0087
0088 bool analyzeMuon(edm::Event const&, math::XYZPoint&);
0089 bool analyzeHadron(edm::Event const&, math::XYZPoint&);
0090 bool analyzeTracks(const reco::Track*, math::XYZPoint&, int, std::vector<spr::propagatedTrackID>&, bool);
0091 void clearVectors();
0092 int matchId(const HcalDetId&, const HcalDetId&);
0093 double activeLength(const DetId&);
0094 bool isGoodVertex(const reco::Vertex&);
0095 double respCorr(const DetId&);
0096 double gainFactor(const HcalDbService*, const HcalDetId&);
0097 int depth16HE(int, int);
0098 bool goodCell(const HcalDetId&, const reco::Track*, const CaloGeometry*, const MagneticField*);
0099 void fillTrackParameters(const reco::Track*, math::XYZPoint);
0100
0101
0102 const edm::InputTag labelEBRecHit_, labelEERecHit_, labelHBHERecHit_;
0103 const std::string labelVtx_, labelMuon_, labelGenTrack_;
0104 const double etaMin_, emaxNearPThr_;
0105 const bool analyzeMuon_, unCorrect_, collapseDepth_, isItPlan1_, getCharge_;
0106 const int useRaw_, verbosity_;
0107 const std::string theTrackQuality_, fileInCorr_;
0108 const bool ignoreHECorr_, isItPreRecHit_, writeRespCorr_;
0109 const int maxDepth_;
0110
0111 const edm::EDGetTokenT<reco::VertexCollection> tok_Vtx_;
0112 const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0113 const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0114 const edm::EDGetTokenT<HBHERecHitCollection> tok_HBHE_;
0115 const edm::EDGetTokenT<reco::MuonCollection> tok_Muon_;
0116 const edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0117
0118 const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec_;
0119 const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
0120 const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_respcorr_;
0121 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0122 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0123 const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_chan_;
0124 const edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0125 const edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_topo_;
0126 const edm::ESGetToken<HcalDbService, HcalDbRecord> tok_dbservice_;
0127
0128 bool mergedDepth_, useMyCorr_;
0129 int kount_;
0130 spr::trackSelectionParameters selectionParameter_;
0131
0132 const HcalDDDRecConstants* hdc_;
0133 const HcalTopology* theHBHETopology_;
0134 const CaloGeometry* geo_;
0135 HcalRespCorrs* respCorrs_;
0136 const MagneticField* bField_;
0137 const EcalChannelStatus* theEcalChStatus_;
0138 const EcalSeverityLevelAlgo* sevlv_;
0139 const CaloTopology* caloTopology_;
0140 const HcalDbService* conditions_;
0141
0142 edm::Handle<EcalRecHitCollection> barrelRecHitsHandle_;
0143 edm::Handle<EcalRecHitCollection> endcapRecHitsHandle_;
0144 edm::Handle<HBHERecHitCollection> hbhe_;
0145
0146
0147 static const int depthMax_ = 7;
0148 TTree* tree_;
0149 unsigned int runNumber_, eventNumber_, goodVertex_;
0150 std::vector<bool> mediumMuon_;
0151 std::vector<double> ptGlob_, etaGlob_, phiGlob_, energyMuon_, pMuon_;
0152 std::vector<double> isolationR04_, isolationR03_;
0153 std::vector<double> ecalEnergy_, hcalEnergy_, hoEnergy_;
0154 std::vector<bool> matchedId_, hcalHot_;
0155 std::vector<double> ecal3x3Energy_, hcal1x1Energy_;
0156 std::vector<unsigned int> ecalDetId_, hcalDetId_, ehcalDetId_;
0157 std::vector<int> hcal_ieta_, hcal_iphi_;
0158 std::vector<double> hcalDepthEnergy_[depthMax_];
0159 std::vector<double> hcalDepthActiveLength_[depthMax_];
0160 std::vector<double> hcalDepthEnergyHot_[depthMax_];
0161 std::vector<double> hcalDepthActiveLengthHot_[depthMax_];
0162 std::vector<double> hcalDepthChargeHot_[depthMax_];
0163 std::vector<double> hcalDepthChargeHotBG_[depthMax_];
0164 std::vector<double> hcalDepthEnergyCorr_[depthMax_];
0165 std::vector<double> hcalDepthEnergyHotCorr_[depthMax_];
0166 std::vector<bool> hcalDepthMatch_[depthMax_];
0167 std::vector<bool> hcalDepthMatchHot_[depthMax_];
0168 std::vector<double> hcalActiveLength_, hcalActiveLengthHot_;
0169 std::vector<double> emaxNearP_, trackDz_;
0170 std::vector<int> trackLayerCrossed_, trackOuterHit_;
0171 std::vector<int> trackMissedInnerHits_, trackMissedOuterHits_;
0172
0173 std::vector<HcalDDDRecConstants::HcalActiveLength> actHB, actHE;
0174 std::map<DetId, double> corrValue_;
0175
0176 };
0177
0178 HcalHBHEMuonHighEtaAnalyzer::HcalHBHEMuonHighEtaAnalyzer(const edm::ParameterSet& iConfig)
0179 : labelEBRecHit_(iConfig.getParameter<edm::InputTag>("labelEBRecHit")),
0180 labelEERecHit_(iConfig.getParameter<edm::InputTag>("labelEERecHit")),
0181 labelHBHERecHit_(iConfig.getParameter<edm::InputTag>("labelHBHERecHit")),
0182 labelVtx_(iConfig.getParameter<std::string>("labelVertex")),
0183 labelMuon_(iConfig.getParameter<std::string>("labelMuon")),
0184 labelGenTrack_(iConfig.getParameter<std::string>("labelTrack")),
0185 etaMin_(iConfig.getParameter<double>("etaMin")),
0186 emaxNearPThr_(iConfig.getParameter<double>("emaxNearPThreshold")),
0187 analyzeMuon_(iConfig.getParameter<bool>("analyzeMuon")),
0188 unCorrect_(iConfig.getParameter<bool>("unCorrect")),
0189 collapseDepth_(iConfig.getParameter<bool>("collapseDepth")),
0190 isItPlan1_(iConfig.getParameter<bool>("isItPlan1")),
0191 getCharge_(iConfig.getParameter<bool>("getCharge")),
0192 useRaw_(iConfig.getParameter<int>("useRaw")),
0193 verbosity_(iConfig.getParameter<int>("verbosity")),
0194 theTrackQuality_(iConfig.getUntrackedParameter<std::string>("trackQuality")),
0195 fileInCorr_(iConfig.getUntrackedParameter<std::string>("fileInCorr", "")),
0196 ignoreHECorr_(iConfig.getUntrackedParameter<bool>("ignoreHECorr", false)),
0197 isItPreRecHit_(iConfig.getUntrackedParameter<bool>("isItPreRecHit", false)),
0198 writeRespCorr_(iConfig.getUntrackedParameter<bool>("writeRespCorr", false)),
0199 maxDepth_(iConfig.getUntrackedParameter<int>("maxDepth", 7)),
0200 tok_Vtx_(consumes<reco::VertexCollection>(labelVtx_)),
0201 tok_EB_(consumes<EcalRecHitCollection>(labelEBRecHit_)),
0202 tok_EE_(consumes<EcalRecHitCollection>(labelEERecHit_)),
0203 tok_HBHE_(consumes<HBHERecHitCollection>(labelHBHERecHit_)),
0204 tok_Muon_(consumes<reco::MuonCollection>(labelMuon_)),
0205 tok_genTrack_(consumes<reco::TrackCollection>(labelGenTrack_)),
0206 tok_ddrec_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0207 tok_htopo_(esConsumes<HcalTopology, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0208 tok_respcorr_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd, edm::Transition::BeginRun>()),
0209 tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0210 tok_magField_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0211 tok_chan_(esConsumes<EcalChannelStatus, EcalChannelStatusRcd>()),
0212 tok_sevlv_(esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>()),
0213 tok_topo_(esConsumes<CaloTopology, CaloTopologyRecord>()),
0214 tok_dbservice_(esConsumes<HcalDbService, HcalDbRecord>()),
0215 hdc_(nullptr),
0216 theHBHETopology_(nullptr),
0217 respCorrs_(nullptr),
0218 tree_(nullptr) {
0219 usesResource(TFileService::kSharedResource);
0220
0221 kount_ = 0;
0222
0223 reco::TrackBase::TrackQuality trackQuality = reco::TrackBase::qualityByName(theTrackQuality_);
0224 selectionParameter_.minPt = iConfig.getUntrackedParameter<double>("minTrackPt");
0225 selectionParameter_.minQuality = trackQuality;
0226 selectionParameter_.maxDxyPV = iConfig.getUntrackedParameter<double>("maxDxyPV");
0227 selectionParameter_.maxDzPV = iConfig.getUntrackedParameter<double>("maxDzPV");
0228 selectionParameter_.maxChi2 = iConfig.getUntrackedParameter<double>("maxChi2");
0229 selectionParameter_.maxDpOverP = iConfig.getUntrackedParameter<double>("maxDpOverP");
0230 selectionParameter_.minOuterHit = selectionParameter_.minLayerCrossed = 0;
0231 selectionParameter_.maxInMiss = selectionParameter_.maxOutMiss = 2;
0232
0233 mergedDepth_ = (!isItPreRecHit_) || (collapseDepth_);
0234 edm::LogVerbatim("HBHEMuon") << "Labels used: Track " << labelGenTrack_ << " Vtx " << labelVtx_ << " EB "
0235 << labelEBRecHit_ << " EE " << labelEERecHit_ << " HBHE " << labelHBHERecHit_ << " MU "
0236 << labelMuon_;
0237
0238 if (!fileInCorr_.empty()) {
0239 std::ifstream infile(fileInCorr_.c_str());
0240 if (infile.is_open()) {
0241 while (true) {
0242 unsigned int id;
0243 double cfac;
0244 infile >> id >> cfac;
0245 if (!infile.good())
0246 break;
0247 corrValue_[DetId(id)] = cfac;
0248 }
0249 infile.close();
0250 }
0251 }
0252 useMyCorr_ = (!corrValue_.empty());
0253 edm::LogVerbatim("HBHEMuon") << "Flags used: UseRaw " << useRaw_ << " GetCharge " << getCharge_ << " UnCorrect "
0254 << unCorrect_ << " IgnoreHECorr " << ignoreHECorr_ << " CollapseDepth " << collapseDepth_
0255 << ":" << mergedDepth_ << " IsItPlan1 " << isItPlan1_ << " IsItPreRecHit "
0256 << isItPreRecHit_ << " UseMyCorr " << useMyCorr_;
0257 }
0258
0259
0260
0261
0262
0263
0264 void HcalHBHEMuonHighEtaAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0265 edm::ParameterSetDescription desc;
0266 desc.add<edm::InputTag>("labelEBRecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0267 desc.add<edm::InputTag>("labelEERecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0268 desc.add<edm::InputTag>("labelHBHERecHit", edm::InputTag("hbhereco"));
0269 desc.add<std::string>("labelVertex", "offlinePrimaryVertices");
0270 desc.add<std::string>("labelMuon", "muons");
0271 desc.add<std::string>("labelTrack", "generalTracks");
0272 desc.add<double>("etaMin", 2.0);
0273 desc.add<double>("emaxNearPThreshold", 10.0);
0274 desc.add<bool>("analyzeMuon", true);
0275 desc.add<bool>("unCorrect", true);
0276 desc.add<bool>("collapseDepth", false);
0277 desc.add<bool>("isItPlan1", false);
0278 desc.add<bool>("getCharge", true);
0279 desc.add<int>("useRaw", 0);
0280 desc.add<int>("verbosity", 0);
0281 desc.addUntracked<std::string>("fileInCorr", "");
0282 desc.addUntracked<std::string>("trackQuality", "highPurity");
0283 desc.addUntracked<double>("minTrackPt", 1.0);
0284 desc.addUntracked<double>("maxDxyPV", 0.02);
0285 desc.addUntracked<double>("maxDzPV", 100.0);
0286 desc.addUntracked<double>("maxChi2", 5.0);
0287 desc.addUntracked<double>("maxDpOverP", 0.1);
0288 desc.addUntracked<bool>("ignoreHECorr", false);
0289 desc.addUntracked<bool>("isItPreRecHit", false);
0290 desc.addUntracked<bool>("writeRespCorr", false);
0291 desc.addUntracked<int>("maxDepth", 7);
0292 descriptions.add("hcalHBHEMuonHighEta", desc);
0293 }
0294
0295
0296 void HcalHBHEMuonHighEtaAnalyzer::beginJob() {
0297 edm::Service<TFileService> fs;
0298 tree_ = fs->make<TTree>("HBHEMuonHighEta", "HBHEMuonHighEta");
0299 tree_->Branch("pt_of_muon", &ptGlob_);
0300 tree_->Branch("eta_of_muon", &etaGlob_);
0301 tree_->Branch("phi_of_muon", &phiGlob_);
0302 tree_->Branch("energy_of_muon", &energyMuon_);
0303 tree_->Branch("p_of_muon", &pMuon_);
0304 tree_->Branch("MediumMuon", &mediumMuon_);
0305 tree_->Branch("IsolationR04", &isolationR04_);
0306 tree_->Branch("IsolationR03", &isolationR03_);
0307 tree_->Branch("ecal_3into3", &ecalEnergy_);
0308 tree_->Branch("hcal_3into3", &hcalEnergy_);
0309 tree_->Branch("ho_3into3", &hoEnergy_);
0310 tree_->Branch("emaxNearP", &emaxNearP_);
0311
0312 tree_->Branch("Run_No", &runNumber_);
0313 tree_->Branch("Event_No", &eventNumber_);
0314 tree_->Branch("GoodVertex", &goodVertex_);
0315 tree_->Branch("matchedId", &matchedId_);
0316 tree_->Branch("hcal_cellHot", &hcalHot_);
0317 tree_->Branch("ecal_3x3", &ecal3x3Energy_);
0318 tree_->Branch("hcal_1x1", &hcal1x1Energy_);
0319 tree_->Branch("ecal_detID", &ecalDetId_);
0320 tree_->Branch("hcal_detID", &hcalDetId_);
0321 tree_->Branch("ehcal_detID", &ehcalDetId_);
0322 tree_->Branch("hcal_ieta", &hcal_ieta_);
0323 tree_->Branch("hcal_iphi", &hcal_iphi_);
0324
0325 char name[100];
0326 for (int k = 0; k < maxDepth_; ++k) {
0327 sprintf(name, "hcal_edepth%d", (k + 1));
0328 tree_->Branch(name, &hcalDepthEnergy_[k]);
0329 sprintf(name, "hcal_activeL%d", (k + 1));
0330 tree_->Branch(name, &hcalDepthActiveLength_[k]);
0331 sprintf(name, "hcal_edepthHot%d", (k + 1));
0332 tree_->Branch(name, &hcalDepthEnergyHot_[k]);
0333 sprintf(name, "hcal_activeHotL%d", (k + 1));
0334 tree_->Branch(name, &hcalDepthActiveLengthHot_[k]);
0335 sprintf(name, "hcal_cdepthHot%d", (k + 1));
0336 tree_->Branch(name, &hcalDepthChargeHot_[k]);
0337 sprintf(name, "hcal_cdepthHotBG%d", (k + 1));
0338 tree_->Branch(name, &hcalDepthChargeHotBG_[k]);
0339 sprintf(name, "hcal_edepthCorrect%d", (k + 1));
0340 tree_->Branch(name, &hcalDepthEnergyCorr_[k]);
0341 sprintf(name, "hcal_edepthHotCorrect%d", (k + 1));
0342 tree_->Branch(name, &hcalDepthEnergyHotCorr_[k]);
0343 sprintf(name, "hcal_depthMatch%d", (k + 1));
0344 tree_->Branch(name, &hcalDepthMatch_[k]);
0345 sprintf(name, "hcal_depthMatchHot%d", (k + 1));
0346 tree_->Branch(name, &hcalDepthMatchHot_[k]);
0347 }
0348 tree_->Branch("activeLength", &hcalActiveLength_);
0349 tree_->Branch("activeLengthHot", &hcalActiveLengthHot_);
0350 tree_->Branch("trackDz", &trackDz_);
0351 tree_->Branch("trackLayerCrossed", &trackLayerCrossed_);
0352 tree_->Branch("trackOuterHit", &trackOuterHit_);
0353 tree_->Branch("trackMissedInnerHits", &trackMissedInnerHits_);
0354 tree_->Branch("trackMissedOuterHits", &trackMissedOuterHits_);
0355 }
0356
0357
0358 void HcalHBHEMuonHighEtaAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0359 ++kount_;
0360 clearVectors();
0361 runNumber_ = iEvent.id().run();
0362 eventNumber_ = iEvent.id().event();
0363 #ifdef EDM_ML_DEBUG
0364 edm::LogVerbatim("HBHEMuon") << "Run " << runNumber_ << " Event " << eventNumber_;
0365 #endif
0366
0367
0368 bField_ = &iSetup.getData(tok_magField_);
0369 theEcalChStatus_ = &iSetup.getData(tok_chan_);
0370 sevlv_ = &iSetup.getData(tok_sevlv_);
0371 caloTopology_ = &iSetup.getData(tok_topo_);
0372 conditions_ = &iSetup.getData(tok_dbservice_);
0373
0374
0375 const edm::Handle<reco::VertexCollection>& vtx = iEvent.getHandle(tok_Vtx_);
0376
0377 iEvent.getByToken(tok_EB_, barrelRecHitsHandle_);
0378 iEvent.getByToken(tok_EE_, endcapRecHitsHandle_);
0379 iEvent.getByToken(tok_HBHE_, hbhe_);
0380
0381
0382 math::XYZPoint pvx;
0383 goodVertex_ = 0;
0384 if (!vtx.isValid()) {
0385 #ifdef EDM_ML_DEBUG
0386 edm::LogVerbatim("HBHEMuon") << "No Good Vertex found == Reject\n";
0387 #endif
0388 return;
0389 }
0390
0391 reco::VertexCollection::const_iterator firstGoodVertex = vtx->end();
0392 for (reco::VertexCollection::const_iterator it = vtx->begin(); it != vtx->end(); it++) {
0393 if (isGoodVertex(*it)) {
0394 if (firstGoodVertex == vtx->end())
0395 firstGoodVertex = it;
0396 ++goodVertex_;
0397 }
0398 }
0399 if (firstGoodVertex != vtx->end())
0400 pvx = firstGoodVertex->position();
0401
0402 bool accept(false);
0403 if (barrelRecHitsHandle_.isValid() && endcapRecHitsHandle_.isValid() && hbhe_.isValid()) {
0404 accept = analyzeMuon_ ? analyzeMuon(iEvent, pvx) : analyzeHadron(iEvent, pvx);
0405 }
0406 if (accept) {
0407 #ifdef EDM_ML_DEBUG
0408 edm::LogVerbatim("HBHEMuon") << "Total of " << hcal_ieta_.size() << " propagated points";
0409 for (unsigned int i = 0; i < hcal_ieta_.size(); ++i)
0410 edm::LogVerbatim("HBHEMuon") << "[" << i << "] ieta/iphi for entry to "
0411 << "HCAL has value of " << hcal_ieta_[i] << ":" << hcal_iphi_[i];
0412 if ((verbosity_ / 100) % 10 > 0) {
0413 edm::LogVerbatim("HBHEMuon") << "Sizes:: ptGlob:" << ptGlob_.size() << " etaGlob:" << etaGlob_.size()
0414 << " phiGlob:" << phiGlob_.size() << " energyMuon:" << energyMuon_.size()
0415 << " pMuon:" << pMuon_.size() << " mediumMuon: " << mediumMuon_.size()
0416 << " isolation:" << isolationR04_.size() << ":" << isolationR03_.size()
0417 << " e|h|ho energy: " << ecalEnergy_.size() << ":" << hcalEnergy_.size() << ":"
0418 << hoEnergy_.size();
0419 edm::LogVerbatim("HBHEMuon") << " matchedId:" << matchedId_.size() << " hcalHot:" << hcalHot_.size()
0420 << " 3x3|1x1 energy:" << ecal3x3Energy_.size() << ":" << hcal1x1Energy_.size()
0421 << " detId:" << ecalDetId_.size() << ":" << hcalDetId_.size() << ":"
0422 << ehcalDetId_.size() << " eta|phi:" << hcal_ieta_.size() << ":"
0423 << hcal_iphi_.size();
0424 edm::LogVerbatim("HBHEMuon") << " activeLength:" << hcalActiveLength_.size() << ":"
0425 << hcalActiveLengthHot_.size() << " emaxNearP:" << emaxNearP_.size()
0426 << " trackDz: " << trackDz_.size() << " tracks:" << trackLayerCrossed_.size() << ":"
0427 << trackOuterHit_.size() << ":" << trackMissedInnerHits_.size() << ":"
0428 << trackMissedOuterHits_.size();
0429 for (unsigned int i = 0; i < depthMax_; ++i)
0430 edm::LogVerbatim("HBHEMuon")
0431 << "Depth " << i
0432 << " Energy|Length|EnergyHot|LengthHot|Charge|ChargeBG|EnergyCorr|EnergyHotCorr|Match|MatchHot:"
0433 << hcalDepthEnergy_[i].size() << ":" << hcalDepthActiveLength_[i].size() << ":"
0434 << hcalDepthEnergyHot_[i].size() << ":" << hcalDepthActiveLengthHot_[i].size() << ":"
0435 << hcalDepthChargeHot_[i].size() << ":" << hcalDepthChargeHotBG_[i].size() << ":"
0436 << hcalDepthEnergyCorr_[i].size() << ":" << hcalDepthEnergyHotCorr_[i].size() << ":"
0437 << hcalDepthMatch_[i].size() << ":" << hcalDepthMatchHot_[i].size();
0438 }
0439 #endif
0440 tree_->Fill();
0441 }
0442 }
0443
0444
0445 void HcalHBHEMuonHighEtaAnalyzer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0446 hdc_ = &iSetup.getData(tok_ddrec_);
0447 actHB.clear();
0448 actHE.clear();
0449 actHB = hdc_->getThickActive(0);
0450 actHE = hdc_->getThickActive(1);
0451 #ifdef EDM_ML_DEBUG
0452 if (verbosity_ % 10 > 0) {
0453 unsigned int k1(0), k2(0);
0454 edm::LogVerbatim("HBHEMuon") << actHB.size() << " Active Length for HB";
0455 for (const auto& act : actHB) {
0456 edm::LogVerbatim("HBHEMuon") << "[" << k1 << "] ieta " << act.ieta << " depth " << act.depth << " zside "
0457 << act.zside << " type " << act.stype << " phi " << act.iphis.size() << ":"
0458 << act.iphis[0] << " L " << act.thick;
0459 HcalDetId hcid1(HcalBarrel, (act.ieta) * (act.zside), act.iphis[0], act.depth);
0460 HcalDetId hcid2 = mergedDepth_ ? hdc_->mergedDepthDetId(hcid1) : hcid1;
0461 edm::LogVerbatim("HBHEMuon") << hcid1 << " | " << hcid2 << " L " << activeLength(DetId(hcid2));
0462 ++k1;
0463 }
0464 edm::LogVerbatim("HBHEMuon") << actHE.size() << " Active Length for HE";
0465 for (const auto& act : actHE) {
0466 edm::LogVerbatim("HBHEMuon") << "[" << k2 << "] ieta " << act.ieta << " depth " << act.depth << " zside "
0467 << act.zside << " type " << act.stype << " phi " << act.iphis.size() << ":"
0468 << act.iphis[0] << " L " << act.thick;
0469 HcalDetId hcid1(HcalEndcap, (act.ieta) * (act.zside), act.iphis[0], act.depth);
0470 HcalDetId hcid2 = mergedDepth_ ? hdc_->mergedDepthDetId(hcid1) : hcid1;
0471 edm::LogVerbatim("HBHEMuon") << hcid1 << " | " << hcid2 << " L " << activeLength(DetId(hcid2));
0472 ++k2;
0473 }
0474 }
0475 #endif
0476
0477 theHBHETopology_ = &iSetup.getData(tok_htopo_);
0478 const HcalRespCorrs* resp = &iSetup.getData(tok_respcorr_);
0479 respCorrs_ = new HcalRespCorrs(*resp);
0480 respCorrs_->setTopo(theHBHETopology_);
0481 geo_ = &iSetup.getData(tok_geom_);
0482
0483
0484 if (writeRespCorr_) {
0485 const HcalGeometry* gHcal = (const HcalGeometry*)(geo_->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0486 const std::vector<DetId>& ids = gHcal->getValidDetIds(DetId::Hcal, 0);
0487 edm::LogVerbatim("HBHEMuon") << "\nTable of Correction Factors for Run " << iRun.run() << "\n";
0488 for (auto const& id : ids) {
0489 if ((id.det() == DetId::Hcal) && ((id.subdetId() == HcalBarrel) || (id.subdetId() == HcalEndcap))) {
0490 edm::LogVerbatim("HBHEMuon") << HcalDetId(id) << " " << id.rawId() << " "
0491 << (respCorrs_->getValues(id))->getValue();
0492 }
0493 }
0494 }
0495 }
0496
0497 bool HcalHBHEMuonHighEtaAnalyzer::analyzeMuon(const edm::Event& iEvent, math::XYZPoint& leadPV) {
0498 const edm::Handle<reco::MuonCollection>& _Muon = iEvent.getHandle(tok_Muon_);
0499 bool accept = false;
0500
0501 if (_Muon.isValid()) {
0502 int nTrack(0);
0503 std::vector<spr::propagatedTrackID> trkCaloDets;
0504 for (const auto& RecMuon : (*(_Muon.product()))) {
0505 if (RecMuon.innerTrack().isNonnull()) {
0506 const reco::Track* pTrack = (RecMuon.innerTrack()).get();
0507 if (std::abs(pTrack->eta()) > etaMin_) {
0508 if (analyzeTracks(pTrack, leadPV, nTrack, trkCaloDets, false)) {
0509 accept = true;
0510 ptGlob_.emplace_back(RecMuon.pt());
0511 etaGlob_.emplace_back(RecMuon.eta());
0512 phiGlob_.emplace_back(RecMuon.phi());
0513 energyMuon_.push_back(RecMuon.energy());
0514 pMuon_.emplace_back(RecMuon.p());
0515 bool mediumMuon = (((RecMuon.isPFMuon()) && (RecMuon.isGlobalMuon() || RecMuon.isTrackerMuon())) &&
0516 (RecMuon.innerTrack()->validFraction() > 0.49));
0517 if (mediumMuon) {
0518 double chiGlobal = ((RecMuon.globalTrack().isNonnull()) ? RecMuon.globalTrack()->normalizedChi2() : 999);
0519 bool goodGlob =
0520 (RecMuon.isGlobalMuon() && chiGlobal < 3 && RecMuon.combinedQuality().chi2LocalPosition < 12 &&
0521 RecMuon.combinedQuality().trkKink < 20);
0522 mediumMuon = muon::segmentCompatibility(RecMuon) > (goodGlob ? 0.303 : 0.451);
0523 }
0524 mediumMuon_.emplace_back(mediumMuon);
0525 bool isoR03 =
0526 ((RecMuon.pfIsolationR03().sumChargedHadronPt +
0527 std::max(0.,
0528 RecMuon.pfIsolationR03().sumNeutralHadronEt + RecMuon.pfIsolationR03().sumPhotonEt -
0529 (0.5 * RecMuon.pfIsolationR03().sumPUPt))) /
0530 RecMuon.pt());
0531 bool isoR04 =
0532 ((RecMuon.pfIsolationR04().sumChargedHadronPt +
0533 std::max(0.,
0534 RecMuon.pfIsolationR04().sumNeutralHadronEt + RecMuon.pfIsolationR04().sumPhotonEt -
0535 (0.5 * RecMuon.pfIsolationR04().sumPUPt))) /
0536 RecMuon.pt());
0537 isolationR03_.emplace_back(isoR03);
0538 isolationR04_.emplace_back(isoR04);
0539
0540 ecalEnergy_.emplace_back(RecMuon.calEnergy().emS9);
0541 hcalEnergy_.emplace_back(RecMuon.calEnergy().hadS9);
0542 hoEnergy_.emplace_back(RecMuon.calEnergy().hoS9);
0543 #ifdef EDM_ML_DEBUG
0544 if ((verbosity_ / 100) % 10 > 0)
0545 edm::LogVerbatim("HBHEMuon")
0546 << "Muon[" << ptGlob_.size() << "] pt:eta:phi:p " << ptGlob_.back() << ":" << etaGlob_.back() << ":"
0547 << phiGlob_.back() << ":" << energyMuon_.back() << ":" << pMuon_.back() << ":"
0548 << " Medium:i3:i4 " << mediumMuon_.back() << ":" << isolationR03_.back() << ":"
0549 << isolationR04_.back() << ":"
0550 << " Energy EC:HC:HO " << ecalEnergy_.back() << ":" << hcalEnergy_.back() << ":" << hoEnergy_.back();
0551 #endif
0552 }
0553 }
0554 }
0555 }
0556 }
0557 return accept;
0558 }
0559
0560 bool HcalHBHEMuonHighEtaAnalyzer::analyzeHadron(const edm::Event& iEvent, math::XYZPoint& leadPV) {
0561
0562 edm::Handle<reco::TrackCollection> trkCollection = iEvent.getHandle(tok_genTrack_);
0563 bool accept = false;
0564
0565 if (!trkCollection.isValid()) {
0566 std::vector<spr::propagatedTrackID> trkCaloDets;
0567 spr::propagateCALO(trkCollection, geo_, bField_, theTrackQuality_, trkCaloDets, false);
0568 int nTrack(0);
0569 std::vector<spr::propagatedTrackID>::const_iterator trkDetItr;
0570 for (trkDetItr = trkCaloDets.begin(), nTrack = 0; trkDetItr != trkCaloDets.end(); trkDetItr++, nTrack++) {
0571 const reco::Track* pTrack = &(*(trkDetItr->trkItr));
0572 if (std::abs(pTrack->eta()) > etaMin_) {
0573 accept = analyzeTracks(pTrack, leadPV, nTrack, trkCaloDets, true);
0574 }
0575 }
0576 }
0577 return accept;
0578 }
0579
0580 bool HcalHBHEMuonHighEtaAnalyzer::analyzeTracks(const reco::Track* pTrack,
0581 math::XYZPoint& leadPV,
0582 int nTrack,
0583 std::vector<spr::propagatedTrackID>& trkCaloDets,
0584 bool ifHadron) {
0585 bool accept(false);
0586
0587 if (spr::goodTrack(pTrack, leadPV, selectionParameter_, false)) {
0588 spr::propagatedTrackID trackID = spr::propagateCALO(pTrack, geo_, bField_, false);
0589
0590 if (trackID.okECAL && trackID.okHCAL) {
0591 double emaxNearP = (ifHadron) ? spr::chargeIsolationEcal(nTrack, trkCaloDets, geo_, caloTopology_, 15, 15) : 0;
0592 if (emaxNearP < emaxNearPThr_) {
0593 double eEcal(0), eHcal(0), activeLengthTot(0), activeLengthHotTot(0);
0594 double eHcalDepth[depthMax_], eHcalDepthHot[depthMax_];
0595 double eHcalDepthC[depthMax_], eHcalDepthHotC[depthMax_];
0596 double cHcalDepthHot[depthMax_], cHcalDepthHotBG[depthMax_];
0597 double activeL[depthMax_], activeHotL[depthMax_];
0598 bool matchDepth[depthMax_], matchDepthHot[depthMax_];
0599 HcalDetId eHcalDetId[depthMax_];
0600 unsigned int isHot(0);
0601 bool tmpmatch(false);
0602 int ieta(-1000), iphi(-1000);
0603 for (int i = 0; i < depthMax_; ++i) {
0604 eHcalDepth[i] = eHcalDepthHot[i] = 0;
0605 eHcalDepthC[i] = eHcalDepthHotC[i] = 0;
0606 cHcalDepthHot[i] = cHcalDepthHotBG[i] = 0;
0607 activeL[i] = activeHotL[i] = 0;
0608 matchDepth[i] = matchDepthHot[i] = true;
0609 }
0610
0611 HcalDetId check;
0612 std::pair<bool, HcalDetId> info = spr::propagateHCALBack(pTrack, geo_, bField_, false);
0613 if (info.first)
0614 check = info.second;
0615
0616 const DetId isoCell(trackID.detIdECAL);
0617 std::pair<double, bool> e3x3 = spr::eECALmatrix(isoCell,
0618 barrelRecHitsHandle_,
0619 endcapRecHitsHandle_,
0620 *theEcalChStatus_,
0621 geo_,
0622 caloTopology_,
0623 sevlv_,
0624 1,
0625 1,
0626 -100.0,
0627 -100.0,
0628 -500.0,
0629 500.0,
0630 false);
0631 eEcal = e3x3.first;
0632 #ifdef EDM_ML_DEBUG
0633 if (verbosity_ % 10 > 0)
0634 edm::LogVerbatim("HBHEMuon") << "Propagate Track to ECAL: " << e3x3.second << ":" << trackID.okECAL << " E "
0635 << eEcal;
0636 #endif
0637
0638 DetId closestCell(trackID.detIdHCAL);
0639 HcalDetId hcidt(closestCell.rawId());
0640 if ((hcidt.ieta() == check.ieta()) && (hcidt.iphi() == check.iphi()))
0641 tmpmatch = true;
0642 #ifdef EDM_ML_DEBUG
0643 if (verbosity_ % 10 > 0)
0644 edm::LogVerbatim("HBHEMuon") << "Front " << hcidt << " Back " << info.first << ":" << check << " Match "
0645 << tmpmatch;
0646 #endif
0647
0648 HcalSubdetector subdet = hcidt.subdet();
0649 ieta = hcidt.ieta();
0650 iphi = hcidt.iphi();
0651 bool hborhe = (std::abs(ieta) == 16);
0652
0653 eHcal = spr::eHCALmatrix(theHBHETopology_,
0654 closestCell,
0655 hbhe_,
0656 0,
0657 0,
0658 false,
0659 true,
0660 -100.0,
0661 -100.0,
0662 -100.0,
0663 -100.0,
0664 -500.,
0665 500.,
0666 useRaw_);
0667 std::vector<std::pair<double, int>> ehdepth;
0668 spr::energyHCALCell((HcalDetId)closestCell,
0669 hbhe_,
0670 ehdepth,
0671 depthMax_,
0672 -100.0,
0673 -100.0,
0674 -100.0,
0675 -100.0,
0676 -500.0,
0677 500.0,
0678 useRaw_,
0679 depth16HE(ieta, iphi),
0680 false);
0681 for (int i = 0; i < depthMax_; ++i)
0682 eHcalDetId[i] = HcalDetId();
0683 for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0684 HcalSubdetector subdet0 =
0685 (hborhe) ? ((ehdepth[i].second >= depth16HE(ieta, iphi)) ? HcalEndcap : HcalBarrel) : subdet;
0686 HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
0687 double actL = activeLength(DetId(hcid0));
0688 double ene = ehdepth[i].first;
0689 bool tmpC(false);
0690 if (ene > 0.0) {
0691 if (!(theHBHETopology_->validHcal(hcid0))) {
0692 edm::LogWarning("HBHEMuon") << "(1) Invalid ID " << hcid0 << " with E = " << ene;
0693 edm::LogWarning("HBHEMuon") << HcalDetId(closestCell) << " with " << ehdepth.size() << " depths:";
0694 for (const auto& ehd : ehdepth)
0695 edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0696 } else {
0697 tmpC = goodCell(hcid0, pTrack, geo_, bField_);
0698 double enec(ene);
0699 if (unCorrect_) {
0700 double corr = (ignoreHECorr_ && (subdet0 == HcalEndcap)) ? 1.0 : respCorr(DetId(hcid0));
0701 if (corr != 0)
0702 ene /= corr;
0703 #ifdef EDM_ML_DEBUG
0704 if (verbosity_ % 10 > 0) {
0705 HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc_->mergedDepthDetId(hcid0) : hcid0;
0706 edm::LogVerbatim("HBHEMuon") << hcid0 << ":" << id << " Corr " << corr;
0707 }
0708 #endif
0709 }
0710 int depth = ehdepth[i].second - 1;
0711 if (collapseDepth_) {
0712 HcalDetId id = hdc_->mergedDepthDetId(hcid0);
0713 depth = id.depth() - 1;
0714 }
0715 eHcalDepth[depth] += ene;
0716 eHcalDepthC[depth] += enec;
0717 activeL[depth] += actL;
0718 activeLengthTot += actL;
0719 matchDepth[depth] = (matchDepth[depth] && tmpC);
0720 #ifdef EDM_ML_DEBUG
0721 if ((verbosity_ / 10) % 10 > 0)
0722 edm::LogVerbatim("HBHEMuon")
0723 << hcid0 << " E " << ene << ":" << enec << " L " << actL << " Match " << tmpC;
0724 #endif
0725 }
0726 }
0727 }
0728 #ifdef EDM_ML_DEBUG
0729 if ((verbosity_ / 10) % 10 > 0) {
0730 edm::LogVerbatim("HBHEMuon") << hcidt << " Match " << tmpmatch << " Depths " << ehdepth.size();
0731 for (unsigned int k = 0; k < ehdepth.size(); ++k)
0732 edm::LogVerbatim("HBHEMuon") << " [" << k << ":" << ehdepth[k].second << "] " << matchDepth[k];
0733 }
0734 #endif
0735 HcalDetId hotCell;
0736 spr::eHCALmatrix(geo_, theHBHETopology_, closestCell, hbhe_, 1, 1, hotCell, false, useRaw_, false);
0737 isHot = matchId(closestCell, hotCell);
0738 if (hotCell != HcalDetId()) {
0739 subdet = HcalDetId(hotCell).subdet();
0740 ieta = HcalDetId(hotCell).ieta();
0741 iphi = HcalDetId(hotCell).iphi();
0742 hborhe = (std::abs(ieta) == 16);
0743 std::vector<std::pair<double, int>> ehdepth;
0744 spr::energyHCALCell(hotCell,
0745 hbhe_,
0746 ehdepth,
0747 depthMax_,
0748 -100.0,
0749 -100.0,
0750 -100.0,
0751 -100.0,
0752 -500.0,
0753 500.0,
0754 useRaw_,
0755 depth16HE(ieta, iphi),
0756 false);
0757 for (int i = 0; i < depthMax_; ++i)
0758 eHcalDetId[i] = HcalDetId();
0759 for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0760 HcalSubdetector subdet0 =
0761 (hborhe) ? ((ehdepth[i].second >= depth16HE(ieta, iphi)) ? HcalEndcap : HcalBarrel) : subdet;
0762 HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
0763 double actL = activeLength(DetId(hcid0));
0764 double ene = ehdepth[i].first;
0765 bool tmpC(false);
0766 if (ene > 0.0) {
0767 if (!(theHBHETopology_->validHcal(hcid0))) {
0768 edm::LogWarning("HBHEMuon") << "(2) Invalid ID " << hcid0 << " with E = " << ene;
0769 edm::LogWarning("HBHEMuon") << HcalDetId(hotCell) << " with " << ehdepth.size() << " depths:";
0770 for (const auto& ehd : ehdepth)
0771 edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0772 } else {
0773 tmpC = goodCell(hcid0, pTrack, geo_, bField_);
0774 double chg(ene), enec(ene);
0775 if (unCorrect_) {
0776 double corr = (ignoreHECorr_ && (subdet0 == HcalEndcap)) ? 1.0 : respCorr(DetId(hcid0));
0777 if (corr != 0)
0778 ene /= corr;
0779 #ifdef EDM_ML_DEBUG
0780 if (verbosity_ % 10 > 0) {
0781 HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc_->mergedDepthDetId(hcid0) : hcid0;
0782 edm::LogVerbatim("HBHEMuon")
0783 << hcid0 << ":" << id << " Corr " << corr << " E " << ene << ":" << enec;
0784 }
0785 #endif
0786 }
0787 if (getCharge_) {
0788 double gain = gainFactor(conditions_, hcid0);
0789 if (gain != 0)
0790 chg /= gain;
0791 #ifdef EDM_ML_DEBUG
0792 if (verbosity_ % 10 > 0)
0793 edm::LogVerbatim("HBHEMuon") << hcid0 << " Gain " << gain << " C " << chg;
0794 #endif
0795 }
0796 int depth = ehdepth[i].second - 1;
0797 if (collapseDepth_) {
0798 HcalDetId id = hdc_->mergedDepthDetId(hcid0);
0799 depth = id.depth() - 1;
0800 }
0801 eHcalDepthHot[depth] += ene;
0802 eHcalDepthHotC[depth] += enec;
0803 cHcalDepthHot[depth] += chg;
0804 activeHotL[depth] += actL;
0805 activeLengthHotTot += actL;
0806 matchDepthHot[depth] = (matchDepthHot[depth] && tmpC);
0807 #ifdef EDM_ML_DEBUG
0808 if ((verbosity_ / 10) % 10 > 0)
0809 edm::LogVerbatim("HBHEMuon") << hcid0 << " depth " << depth << " E " << ene << ":" << enec << " C "
0810 << chg << " L " << actL << " Match " << tmpC;
0811 #endif
0812 }
0813 }
0814 }
0815 }
0816 #ifdef EDM_ML_DEBUG
0817 edm::LogVerbatim("HBHEMuon") << "Propagate Track to HCAL: " << trackID.okHCAL << " Match " << tmpmatch
0818 << " Hot " << isHot << " Energy " << eHcal;
0819 #endif
0820
0821 accept = true;
0822 ecalDetId_.emplace_back((trackID.detIdECAL)());
0823 hcalDetId_.emplace_back((trackID.detIdHCAL)());
0824 ehcalDetId_.emplace_back((trackID.detIdEHCAL)());
0825 emaxNearP_.emplace_back(emaxNearP);
0826 matchedId_.emplace_back(tmpmatch);
0827 ecal3x3Energy_.emplace_back(eEcal);
0828 hcal1x1Energy_.emplace_back(eHcal);
0829 hcal_ieta_.emplace_back(ieta);
0830 hcal_iphi_.emplace_back(iphi);
0831 for (int i = 0; i < maxDepth_; ++i) {
0832 hcalDepthEnergy_[i].emplace_back(eHcalDepth[i]);
0833 hcalDepthActiveLength_[i].emplace_back(activeL[i]);
0834 hcalDepthEnergyHot_[i].emplace_back(eHcalDepthHot[i]);
0835 hcalDepthActiveLengthHot_[i].emplace_back(activeHotL[i]);
0836 hcalDepthEnergyCorr_[i].emplace_back(eHcalDepthC[i]);
0837 hcalDepthEnergyHotCorr_[i].emplace_back(eHcalDepthHotC[i]);
0838 hcalDepthChargeHot_[i].emplace_back(cHcalDepthHot[i]);
0839 hcalDepthChargeHotBG_[i].emplace_back(cHcalDepthHotBG[i]);
0840 hcalDepthMatch_[i].emplace_back(matchDepth[i]);
0841 hcalDepthMatchHot_[i].emplace_back(matchDepthHot[i]);
0842 }
0843 hcalActiveLength_.emplace_back(activeLengthTot);
0844 hcalHot_.emplace_back(isHot);
0845 hcalActiveLengthHot_.emplace_back(activeLengthHotTot);
0846 #ifdef EDM_ML_DEBUG
0847 if ((verbosity_ / 100) % 10 > 0) {
0848 edm::LogVerbatim("HBHEMuon") << "Track " << std::hex << ecalDetId_.back() << ":" << hcalDetId_.back() << ":"
0849 << ehcalDetId_.back() << std::dec << ":" << emaxNearP_.back() << ":"
0850 << matchedId_.back() << ":" << ecal3x3Energy_.back() << ":"
0851 << hcal1x1Energy_.back() << ":" << hcal_ieta_.back() << ":" << hcal_iphi_.back()
0852 << ":" << hcalActiveLength_.back() << ":" << hcalHot_.back() << ":"
0853 << hcalActiveLengthHot_.back();
0854 for (int i = 0; i < maxDepth_; ++i) {
0855 edm::LogVerbatim("HBHEMuon") << "Depth[" << i << "] " << hcalDepthEnergy_[i].back() << ":"
0856 << hcalDepthActiveLength_[i].back() << ":" << hcalDepthEnergyHot_[i].back()
0857 << ":" << hcalDepthActiveLengthHot_[i].back() << ":"
0858 << hcalDepthEnergyCorr_[i].back() << ":" << hcalDepthEnergyHotCorr_[i].back()
0859 << ":" << hcalDepthChargeHot_[i].back() << ":"
0860 << hcalDepthChargeHotBG_[i].back() << ":" << hcalDepthMatch_[i].back() << ":"
0861 << hcalDepthMatchHot_[i].back();
0862 }
0863 }
0864 #endif
0865 fillTrackParameters(pTrack, leadPV);
0866 }
0867 }
0868 }
0869 return accept;
0870 }
0871
0872 void HcalHBHEMuonHighEtaAnalyzer::clearVectors() {
0873
0874 eventNumber_ = -99999;
0875 runNumber_ = -99999;
0876 goodVertex_ = -99999;
0877
0878 mediumMuon_.clear();
0879 ptGlob_.clear();
0880 etaGlob_.clear();
0881 phiGlob_.clear();
0882 energyMuon_.clear();
0883 pMuon_.clear();
0884 isolationR04_.clear();
0885 isolationR03_.clear();
0886 ecalEnergy_.clear();
0887 hcalEnergy_.clear();
0888 hoEnergy_.clear();
0889
0890 matchedId_.clear();
0891 hcalHot_.clear();
0892 ecal3x3Energy_.clear();
0893 hcal1x1Energy_.clear();
0894 ecalDetId_.clear();
0895 hcalDetId_.clear();
0896 ehcalDetId_.clear();
0897 hcal_ieta_.clear();
0898 hcal_iphi_.clear();
0899 for (int i = 0; i < depthMax_; ++i) {
0900 hcalDepthEnergy_[i].clear();
0901 hcalDepthActiveLength_[i].clear();
0902 hcalDepthEnergyHot_[i].clear();
0903 hcalDepthActiveLengthHot_[i].clear();
0904 hcalDepthChargeHot_[i].clear();
0905 hcalDepthChargeHotBG_[i].clear();
0906 hcalDepthEnergyCorr_[i].clear();
0907 hcalDepthEnergyHotCorr_[i].clear();
0908 hcalDepthMatch_[i].clear();
0909 hcalDepthMatchHot_[i].clear();
0910 }
0911 hcalActiveLength_.clear();
0912 hcalActiveLengthHot_.clear();
0913
0914 emaxNearP_.clear();
0915 trackDz_.clear();
0916 trackLayerCrossed_.clear();
0917 trackOuterHit_.clear();
0918 trackMissedInnerHits_.clear();
0919 trackMissedOuterHits_.clear();
0920 }
0921
0922 int HcalHBHEMuonHighEtaAnalyzer::matchId(const HcalDetId& id1, const HcalDetId& id2) {
0923 HcalDetId kd1(id1.subdet(), id1.ieta(), id1.iphi(), 1);
0924 HcalDetId kd2(id1.subdet(), id2.ieta(), id2.iphi(), 1);
0925 int match = ((kd1 == kd2) ? 1 : 0);
0926 return match;
0927 }
0928
0929 double HcalHBHEMuonHighEtaAnalyzer::activeLength(const DetId& hid) {
0930 HcalDetId id(hid);
0931 int ieta = id.ietaAbs();
0932 int zside = id.zside();
0933 int iphi = id.iphi();
0934 std::vector<int> dpths;
0935 if (mergedDepth_) {
0936 std::vector<HcalDetId> ids;
0937 hdc_->unmergeDepthDetId(id, ids);
0938 for (auto idh : ids)
0939 dpths.emplace_back(idh.depth());
0940 } else {
0941 dpths.emplace_back(id.depth());
0942 }
0943 double lx(0);
0944 if (id.subdet() == HcalBarrel) {
0945 for (unsigned int i = 0; i < actHB.size(); ++i) {
0946 if ((ieta == actHB[i].ieta) && (zside == actHB[i].zside) &&
0947 (std::find(dpths.begin(), dpths.end(), actHB[i].depth) != dpths.end()) &&
0948 (std::find(actHB[i].iphis.begin(), actHB[i].iphis.end(), iphi) != actHB[i].iphis.end())) {
0949 lx += actHB[i].thick;
0950 }
0951 }
0952 } else {
0953 for (unsigned int i = 0; i < actHE.size(); ++i) {
0954 if ((ieta == actHE[i].ieta) && (zside == actHE[i].zside) &&
0955 (std::find(dpths.begin(), dpths.end(), actHE[i].depth) != dpths.end()) &&
0956 (std::find(actHE[i].iphis.begin(), actHE[i].iphis.end(), iphi) != actHE[i].iphis.end())) {
0957 lx += actHE[i].thick;
0958 }
0959 }
0960 }
0961 return lx;
0962 }
0963
0964 bool HcalHBHEMuonHighEtaAnalyzer::isGoodVertex(const reco::Vertex& vtx) {
0965 if (vtx.isFake())
0966 return false;
0967 if (vtx.ndof() < 4)
0968 return false;
0969 if (vtx.position().Rho() > 2.)
0970 return false;
0971 if (fabs(vtx.position().Z()) > 24.)
0972 return false;
0973 return true;
0974 }
0975
0976 double HcalHBHEMuonHighEtaAnalyzer::respCorr(const DetId& id) {
0977 double cfac(1.0);
0978 if (useMyCorr_) {
0979 auto itr = corrValue_.find(id);
0980 if (itr != corrValue_.end())
0981 cfac = itr->second;
0982 } else if (respCorrs_ != nullptr) {
0983 cfac = (respCorrs_->getValues(id))->getValue();
0984 }
0985 return cfac;
0986 }
0987
0988 double HcalHBHEMuonHighEtaAnalyzer::gainFactor(const HcalDbService* conditions, const HcalDetId& id) {
0989 double gain(0.0);
0990 const HcalCalibrations& calibs = conditions->getHcalCalibrations(id);
0991 for (int capid = 0; capid < 4; ++capid)
0992 gain += (0.25 * calibs.respcorrgain(capid));
0993 return gain;
0994 }
0995
0996 int HcalHBHEMuonHighEtaAnalyzer::depth16HE(int ieta, int iphi) {
0997
0998
0999
1000 int zside = (ieta > 0) ? 1 : -1;
1001 int depth = theHBHETopology_->dddConstants()->getMinDepth(1, 16, iphi, zside);
1002 if (isItPlan1_ && (!isItPreRecHit_))
1003 depth = 3;
1004 #ifdef EDM_ML_DEBUG
1005 if (verbosity_ % 10 > 0)
1006 edm::LogVerbatim("HBHEMuon") << "Plan1 " << isItPlan1_ << " PreRecHit " << isItPreRecHit_ << " phi " << iphi
1007 << " depth " << depth;
1008 #endif
1009 return depth;
1010 }
1011
1012 bool HcalHBHEMuonHighEtaAnalyzer::goodCell(const HcalDetId& hcid,
1013 const reco::Track* pTrack,
1014 const CaloGeometry* geo,
1015 const MagneticField* bField) {
1016 std::pair<double, double> rz = hdc_->getRZ(hcid);
1017 bool typeRZ = (hcid.subdet() == HcalEndcap) ? false : true;
1018 bool match = spr::propagateHCAL(pTrack, geo, bField, typeRZ, rz, false);
1019 return match;
1020 }
1021
1022 void HcalHBHEMuonHighEtaAnalyzer::fillTrackParameters(const reco::Track* pTrack, math::XYZPoint leadPV) {
1023 trackDz_.emplace_back(pTrack->dz(leadPV));
1024 const reco::HitPattern& hitp = pTrack->hitPattern();
1025 trackLayerCrossed_.emplace_back(hitp.trackerLayersWithMeasurement());
1026 trackOuterHit_.emplace_back(hitp.stripTOBLayersWithMeasurement() + hitp.stripTECLayersWithMeasurement());
1027 trackMissedInnerHits_.emplace_back(hitp.trackerLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1028 trackMissedOuterHits_.emplace_back(hitp.trackerLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1029 }
1030
1031
1032 #include "FWCore/Framework/interface/MakerMacros.h"
1033
1034 DEFINE_FWK_MODULE(HcalHBHEMuonHighEtaAnalyzer);
