Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-11-28 03:03:21

0001 #include <fstream>
0002 #include <vector>
0003 #include <TTree.h>
0004 
0005 // user include files
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 //////////////trigger info////////////////////////////////////
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 //#define EDM_ML_DEBUG
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   // ----------member data ---------------------------
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   //now do what ever initialization is needed
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 // member functions
0257 //
0258 
0259 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
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 // ------------ method called once each job just before starting event loop  ------------
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 // ------------ method called for each event  ------------
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   // get handles to calogeometry and calotopology
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   // Relevant blocks from iEvent
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   // require a good vertex
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 // ------------ method called when starting to processes a run  ------------
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   // Write correction factors for all HB/HE events
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   //Get track collection
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   ///clearing vectots
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   // Transition between HB/HE is special
0994   // For Run 1 or for Plan1 standard reconstruction it is 3
0995   // For runs beyond 2018 or in Plan1 for HEP17 it is 4
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 //define this as a plug-in
1028 #include "FWCore/Framework/interface/MakerMacros.h"
1029 
1030 DEFINE_FWK_MODULE(HcalHBHEMuonHighEtaAnalyzer);