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