Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-28 01:32:02

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 // user include files
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 //////////////trigger info////////////////////////////////////
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 //#define EDM_ML_DEBUG
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   // ----------member data ---------------------------
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   //now do what ever initialization is needed
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 // member functions
0261 //
0262 
0263 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
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", false);
0276   desc.add<bool>("collapseDepth", false);
0277   desc.add<bool>("isItPlan1", false);
0278   desc.add<bool>("getCharge", false);
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 // ------------ method called once each job just before starting event loop  ------------
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 // ------------ method called for each event  ------------
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   // get handles to calogeometry and calotopology
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   // Relevant blocks from iEvent
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   // require a good vertex
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 // ------------ method called when starting to processes a run  ------------
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   // Write correction factors for all HB/HE events
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   //Get track collection
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   ///clearing vectots
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   // Transition between HB/HE is special
0998   // For Run 1 or for Plan1 standard reconstruction it is 3
0999   // For runs beyond 2018 or in Plan1 for HEP17 it is 4
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 //define this as a plug-in
1032 #include "FWCore/Framework/interface/MakerMacros.h"
1033 
1034 DEFINE_FWK_MODULE(HcalHBHEMuonHighEtaAnalyzer);