Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-03 00:58:37

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/stream/EDProducer.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 
0019 #include "DataFormats/HcalCalibObjects/interface/HcalHBHEMuonVariables.h"
0020 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0021 #include "DataFormats/MuonReco/interface/Muon.h"
0022 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0023 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0024 #include "DataFormats/VertexReco/interface/Vertex.h"
0025 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 
0028 //////////////trigger info////////////////////////////////////
0029 
0030 #include "DataFormats/Common/interface/TriggerResults.h"
0031 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0032 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0033 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0034 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0035 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0036 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0037 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0038 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0039 
0040 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0041 #include "HLTrigger/HLTcore/interface/HLTConfigData.h"
0042 
0043 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0044 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0045 
0046 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0047 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0048 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0049 
0050 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0051 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0052 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0053 
0054 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0055 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0056 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0057 
0058 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0059 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0060 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0061 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0062 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0063 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0064 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0065 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0066 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0067 #include "MagneticField/Engine/interface/MagneticField.h"
0068 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0069 
0070 //#define EDM_ML_DEBUG
0071 
0072 namespace alcaHcalHBHEMuon {
0073   struct Counters {
0074     Counters() : nAll_(0), nGood_(0) {}
0075     mutable std::atomic<unsigned int> nAll_, nGood_;
0076   };
0077 }  // namespace alcaHcalHBHEMuon
0078 
0079 class AlCaHcalHBHEMuonProducer : public edm::stream::EDProducer<edm::GlobalCache<alcaHcalHBHEMuon::Counters>> {
0080 public:
0081   explicit AlCaHcalHBHEMuonProducer(const edm::ParameterSet&, const alcaHcalHBHEMuon::Counters*);
0082   ~AlCaHcalHBHEMuonProducer() override = default;
0083 
0084   static std::unique_ptr<alcaHcalHBHEMuon::Counters> initializeGlobalCache(edm::ParameterSet const&) {
0085     return std::make_unique<alcaHcalHBHEMuon::Counters>();
0086   }
0087 
0088   void produce(edm::Event&, const edm::EventSetup&) override;
0089 
0090   void endStream() override;
0091 
0092   static void globalEndJob(const alcaHcalHBHEMuon::Counters* counters);
0093 
0094   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0095 
0096 private:
0097   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0098   void endRun(edm::Run const&, edm::EventSetup const&) override;
0099   int matchId(const HcalDetId&, const HcalDetId&);
0100   double activeLength(const DetId&, const HcalDDDRecConstants* hdc);
0101   bool isGoodVertex(const reco::Vertex& vtx);
0102   double respCorr(const DetId& id, const HcalRespCorrs* respCorrs);
0103   double gainFactor(const HcalDbService* dbserv, const HcalDetId& id);
0104   int depth16HE(int ieta, int iphi, const HcalTopology* theHBHETopology);
0105   bool goodCell(const HcalDetId& hcid,
0106                 const reco::Track* pTrack,
0107                 const CaloGeometry* geo,
0108                 const MagneticField* bField,
0109                 const HcalDDDRecConstants* hdc);
0110 
0111   // ----------member data ---------------------------
0112   HLTConfigProvider hltConfig_;
0113   const std::vector<std::string> trigNames_;
0114   const std::string processName_;
0115   const edm::InputTag triggerResults_;
0116   const edm::InputTag labelEBRecHit_, labelEERecHit_, labelHBHERecHit_;
0117   const std::string labelVtx_, labelMuon_, labelHBHEMuon_;
0118   const bool collapseDepth_, isItPlan1_;
0119   const int verbosity_;
0120   const bool isItPreRecHit_, writeRespCorr_;
0121   const std::string fileInCorr_;
0122   const int maxDepth_;
0123   const bool mergedDepth_;
0124 
0125   bool useMyCorr_;
0126   int nRun_, nAll_, nGood_;
0127 
0128   const edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0129   const edm::EDGetTokenT<reco::VertexCollection> tok_Vtx_;
0130   const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0131   const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0132   const edm::EDGetTokenT<HBHERecHitCollection> tok_HBHE_;
0133   const edm::EDGetTokenT<reco::MuonCollection> tok_Muon_;
0134 
0135   const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec0_;
0136   const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec1_;
0137   const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_respcorr0_;
0138   const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_respcorr1_;
0139   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom0_;
0140   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom1_;
0141   const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo0_;
0142   const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo1_;
0143   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0144   const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_chan_;
0145   const edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0146   const edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_topo_;
0147   const edm::ESGetToken<HcalDbService, HcalDbRecord> tok_dbservice_;
0148 
0149   //////////////////////////////////////////////////////
0150   static const int depthMax_ = 7;
0151   std::vector<std::string> all_triggers_;
0152   std::vector<int> hltresults_;
0153 
0154   std::vector<HcalDDDRecConstants::HcalActiveLength> actHB, actHE;
0155   std::map<DetId, double> corrValue_;
0156   ////////////////////////////////////////////////////////////
0157 };
0158 
0159 AlCaHcalHBHEMuonProducer::AlCaHcalHBHEMuonProducer(const edm::ParameterSet& iConfig, const alcaHcalHBHEMuon::Counters*)
0160     : trigNames_(iConfig.getParameter<std::vector<std::string>>("triggers")),
0161       processName_(iConfig.getParameter<std::string>("processName")),
0162       triggerResults_(iConfig.getParameter<edm::InputTag>("triggerResults")),
0163       labelEBRecHit_(iConfig.getParameter<edm::InputTag>("labelEBRecHit")),
0164       labelEERecHit_(iConfig.getParameter<edm::InputTag>("labelEERecHit")),
0165       labelHBHERecHit_(iConfig.getParameter<edm::InputTag>("labelHBHERecHit")),
0166       labelVtx_(iConfig.getParameter<std::string>("labelVertex")),
0167       labelMuon_(iConfig.getParameter<std::string>("labelMuon")),
0168       labelHBHEMuon_(iConfig.getParameter<std::string>("labelHBHEMuon")),
0169       collapseDepth_(iConfig.getParameter<bool>("collapseDepth")),
0170       isItPlan1_(iConfig.getParameter<bool>("isItPlan1")),
0171       verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0172       isItPreRecHit_(iConfig.getUntrackedParameter<bool>("isItPreRecHit", false)),
0173       writeRespCorr_(iConfig.getUntrackedParameter<bool>("writeRespCorr", false)),
0174       fileInCorr_(iConfig.getUntrackedParameter<std::string>("fileInCorr", "")),
0175       maxDepth_(iConfig.getUntrackedParameter<int>("maxDepth", 4)),
0176       mergedDepth_((!isItPreRecHit_) || (collapseDepth_)),
0177       nRun_(0),
0178       nAll_(0),
0179       nGood_(0),
0180       tok_trigRes_(consumes<edm::TriggerResults>(triggerResults_)),
0181       tok_Vtx_(consumes<reco::VertexCollection>(labelVtx_)),
0182       tok_EB_(consumes<EcalRecHitCollection>(labelEBRecHit_)),
0183       tok_EE_(consumes<EcalRecHitCollection>(labelEERecHit_)),
0184       tok_HBHE_(consumes<HBHERecHitCollection>(labelHBHERecHit_)),
0185       tok_Muon_(consumes<reco::MuonCollection>(labelMuon_)),
0186       tok_ddrec0_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0187       tok_ddrec1_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>()),
0188       tok_respcorr0_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd, edm::Transition::BeginRun>()),
0189       tok_respcorr1_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd>()),
0190       tok_geom0_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0191       tok_geom1_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0192       tok_htopo0_(esConsumes<HcalTopology, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0193       tok_htopo1_(esConsumes<HcalTopology, HcalRecNumberingRecord>()),
0194       tok_magField_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0195       tok_chan_(esConsumes<EcalChannelStatus, EcalChannelStatusRcd>()),
0196       tok_sevlv_(esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>()),
0197       tok_topo_(esConsumes<CaloTopology, CaloTopologyRecord>()),
0198       tok_dbservice_(esConsumes<HcalDbService, HcalDbRecord>()) {
0199   //now do what ever initialization is needed
0200   edm::LogVerbatim("HBHEMuon") << "Labels used: Trig " << triggerResults_ << " Vtx " << labelVtx_ << " EB "
0201                                << labelEBRecHit_ << " EE " << labelEERecHit_ << " HBHE " << labelHBHERecHit_ << " MU "
0202                                << labelMuon_;
0203 
0204   if (!fileInCorr_.empty()) {
0205     std::ifstream infile(fileInCorr_.c_str());
0206     if (infile.is_open()) {
0207       while (true) {
0208         unsigned int id;
0209         double cfac;
0210         infile >> id >> cfac;
0211         if (!infile.good())
0212           break;
0213         corrValue_[DetId(id)] = cfac;
0214       }
0215       infile.close();
0216     }
0217   }
0218   useMyCorr_ = (!corrValue_.empty());
0219   edm::LogVerbatim("HBHEMuon") << "Flags used: ollapseDepth " << collapseDepth_ << ":" << mergedDepth_ << " IsItPlan1 "
0220                                << isItPlan1_ << " IsItPreRecHit " << isItPreRecHit_ << " UseMyCorr " << useMyCorr_;
0221 
0222   //create the objects for HcalHBHEMuonVariables which has information of isolated muons
0223   produces<HcalHBHEMuonVariablesCollection>(labelHBHEMuon_);
0224   edm::LogVerbatim("HcalIsoTrack") << " Expected to produce the collections:\n"
0225                                    << "HcalHBHEMuonVariablesCollection with label " << labelHBHEMuon_;
0226 }
0227 
0228 //
0229 // member functions
0230 //
0231 
0232 // ------------ method called for each event  ------------
0233 void AlCaHcalHBHEMuonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0234   ++nAll_;
0235   auto outputHcalHBHEMuonColl = std::make_unique<HcalHBHEMuonVariablesCollection>();
0236 
0237   unsigned int runNumber = iEvent.id().run();
0238   unsigned int eventNumber = iEvent.id().event();
0239   unsigned int lumiNumber = iEvent.id().luminosityBlock();
0240   unsigned int bxNumber = iEvent.bunchCrossing();
0241 #ifdef EDM_ML_DEBUG
0242   edm::LogVerbatim("HBHEMuon") << "Run " << runNumber << " Event " << eventNumber << " Lumi " << lumiNumber << " BX "
0243                                << bxNumber << std::endl;
0244 #endif
0245 
0246   //Step1: Find if the event passes one of the chosen triggers
0247   bool ok(false);
0248   /////////////////////////////TriggerResults
0249   if (trigNames_.empty()) {
0250     ok = true;
0251   } else {
0252     auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
0253     if (triggerResults.isValid()) {
0254       std::vector<std::string> modules;
0255       const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0256       const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
0257       for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0258         int hlt = triggerResults->accept(iHLT);
0259         for (unsigned int i = 0; i < trigNames_.size(); ++i) {
0260           if (triggerNames_[iHLT].find(trigNames_[i]) != std::string::npos) {
0261             if (hlt > 0) {
0262               ok = true;
0263             }
0264 #ifdef EDM_ML_DEBUG
0265             edm::LogVerbatim("HBHEMuon") << "AlCaHcalHBHEMuonFilter::Trigger " << triggerNames_[iHLT] << " Flag " << hlt
0266                                          << ":" << ok << std::endl;
0267 #endif
0268           }
0269         }
0270       }
0271     }
0272   }
0273 
0274   // get handles to calogeometry and calotopology
0275   const HcalDDDRecConstants* hdc = &iSetup.getData(tok_ddrec1_);
0276   const HcalRespCorrs* resp = &iSetup.getData(tok_respcorr1_);
0277   const CaloGeometry* geo = &iSetup.getData(tok_geom1_);
0278   const HcalTopology* theHBHETopology = &iSetup.getData(tok_htopo1_);
0279   const MagneticField* bField = &iSetup.getData(tok_magField_);
0280   const EcalChannelStatus* theEcalChStatus = &iSetup.getData(tok_chan_);
0281   const EcalSeverityLevelAlgo* sevlv = &iSetup.getData(tok_sevlv_);
0282   const CaloTopology* caloTopology = &iSetup.getData(tok_topo_);
0283   const HcalDbService* conditions = &iSetup.getData(tok_dbservice_);
0284   HcalRespCorrs respCorrsObj(*resp);
0285   HcalRespCorrs* respCorrs = &respCorrsObj;
0286   respCorrs->setTopo(theHBHETopology);
0287 
0288   // Relevant blocks from iEvent
0289   auto const& vtx = iEvent.getHandle(tok_Vtx_);
0290   auto barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0291   auto endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0292   auto hbhe = iEvent.getHandle(tok_HBHE_);
0293   auto const& muons = iEvent.getHandle(tok_Muon_);
0294 
0295   // require a good vertex
0296   math::XYZPoint pvx;
0297   unsigned int goodVertex = 0;
0298   reco::VertexCollection::const_iterator firstGoodVertex;
0299   if (!vtx.isValid()) {
0300 #ifdef EDM_ML_DEBUG
0301     edm::LogVerbatim("HBHEMuon") << "No Good Vertex found == Reject";
0302 #endif
0303   } else {
0304     firstGoodVertex = vtx->end();
0305     for (reco::VertexCollection::const_iterator it = vtx->begin(); it != vtx->end(); it++) {
0306       if (isGoodVertex(*it)) {
0307         if (firstGoodVertex == vtx->end())
0308           firstGoodVertex = it;
0309         ++goodVertex;
0310       }
0311     }
0312     if (firstGoodVertex != vtx->end())
0313       pvx = firstGoodVertex->position();
0314   }
0315 
0316   if (ok && (goodVertex > 0) && muons.isValid() && barrelRecHitsHandle.isValid() && endcapRecHitsHandle.isValid() &&
0317       hbhe.isValid()) {
0318     for (reco::MuonCollection::const_iterator recMuon = muons->begin(); recMuon != muons->end(); ++recMuon) {
0319       HcalHBHEMuonVariables hbheMuon;
0320       hbheMuon.runNumber_ = runNumber;
0321       hbheMuon.eventNumber_ = eventNumber;
0322       hbheMuon.lumiNumber_ = lumiNumber;
0323       hbheMuon.bxNumber_ = bxNumber;
0324       hbheMuon.goodVertex_ = goodVertex;
0325       hbheMuon.muonGood_ = (recMuon->isPFMuon());
0326       hbheMuon.muonGlobal_ = (recMuon->isGlobalMuon());
0327       hbheMuon.muonTracker_ = (recMuon->isTrackerMuon());
0328       hbheMuon.ptGlob_ = ((recMuon)->pt());
0329       hbheMuon.etaGlob_ = (recMuon->eta());
0330       hbheMuon.phiGlob_ = (recMuon->phi());
0331       hbheMuon.energyMuon_ = (recMuon->energy());
0332       hbheMuon.pMuon_ = (recMuon->p());
0333 #ifdef EDM_ML_DEBUG
0334       edm::LogVerbatim("HBHEMuon") << "Energy:" << recMuon->energy() << " P:" << recMuon->p();
0335 #endif
0336       hbheMuon.muonTight_ = (muon::isTightMuon(*recMuon, *firstGoodVertex));
0337       hbheMuon.muonMedium_ = (muon::isMediumMuon(*recMuon));
0338       hbheMuon.muonTrkKink_ = (recMuon->combinedQuality().trkKink);
0339       hbheMuon.muonChi2LocalPosition_ = (recMuon->combinedQuality().chi2LocalPosition);
0340       hbheMuon.muonSegComp_ = (muon::segmentCompatibility(*recMuon));
0341       // acessing tracker hits info
0342       if (recMuon->track().isNonnull()) {
0343         hbheMuon.trackerLayer_ = (recMuon->track()->hitPattern().trackerLayersWithMeasurement());
0344       } else {
0345         hbheMuon.trackerLayer_ = -1;
0346       }
0347       if (recMuon->innerTrack().isNonnull()) {
0348         hbheMuon.innerTrack_ = true;
0349         hbheMuon.numPixelLayers_ = (recMuon->innerTrack()->hitPattern().pixelLayersWithMeasurement());
0350         hbheMuon.chiTracker_ = (recMuon->innerTrack()->normalizedChi2());
0351         hbheMuon.dxyTracker_ = (fabs(recMuon->innerTrack()->dxy(pvx)));
0352         hbheMuon.dzTracker_ = (fabs(recMuon->innerTrack()->dz(pvx)));
0353         hbheMuon.innerTrackPt_ = (recMuon->innerTrack()->pt());
0354         hbheMuon.innerTrackEta_ = (recMuon->innerTrack()->eta());
0355         hbheMuon.innerTrackPhi_ = (recMuon->innerTrack()->phi());
0356         hbheMuon.tightPixelHits_ = (recMuon->innerTrack()->hitPattern().numberOfValidPixelHits());
0357         hbheMuon.tightValidFraction_ = (recMuon->innerTrack()->validFraction());
0358       }
0359       // outer track info
0360       if (recMuon->outerTrack().isNonnull()) {
0361         hbheMuon.outerTrack_ = true;
0362         hbheMuon.outerTrackPt_ = (recMuon->outerTrack()->pt());
0363         hbheMuon.outerTrackEta_ = (recMuon->outerTrack()->eta());
0364         hbheMuon.outerTrackPhi_ = (recMuon->outerTrack()->phi());
0365         hbheMuon.outerTrackChi_ = (recMuon->outerTrack()->normalizedChi2());
0366         hbheMuon.outerTrackHits_ = (recMuon->outerTrack()->numberOfValidHits());
0367         hbheMuon.outerTrackRHits_ = (recMuon->outerTrack()->recHitsSize());
0368       }
0369       // Tight Muon cuts
0370       if (recMuon->globalTrack().isNonnull()) {
0371         hbheMuon.globalTrack_ = true;
0372         hbheMuon.chiGlobal_ = (recMuon->globalTrack()->normalizedChi2());
0373         hbheMuon.globalMuonHits_ = (recMuon->globalTrack()->hitPattern().numberOfValidMuonHits());
0374         hbheMuon.matchedStat_ = (recMuon->numberOfMatchedStations());
0375         hbheMuon.globalTrackPt_ = (recMuon->globalTrack()->pt());
0376         hbheMuon.globalTrackEta_ = (recMuon->globalTrack()->eta());
0377         hbheMuon.globalTrackPhi_ = (recMuon->globalTrack()->phi());
0378         hbheMuon.tightTransImpara_ = (fabs(recMuon->muonBestTrack()->dxy(pvx)));
0379         hbheMuon.tightLongPara_ = (fabs(recMuon->muonBestTrack()->dz(pvx)));
0380       }
0381 
0382       hbheMuon.isolationR04_ =
0383           ((recMuon->pfIsolationR04().sumChargedHadronPt +
0384             std::max(0.,
0385                      recMuon->pfIsolationR04().sumNeutralHadronEt + recMuon->pfIsolationR04().sumPhotonEt -
0386                          (0.5 * recMuon->pfIsolationR04().sumPUPt))) /
0387            recMuon->pt());
0388 
0389       hbheMuon.isolationR03_ =
0390           ((recMuon->pfIsolationR03().sumChargedHadronPt +
0391             std::max(0.,
0392                      recMuon->pfIsolationR03().sumNeutralHadronEt + recMuon->pfIsolationR03().sumPhotonEt -
0393                          (0.5 * recMuon->pfIsolationR03().sumPUPt))) /
0394            recMuon->pt());
0395 
0396       hbheMuon.ecalEnergy_ = (recMuon->calEnergy().emS9);
0397       hbheMuon.hcalEnergy_ = (recMuon->calEnergy().hadS9);
0398       hbheMuon.hoEnergy_ = (recMuon->calEnergy().hoS9);
0399 
0400       if (recMuon->innerTrack().isNonnull()) {
0401         const reco::Track* pTrack = (recMuon->innerTrack()).get();
0402         spr::propagatedTrackID trackID = spr::propagateCALO(pTrack, geo, bField, (((verbosity_ / 100) % 10 > 0)));
0403 
0404         double activeLengthTot(0), activeLengthHotTot(0);
0405         double eHcalDepth[depthMax_], eHcalDepthHot[depthMax_];
0406         double eHcalDepthC[depthMax_], eHcalDepthHotC[depthMax_];
0407         double cHcalDepthHot[depthMax_], cHcalDepthHotBG[depthMax_];
0408         double eHcalDepthRaw[depthMax_], eHcalDepthHotRaw[depthMax_];
0409         double eHcalDepthCRaw[depthMax_], eHcalDepthHotCRaw[depthMax_];
0410         double cHcalDepthHotRaw[depthMax_], cHcalDepthHotBGRaw[depthMax_];
0411         double eHcalDepthAux[depthMax_], eHcalDepthHotAux[depthMax_];
0412         double eHcalDepthCAux[depthMax_], eHcalDepthHotCAux[depthMax_];
0413         double cHcalDepthHotAux[depthMax_], cHcalDepthHotBGAux[depthMax_];
0414         double activeL[depthMax_], activeHotL[depthMax_];
0415         bool matchDepth[depthMax_], matchDepthHot[depthMax_];
0416         HcalDetId eHcalDetId[depthMax_];
0417         unsigned int isHot(0);
0418         int ieta(-1000), iphi(-1000);
0419         for (int i = 0; i < depthMax_; ++i) {
0420           eHcalDepth[i] = eHcalDepthHot[i] = 0;
0421           eHcalDepthC[i] = eHcalDepthHotC[i] = 0;
0422           cHcalDepthHot[i] = cHcalDepthHotBG[i] = 0;
0423           eHcalDepthRaw[i] = eHcalDepthHotRaw[i] = 0;
0424           eHcalDepthCRaw[i] = eHcalDepthHotCRaw[i] = 0;
0425           cHcalDepthHotRaw[i] = cHcalDepthHotBGRaw[i] = 0;
0426           eHcalDepthAux[i] = eHcalDepthHotAux[i] = 0;
0427           eHcalDepthCAux[i] = eHcalDepthHotCAux[i] = 0;
0428           cHcalDepthHotAux[i] = cHcalDepthHotBGAux[i] = 0;
0429           activeL[i] = activeHotL[i] = 0;
0430           matchDepth[i] = matchDepthHot[i] = true;
0431         }
0432 #ifdef EDM_ML_DEBUG
0433         double eHcal(0);
0434 #endif
0435 
0436         hbheMuon.ecalDetId_ = ((trackID.detIdECAL)());
0437         hbheMuon.hcalDetId_ = ((trackID.detIdHCAL)());
0438         hbheMuon.ehcalDetId_ = ((trackID.detIdEHCAL)());
0439 
0440         HcalDetId check(false);
0441         std::pair<bool, HcalDetId> info = spr::propagateHCALBack(pTrack, geo, bField, (((verbosity_ / 100) % 10 > 0)));
0442         if (info.first) {
0443           check = info.second;
0444         }
0445 
0446         bool okE = trackID.okECAL;
0447         if (okE) {
0448           const DetId isoCell(trackID.detIdECAL);
0449           std::pair<double, bool> e3x3 = spr::eECALmatrix(isoCell,
0450                                                           barrelRecHitsHandle,
0451                                                           endcapRecHitsHandle,
0452                                                           *theEcalChStatus,
0453                                                           geo,
0454                                                           caloTopology,
0455                                                           sevlv,
0456                                                           1,
0457                                                           1,
0458                                                           -100.0,
0459                                                           -100.0,
0460                                                           -500.0,
0461                                                           500.0,
0462                                                           false);
0463           hbheMuon.ecal3x3Energy_ = e3x3.first;
0464           okE = e3x3.second;
0465         }
0466 #ifdef EDM_ML_DEBUG
0467         edm::LogVerbatim("HBHEMuon") << "Propagate Track to ECAL: " << okE << ":" << trackID.okECAL << " E "
0468                                      << hbheMuon.ecal3x3Energy_;
0469 #endif
0470 
0471         if (trackID.okHCAL) {
0472           DetId closestCell(trackID.detIdHCAL);
0473           HcalDetId hcidt(closestCell.rawId());
0474           if ((hcidt.ieta() == check.ieta()) && (hcidt.iphi() == check.iphi()))
0475             hbheMuon.matchedId_ = true;
0476 #ifdef EDM_ML_DEBUG
0477           edm::LogVerbatim("HBHEMuon") << "Front " << hcidt << " Back " << info.first << ":" << check << " Match "
0478                                        << hbheMuon.matchedId_;
0479 #endif
0480           HcalSubdetector subdet = hcidt.subdet();
0481           ieta = hcidt.ieta();
0482           iphi = hcidt.iphi();
0483           bool hborhe = (std::abs(ieta) == 16);
0484 
0485           hbheMuon.hcal1x1Energy_ = spr::eHCALmatrix(
0486               theHBHETopology, closestCell, hbhe, 0, 0, false, true, -100.0, -100.0, -100.0, -100.0, -500., 500., 0);
0487           hbheMuon.hcal1x1EnergyAux_ = spr::eHCALmatrix(
0488               theHBHETopology, closestCell, hbhe, 0, 0, false, true, -100.0, -100.0, -100.0, -100.0, -500., 500., 1);
0489           hbheMuon.hcal1x1EnergyRaw_ = spr::eHCALmatrix(
0490               theHBHETopology, closestCell, hbhe, 0, 0, false, true, -100.0, -100.0, -100.0, -100.0, -500., 500., 2);
0491           std::vector<std::pair<double, int>> ehdepth, ehdepthAux, ehdepthRaw;
0492           spr::energyHCALCell((HcalDetId)closestCell,
0493                               hbhe,
0494                               ehdepth,
0495                               maxDepth_,
0496                               -100.0,
0497                               -100.0,
0498                               -100.0,
0499                               -100.0,
0500                               -500.0,
0501                               500.0,
0502                               0,
0503                               depth16HE(ieta, iphi, theHBHETopology),
0504                               (((verbosity_ / 1000) % 10) > 0));
0505           for (int i = 0; i < depthMax_; ++i)
0506             eHcalDetId[i] = HcalDetId();
0507           spr::energyHCALCell((HcalDetId)closestCell,
0508                               hbhe,
0509                               ehdepthAux,
0510                               maxDepth_,
0511                               -100.0,
0512                               -100.0,
0513                               -100.0,
0514                               -100.0,
0515                               -500.0,
0516                               500.0,
0517                               1,
0518                               depth16HE(ieta, iphi, theHBHETopology),
0519                               (((verbosity_ / 1000) % 10) > 0));
0520           spr::energyHCALCell((HcalDetId)closestCell,
0521                               hbhe,
0522                               ehdepthRaw,
0523                               maxDepth_,
0524                               -100.0,
0525                               -100.0,
0526                               -100.0,
0527                               -100.0,
0528                               -500.0,
0529                               500.0,
0530                               2,
0531                               depth16HE(ieta, iphi, theHBHETopology),
0532                               (((verbosity_ / 1000) % 10) > 0));
0533           for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0534             HcalSubdetector subdet0 =
0535                 (hborhe) ? ((ehdepth[i].second >= depth16HE(ieta, iphi, theHBHETopology)) ? HcalEndcap : HcalBarrel)
0536                          : subdet;
0537             HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
0538             double actL = activeLength(DetId(hcid0), hdc);
0539             double ene = ehdepth[i].first;
0540             double eneAux = ehdepthAux[i].first;
0541             double eneRaw = ehdepthRaw[i].first;
0542             bool tmpC(false);
0543             if (ene > 0.0) {
0544               if (!(theHBHETopology->validHcal(hcid0))) {
0545                 edm::LogWarning("HBHEMuon") << "(1) Invalid ID " << hcid0 << " with E = " << ene;
0546                 edm::LogWarning("HBHEMuon") << HcalDetId(closestCell) << " with " << ehdepth.size() << " depths:";
0547                 for (const auto& ehd : ehdepth)
0548                   edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0549               } else {
0550                 tmpC = goodCell(hcid0, pTrack, geo, bField, hdc);
0551                 double enec(ene);
0552                 double corr = respCorr(DetId(hcid0), respCorrs);
0553                 if (corr != 0)
0554                   ene /= corr;
0555 #ifdef EDM_ML_DEBUG
0556                 HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc->mergedDepthDetId(hcid0) : hcid0;
0557                 edm::LogVerbatim("HBHEMuon") << hcid0 << ":" << id << " Corr " << corr;
0558 #endif
0559                 int depth = ehdepth[i].second - 1;
0560                 if (collapseDepth_) {
0561                   HcalDetId id = hdc->mergedDepthDetId(hcid0);
0562                   depth = id.depth() - 1;
0563                 }
0564                 eHcalDepth[depth] += ene;
0565                 eHcalDepthC[depth] += enec;
0566                 activeL[depth] += actL;
0567                 activeLengthTot += actL;
0568                 matchDepth[depth] = (matchDepth[depth] && tmpC);
0569 #ifdef EDM_ML_DEBUG
0570                 if ((verbosity_ % 10) > 0)
0571                   edm::LogVerbatim("HBHEMuon")
0572                       << hcid0 << " E " << ene << ":" << enec << " L " << actL << " Match " << tmpC;
0573 #endif
0574               }
0575             }
0576             if (eneAux > 0.0) {
0577               if (theHBHETopology->validHcal(hcid0)) {
0578                 double enecAux(eneAux);
0579                 double corr = respCorr(DetId(hcid0), respCorrs);
0580                 if (corr != 0)
0581                   eneAux /= corr;
0582                 int depth = ehdepthAux[i].second - 1;
0583                 if (collapseDepth_) {
0584                   HcalDetId id = hdc->mergedDepthDetId(hcid0);
0585                   depth = id.depth() - 1;
0586                 }
0587                 eHcalDepthAux[depth] += eneAux;
0588                 eHcalDepthCAux[depth] += enecAux;
0589 #ifdef EDM_ML_DEBUG
0590                 if ((verbosity_ % 10) > 0)
0591                   edm::LogVerbatim("HBHEMuon")
0592                       << hcid0 << " E " << eneAux << ":" << enecAux << " L " << actL << " Match " << tmpC;
0593 #endif
0594               }
0595             }
0596             if (eneRaw > 0.0) {
0597               if (theHBHETopology->validHcal(hcid0)) {
0598                 double enecRaw(eneRaw);
0599                 double corr = respCorr(DetId(hcid0), respCorrs);
0600                 if (corr != 0)
0601                   eneRaw /= corr;
0602                 int depth = ehdepthRaw[i].second - 1;
0603                 if (collapseDepth_) {
0604                   HcalDetId id = hdc->mergedDepthDetId(hcid0);
0605                   depth = id.depth() - 1;
0606                 }
0607                 eHcalDepthRaw[depth] += eneRaw;
0608                 eHcalDepthCRaw[depth] += enecRaw;
0609 #ifdef EDM_ML_DEBUG
0610                 if ((verbosity_ % 10) > 0)
0611                   edm::LogVerbatim("HBHEMuon")
0612                       << hcid0 << " E " << eneRaw << ":" << enecRaw << " L " << actL << " Match " << tmpC;
0613 #endif
0614               }
0615             }
0616           }
0617 #ifdef EDM_ML_DEBUG
0618           if ((verbosity_ % 10) > 0) {
0619             edm::LogVerbatim("HBHEMuon") << hcidt << " Match " << hbheMuon.matchedId_ << " Depths " << ehdepth.size();
0620             for (unsigned int k = 0; k < ehdepth.size(); ++k)
0621               edm::LogVerbatim("HBHEMuon") << " [" << k << ":" << ehdepth[k].second << "] " << matchDepth[k];
0622           }
0623 #endif
0624           HcalDetId hotCell;
0625           spr::eHCALmatrix(geo, theHBHETopology, closestCell, hbhe, 1, 1, hotCell, false, 0, false);
0626           isHot = matchId(closestCell, hotCell);
0627           if (hotCell != HcalDetId()) {
0628             subdet = HcalDetId(hotCell).subdet();
0629             ieta = HcalDetId(hotCell).ieta();
0630             iphi = HcalDetId(hotCell).iphi();
0631             hborhe = (std::abs(ieta) == 16);
0632             std::vector<std::pair<double, int>> ehdepth, ehdepthAux, ehdepthRaw;
0633             spr::energyHCALCell(hotCell,
0634                                 hbhe,
0635                                 ehdepth,
0636                                 maxDepth_,
0637                                 -100.0,
0638                                 -100.0,
0639                                 -100.0,
0640                                 -100.0,
0641                                 -500.0,
0642                                 500.0,
0643                                 0,
0644                                 depth16HE(ieta, iphi, theHBHETopology),
0645                                 false);
0646             spr::energyHCALCell(hotCell,
0647                                 hbhe,
0648                                 ehdepthAux,
0649                                 maxDepth_,
0650                                 -100.0,
0651                                 -100.0,
0652                                 -100.0,
0653                                 -100.0,
0654                                 -500.0,
0655                                 500.0,
0656                                 1,
0657                                 depth16HE(ieta, iphi, theHBHETopology),
0658                                 false);
0659             spr::energyHCALCell(hotCell,
0660                                 hbhe,
0661                                 ehdepthRaw,
0662                                 maxDepth_,
0663                                 -100.0,
0664                                 -100.0,
0665                                 -100.0,
0666                                 -100.0,
0667                                 -500.0,
0668                                 500.0,
0669                                 2,
0670                                 depth16HE(ieta, iphi, theHBHETopology),
0671                                 false);
0672             for (int i = 0; i < depthMax_; ++i)
0673               eHcalDetId[i] = HcalDetId();
0674             for (unsigned int i = 0; i < ehdepth.size(); ++i) {
0675               HcalSubdetector subdet0 =
0676                   (hborhe) ? ((ehdepth[i].second >= depth16HE(ieta, iphi, theHBHETopology)) ? HcalEndcap : HcalBarrel)
0677                            : subdet;
0678               HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
0679               double actL = activeLength(DetId(hcid0), hdc);
0680               double ene = ehdepth[i].first;
0681               bool tmpC(false);
0682               if (ene > 0.0) {
0683                 if (!(theHBHETopology->validHcal(hcid0))) {
0684                   edm::LogWarning("HBHEMuon") << "(2) Invalid ID " << hcid0 << " with E = " << ene;
0685                   edm::LogWarning("HBHEMuon") << HcalDetId(hotCell) << " with " << ehdepth.size() << " depths:";
0686                   for (const auto& ehd : ehdepth)
0687                     edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0688                 } else {
0689                   tmpC = goodCell(hcid0, pTrack, geo, bField, hdc);
0690                   double chg(ene), enec(ene);
0691                   double corr = respCorr(DetId(hcid0), respCorrs);
0692                   if (corr != 0)
0693                     ene /= corr;
0694 #ifdef EDM_ML_DEBUG
0695                   HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc->mergedDepthDetId(hcid0) : hcid0;
0696                   edm::LogVerbatim("HBHEMuon") << hcid0 << ":" << id << " Corr " << corr << " E " << ene << ":" << enec;
0697 #endif
0698                   double gain = gainFactor(conditions, hcid0);
0699                   if (gain != 0)
0700                     chg /= gain;
0701 #ifdef EDM_ML_DEBUG
0702                   edm::LogVerbatim("HBHEMuon") << hcid0 << " Gain " << gain << " C " << chg;
0703 #endif
0704                   int depth = ehdepth[i].second - 1;
0705                   if (collapseDepth_) {
0706                     HcalDetId id = hdc->mergedDepthDetId(hcid0);
0707                     depth = id.depth() - 1;
0708                   }
0709                   eHcalDepthHot[depth] += ene;
0710                   eHcalDepthHotC[depth] += enec;
0711                   cHcalDepthHot[depth] += chg;
0712                   activeHotL[depth] += actL;
0713                   activeLengthHotTot += actL;
0714                   matchDepthHot[depth] = (matchDepthHot[depth] && tmpC);
0715 #ifdef EDM_ML_DEBUG
0716                   eHcal += ene;
0717                   if ((verbosity_ % 10) > 0)
0718                     edm::LogVerbatim("HBHEMuon") << hcid0 << " depth " << depth << " E " << ene << ":" << enec << " C "
0719                                                  << chg << " L " << actL << " Match " << tmpC;
0720 #endif
0721                 }
0722               }
0723               double eneAux = ehdepthAux[i].first;
0724               if (eneAux > 0.0) {
0725                 if (theHBHETopology->validHcal(hcid0)) {
0726                   double chgAux(eneAux), enecAux(eneAux);
0727                   double corr = respCorr(DetId(hcid0), respCorrs);
0728                   if (corr != 0)
0729                     eneAux /= corr;
0730                   double gain = gainFactor(conditions, hcid0);
0731                   if (gain != 0)
0732                     chgAux /= gain;
0733                   int depth = ehdepthAux[i].second - 1;
0734                   if (collapseDepth_) {
0735                     HcalDetId id = hdc->mergedDepthDetId(hcid0);
0736                     depth = id.depth() - 1;
0737                   }
0738                   eHcalDepthHotAux[depth] += eneAux;
0739                   eHcalDepthHotCAux[depth] += enecAux;
0740                   cHcalDepthHotAux[depth] += chgAux;
0741 #ifdef EDM_ML_DEBUG
0742                   if ((verbosity_ % 10) > 0)
0743                     edm::LogVerbatim("HBHEMuon") << hcid0 << " depth " << depth << " E " << eneAux << ":" << enecAux
0744                                                  << " C " << chgAux << " L " << actL << " Match " << tmpC;
0745 #endif
0746                 }
0747               }
0748               double eneRaw = ehdepthRaw[i].first;
0749               if (eneRaw > 0.0) {
0750                 if (theHBHETopology->validHcal(hcid0)) {
0751                   double chgRaw(eneRaw), enecRaw(eneRaw);
0752                   double corr = respCorr(DetId(hcid0), respCorrs);
0753                   if (corr != 0)
0754                     eneRaw /= corr;
0755                   double gain = gainFactor(conditions, hcid0);
0756                   if (gain != 0)
0757                     chgRaw /= gain;
0758                   int depth = ehdepthRaw[i].second - 1;
0759                   if (collapseDepth_) {
0760                     HcalDetId id = hdc->mergedDepthDetId(hcid0);
0761                     depth = id.depth() - 1;
0762                   }
0763                   eHcalDepthHotRaw[depth] += eneRaw;
0764                   eHcalDepthHotCRaw[depth] += enecRaw;
0765                   cHcalDepthHotRaw[depth] += chgRaw;
0766 #ifdef EDM_ML_DEBUG
0767                   if ((verbosity_ % 10) > 0)
0768                     edm::LogVerbatim("HBHEMuon") << hcid0 << " depth " << depth << " E " << eneRaw << ":" << enecRaw
0769                                                  << " C " << chgRaw << " L " << actL << " Match " << tmpC;
0770 #endif
0771                 }
0772               }
0773             }
0774 
0775             HcalDetId oppCell(subdet, -ieta, iphi, HcalDetId(hotCell).depth());
0776             std::vector<std::pair<double, int>> ehdeptho, ehdepthoAux, ehdepthoRaw;
0777             spr::energyHCALCell(oppCell,
0778                                 hbhe,
0779                                 ehdeptho,
0780                                 maxDepth_,
0781                                 -100.0,
0782                                 -100.0,
0783                                 -100.0,
0784                                 -100.0,
0785                                 -500.0,
0786                                 500.0,
0787                                 0,
0788                                 depth16HE(-ieta, iphi, theHBHETopology),
0789                                 false);
0790             spr::energyHCALCell(oppCell,
0791                                 hbhe,
0792                                 ehdepthoAux,
0793                                 maxDepth_,
0794                                 -100.0,
0795                                 -100.0,
0796                                 -100.0,
0797                                 -100.0,
0798                                 -500.0,
0799                                 500.0,
0800                                 1,
0801                                 depth16HE(-ieta, iphi, theHBHETopology),
0802                                 false);
0803             spr::energyHCALCell(oppCell,
0804                                 hbhe,
0805                                 ehdepthoRaw,
0806                                 maxDepth_,
0807                                 -100.0,
0808                                 -100.0,
0809                                 -100.0,
0810                                 -100.0,
0811                                 -500.0,
0812                                 500.0,
0813                                 2,
0814                                 depth16HE(-ieta, iphi, theHBHETopology),
0815                                 false);
0816             for (unsigned int i = 0; i < ehdeptho.size(); ++i) {
0817               HcalSubdetector subdet0 =
0818                   (hborhe) ? ((ehdeptho[i].second >= depth16HE(-ieta, iphi, theHBHETopology)) ? HcalEndcap : HcalBarrel)
0819                            : subdet;
0820               HcalDetId hcid0(subdet0, -ieta, iphi, ehdeptho[i].second);
0821               double ene = ehdeptho[i].first;
0822               if (ene > 0.0) {
0823                 if (!(theHBHETopology->validHcal(hcid0))) {
0824                   edm::LogWarning("HBHEMuon") << "(3) Invalid ID " << hcid0 << " with E = " << ene;
0825                   edm::LogWarning("HBHEMuon") << oppCell << " with " << ehdeptho.size() << " depths:";
0826                   for (const auto& ehd : ehdeptho)
0827                     edm::LogWarning("HBHEMuon") << " " << ehd.second << ":" << ehd.first;
0828                 } else {
0829                   double chg(ene);
0830                   double corr = respCorr(DetId(hcid0), respCorrs);
0831                   if (corr != 0)
0832                     ene /= corr;
0833 #ifdef EDM_ML_DEBUG
0834                   HcalDetId id = (isItPlan1_ && isItPreRecHit_) ? hdc->mergedDepthDetId(hcid0) : hcid0;
0835                   edm::LogVerbatim("HBHEMuon")
0836                       << hcid0 << ":" << id << " Corr " << corr << " E " << ene << ":" << ehdeptho[i].first;
0837 #endif
0838                   double gain = gainFactor(conditions, hcid0);
0839                   if (gain != 0)
0840                     chg /= gain;
0841 #ifdef EDM_ML_DEBUG
0842                   edm::LogVerbatim("HBHEMuon") << hcid0 << " Gain " << gain << " C " << chg;
0843 #endif
0844                   int depth = ehdeptho[i].second - 1;
0845                   if (collapseDepth_) {
0846                     HcalDetId id = hdc->mergedDepthDetId(hcid0);
0847                     depth = id.depth() - 1;
0848                   }
0849                   cHcalDepthHotBG[depth] += chg;
0850 #ifdef EDM_ML_DEBUG
0851                   if ((verbosity_ % 10) > 0)
0852                     edm::LogVerbatim("HBHEMuon") << hcid0 << " Depth " << depth << " E " << ene << " C " << chg;
0853 #endif
0854                 }
0855               }
0856               double eneAux = ehdepthoAux[i].first;
0857               if (eneAux > 0.0) {
0858                 if (theHBHETopology->validHcal(hcid0)) {
0859                   double chgAux(eneAux);
0860                   double corr = respCorr(DetId(hcid0), respCorrs);
0861                   if (corr != 0)
0862                     eneAux /= corr;
0863                   double gain = gainFactor(conditions, hcid0);
0864                   if (gain != 0)
0865                     chgAux /= gain;
0866                   int depth = ehdepthoAux[i].second - 1;
0867                   if (collapseDepth_) {
0868                     HcalDetId id = hdc->mergedDepthDetId(hcid0);
0869                     depth = id.depth() - 1;
0870                   }
0871                   cHcalDepthHotBGAux[depth] += chgAux;
0872 #ifdef EDM_ML_DEBUG
0873                   if ((verbosity_ % 10) > 0)
0874                     edm::LogVerbatim("HBHEMuon") << hcid0 << " Depth " << depth << " E " << eneAux << " C " << chgAux;
0875 #endif
0876                 }
0877               }
0878               double eneRaw = ehdepthoRaw[i].first;
0879               if (eneRaw > 0.0) {
0880                 if (theHBHETopology->validHcal(hcid0)) {
0881                   double chgRaw(eneRaw);
0882                   double corr = respCorr(DetId(hcid0), respCorrs);
0883                   if (corr != 0)
0884                     eneRaw /= corr;
0885                   double gain = gainFactor(conditions, hcid0);
0886                   if (gain != 0)
0887                     chgRaw /= gain;
0888                   int depth = ehdepthoRaw[i].second - 1;
0889                   if (collapseDepth_) {
0890                     HcalDetId id = hdc->mergedDepthDetId(hcid0);
0891                     depth = id.depth() - 1;
0892                   }
0893                   cHcalDepthHotBGRaw[depth] += chgRaw;
0894 #ifdef EDM_ML_DEBUG
0895                   if ((verbosity_ % 10) > 0)
0896                     edm::LogVerbatim("HBHEMuon") << hcid0 << " Depth " << depth << " E " << eneRaw << " C " << chgRaw;
0897 #endif
0898                 }
0899               }
0900             }
0901           }
0902         }
0903 #ifdef EDM_ML_DEBUG
0904         edm::LogVerbatim("HBHEMuon") << "Propagate Track to HCAL: " << trackID.okHCAL << " Match "
0905                                      << hbheMuon.matchedId_ << " Hot " << isHot << " Energy " << eHcal;
0906 #endif
0907         hbheMuon.hcalIeta_ = ieta;
0908         hbheMuon.hcalIphi_ = iphi;
0909         for (int i = 0; i < depthMax_; ++i) {
0910           hbheMuon.hcalDepthEnergy_.push_back(eHcalDepth[i]);
0911           hbheMuon.hcalDepthActiveLength_.push_back(activeL[i]);
0912           hbheMuon.hcalDepthEnergyHot_.push_back(eHcalDepthHot[i]);
0913           hbheMuon.hcalDepthActiveLengthHot_.push_back(activeHotL[i]);
0914           hbheMuon.hcalDepthEnergyCorr_.push_back(eHcalDepthC[i]);
0915           hbheMuon.hcalDepthEnergyHotCorr_.push_back(eHcalDepthHotC[i]);
0916           hbheMuon.hcalDepthChargeHot_.push_back(cHcalDepthHot[i]);
0917           hbheMuon.hcalDepthChargeHotBG_.push_back(cHcalDepthHotBG[i]);
0918           hbheMuon.hcalDepthMatch_.push_back(matchDepth[i]);
0919           hbheMuon.hcalDepthMatchHot_.push_back(matchDepthHot[i]);
0920           hbheMuon.hcalDepthEnergyAux_.push_back(eHcalDepthAux[i]);
0921           hbheMuon.hcalDepthEnergyHotAux_.push_back(eHcalDepthHotAux[i]);
0922           hbheMuon.hcalDepthEnergyCorrAux_.push_back(eHcalDepthCAux[i]);
0923           hbheMuon.hcalDepthEnergyHotCorrAux_.push_back(eHcalDepthHotCAux[i]);
0924           hbheMuon.hcalDepthChargeHotAux_.push_back(cHcalDepthHotAux[i]);
0925           hbheMuon.hcalDepthChargeHotBGAux_.push_back(cHcalDepthHotBGAux[i]);
0926           hbheMuon.hcalDepthEnergyRaw_.push_back(eHcalDepthRaw[i]);
0927           hbheMuon.hcalDepthEnergyHotRaw_.push_back(eHcalDepthHotRaw[i]);
0928           hbheMuon.hcalDepthEnergyCorrRaw_.push_back(eHcalDepthCRaw[i]);
0929           hbheMuon.hcalDepthEnergyHotCorrRaw_.push_back(eHcalDepthHotCRaw[i]);
0930           hbheMuon.hcalDepthChargeHotRaw_.push_back(cHcalDepthHotRaw[i]);
0931           hbheMuon.hcalDepthChargeHotBGRaw_.push_back(cHcalDepthHotBGRaw[i]);
0932         }
0933         hbheMuon.hcalActiveLength_ = activeLengthTot;
0934         hbheMuon.hcalHot_ = isHot;
0935         hbheMuon.hcalActiveLengthHot_ = activeLengthHotTot;
0936 
0937         if ((recMuon->p() > 10.0) && (trackID.okHCAL))
0938           outputHcalHBHEMuonColl->emplace_back(hbheMuon);
0939       }
0940     }
0941   }
0942   if (!outputHcalHBHEMuonColl->empty())
0943     ++nGood_;
0944   iEvent.put(std::move(outputHcalHBHEMuonColl), labelHBHEMuon_);
0945 }
0946 
0947 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0948 void AlCaHcalHBHEMuonProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0949   edm::ParameterSetDescription desc;
0950   std::vector<std::string> triggers = {"HLT_IsoMu", "HLT_Mu"};
0951   desc.add<std::vector<std::string>>("triggers", triggers);
0952   desc.add<std::string>("processName", "HLT");
0953   desc.add<edm::InputTag>("triggerResults", edm::InputTag("TriggerResults", "", "HLT"));
0954   desc.add<edm::InputTag>("labelEBRecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0955   desc.add<edm::InputTag>("labelEERecHit", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0956   desc.add<edm::InputTag>("labelHBHERecHit", edm::InputTag("hbhereco"));
0957   desc.add<std::string>("labelVertex", "offlinePrimaryVertices");
0958   desc.add<std::string>("labelMuon", "muons");
0959   desc.add<std::string>("labelHBHEMuon", "hbheMuon");
0960   desc.add<bool>("collapseDepth", false);
0961   desc.add<bool>("isItPlan1", false);
0962   desc.addUntracked<int>("verbosity", 0);
0963   desc.addUntracked<bool>("isItPreRecHit", false);
0964   desc.addUntracked<bool>("writeRespCorr", false);
0965   desc.addUntracked<std::string>("fileInCorr", "");
0966   desc.addUntracked<int>("maxDepth", 4);
0967   descriptions.add("alcaHcalHBHEMuonProducer", desc);
0968 }
0969 
0970 // ------------ method called once each job just after ending the event loop  ------------
0971 void AlCaHcalHBHEMuonProducer::endStream() {
0972   globalCache()->nAll_ += nAll_;
0973   globalCache()->nGood_ += nGood_;
0974 }
0975 
0976 void AlCaHcalHBHEMuonProducer::globalEndJob(const alcaHcalHBHEMuon::Counters* count) {
0977   edm::LogVerbatim("HBHEMuon") << "Selects " << count->nGood_ << " out of " << count->nAll_ << " total # of events\n";
0978 }
0979 
0980 // ------------ method called when starting or ending a run  ------------
0981 void AlCaHcalHBHEMuonProducer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0982   const HcalDDDRecConstants* hdc = &iSetup.getData(tok_ddrec0_);
0983   actHB.clear();
0984   actHE.clear();
0985   actHB = hdc->getThickActive(0);
0986   actHE = hdc->getThickActive(1);
0987 #ifdef EDM_ML_DEBUG
0988   unsigned int k1(0), k2(0);
0989   edm::LogVerbatim("HBHEMuon") << actHB.size() << " Active Length for HB";
0990   for (const auto& act : actHB) {
0991     edm::LogVerbatim("HBHEMuon") << "[" << k1 << "] ieta " << act.ieta << " depth " << act.depth << " zside "
0992                                  << act.zside << " type " << act.stype << " phi " << act.iphis.size() << ":"
0993                                  << act.iphis[0] << " L " << act.thick;
0994     HcalDetId hcid1(HcalBarrel, (act.ieta) * (act.zside), act.iphis[0], act.depth);
0995     HcalDetId hcid2 = mergedDepth_ ? hdc->mergedDepthDetId(hcid1) : hcid1;
0996     edm::LogVerbatim("HBHEMuon") << hcid1 << " | " << hcid2 << " L " << activeLength(DetId(hcid2), hdc);
0997     ++k1;
0998   }
0999   edm::LogVerbatim("HBHEMuon") << actHE.size() << " Active Length for HE";
1000   for (const auto& act : actHE) {
1001     edm::LogVerbatim("HBHEMuon") << "[" << k2 << "] ieta " << act.ieta << " depth " << act.depth << " zside "
1002                                  << act.zside << " type " << act.stype << " phi " << act.iphis.size() << ":"
1003                                  << act.iphis[0] << " L " << act.thick;
1004     HcalDetId hcid1(HcalEndcap, (act.ieta) * (act.zside), act.iphis[0], act.depth);
1005     HcalDetId hcid2 = mergedDepth_ ? hdc->mergedDepthDetId(hcid1) : hcid1;
1006     edm::LogVerbatim("HBHEMuon") << hcid1 << " | " << hcid2 << " L " << activeLength(DetId(hcid2), hdc);
1007     ++k2;
1008   }
1009 #endif
1010 
1011   bool changed = true;
1012   bool flag = hltConfig_.init(iRun, iSetup, processName_, changed);
1013   edm::LogVerbatim("HBHEMuon") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init " << flag << std::endl;
1014 
1015   const HcalRespCorrs* resp = &iSetup.getData(tok_respcorr0_);
1016   const HcalTopology* theHBHETopology = &iSetup.getData(tok_htopo0_);
1017   const CaloGeometry* geo = &iSetup.getData(tok_geom0_);
1018   HcalRespCorrs respCorrsObj(*resp);
1019   HcalRespCorrs* respCorrs = &respCorrsObj;
1020   respCorrs->setTopo(theHBHETopology);
1021 
1022   // Write correction factors for all HB/HE events
1023   if (writeRespCorr_) {
1024     const HcalGeometry* gHcal = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
1025     const std::vector<DetId>& ids = gHcal->getValidDetIds(DetId::Hcal, 0);
1026     edm::LogVerbatim("HBHEMuon") << "\nTable of Correction Factors for Run " << iRun.run() << "\n";
1027     for (auto const& id : ids) {
1028       if ((id.det() == DetId::Hcal) && ((id.subdetId() == HcalBarrel) || (id.subdetId() == HcalEndcap))) {
1029         edm::LogVerbatim("HBHEMuon") << HcalDetId(id) << " " << id.rawId() << " "
1030                                      << (respCorrs->getValues(id))->getValue();
1031       }
1032     }
1033   }
1034 }
1035 
1036 void AlCaHcalHBHEMuonProducer::endRun(edm::Run const& iRun, edm::EventSetup const&) {
1037   edm::LogVerbatim("HBHEMuon") << "endRun[" << nRun_ << "] " << iRun.run() << "\n";
1038   ++nRun_;
1039 }
1040 
1041 // ------------ methods called by produce()  ------------
1042 
1043 int AlCaHcalHBHEMuonProducer::matchId(const HcalDetId& id1, const HcalDetId& id2) {
1044   HcalDetId kd1(id1.subdet(), id1.ieta(), id1.iphi(), 1);
1045   HcalDetId kd2(id1.subdet(), id2.ieta(), id2.iphi(), 1);
1046   int match = ((kd1 == kd2) ? 1 : 0);
1047   return match;
1048 }
1049 
1050 double AlCaHcalHBHEMuonProducer::activeLength(const DetId& id0, const HcalDDDRecConstants* hdc) {
1051   HcalDetId id(id0);
1052   int ieta = id.ietaAbs();
1053   int zside = id.zside();
1054   int iphi = id.iphi();
1055   std::vector<int> dpths;
1056   if (mergedDepth_) {
1057     std::vector<HcalDetId> ids;
1058     hdc->unmergeDepthDetId(id, ids);
1059     for (auto idh : ids)
1060       dpths.emplace_back(idh.depth());
1061   } else {
1062     dpths.emplace_back(id.depth());
1063   }
1064   double lx(0);
1065   if (id.subdet() == HcalBarrel) {
1066     for (unsigned int i = 0; i < actHB.size(); ++i) {
1067       if ((ieta == actHB[i].ieta) && (zside == actHB[i].zside) &&
1068           (std::find(dpths.begin(), dpths.end(), actHB[i].depth) != dpths.end()) &&
1069           (std::find(actHB[i].iphis.begin(), actHB[i].iphis.end(), iphi) != actHB[i].iphis.end())) {
1070         lx += actHB[i].thick;
1071       }
1072     }
1073   } else {
1074     for (unsigned int i = 0; i < actHE.size(); ++i) {
1075       if ((ieta == actHE[i].ieta) && (zside == actHE[i].zside) &&
1076           (std::find(dpths.begin(), dpths.end(), actHE[i].depth) != dpths.end()) &&
1077           (std::find(actHE[i].iphis.begin(), actHE[i].iphis.end(), iphi) != actHE[i].iphis.end())) {
1078         lx += actHE[i].thick;
1079       }
1080     }
1081   }
1082   return lx;
1083 }
1084 
1085 bool AlCaHcalHBHEMuonProducer::isGoodVertex(const reco::Vertex& vtx) {
1086   if (vtx.isFake())
1087     return false;
1088   if (vtx.ndof() < 4)
1089     return false;
1090   if (vtx.position().Rho() > 2.)
1091     return false;
1092   if (fabs(vtx.position().Z()) > 24.)
1093     return false;
1094   return true;
1095 }
1096 
1097 double AlCaHcalHBHEMuonProducer::respCorr(const DetId& id, const HcalRespCorrs* respCorrs) {
1098   double cfac(1.0);
1099   if (useMyCorr_) {
1100     auto itr = corrValue_.find(id);
1101     if (itr != corrValue_.end())
1102       cfac = itr->second;
1103   } else if (respCorrs != nullptr) {
1104     cfac = (respCorrs->getValues(id))->getValue();
1105   }
1106   return cfac;
1107 }
1108 
1109 double AlCaHcalHBHEMuonProducer::gainFactor(const HcalDbService* conditions, const HcalDetId& id) {
1110   double gain(0.0);
1111   const HcalCalibrations& calibs = conditions->getHcalCalibrations(id);
1112   for (int capid = 0; capid < 4; ++capid)
1113     gain += (0.25 * calibs.respcorrgain(capid));
1114   return gain;
1115 }
1116 
1117 int AlCaHcalHBHEMuonProducer::depth16HE(int ieta, int iphi, const HcalTopology* theHBHETopology) {
1118   // Transition between HB/HE is special
1119   // For Run 1 or for Plan1 standard reconstruction it is 3
1120   // For runs beyond 2018 or in Plan1 for HEP17 it is 4
1121   int zside = (ieta > 0) ? 1 : -1;
1122   int depth = theHBHETopology->dddConstants()->getMinDepth(1, 16, iphi, zside);
1123   if (isItPlan1_ && (!isItPreRecHit_))
1124     depth = 3;
1125 #ifdef EDM_ML_DEBUG
1126   edm::LogVerbatim("HBHEMuon") << "Plan1 " << isItPlan1_ << " PreRecHit " << isItPreRecHit_ << " phi " << iphi
1127                                << " depth " << depth;
1128 #endif
1129   return depth;
1130 }
1131 
1132 bool AlCaHcalHBHEMuonProducer::goodCell(const HcalDetId& hcid,
1133                                         const reco::Track* pTrack,
1134                                         const CaloGeometry* geo,
1135                                         const MagneticField* bField,
1136                                         const HcalDDDRecConstants* hdc) {
1137   std::pair<double, double> rz = hdc->getRZ(hcid);
1138   bool typeRZ = (hcid.subdet() == HcalEndcap) ? false : true;
1139   bool match = spr::propagateHCAL(pTrack, geo, bField, typeRZ, rz, (((verbosity_ / 10000) % 10) > 0));
1140   return match;
1141 }
1142 
1143 //define this as a plug-in
1144 #include "FWCore/Framework/interface/MakerMacros.h"
1145 
1146 DEFINE_FWK_MODULE(AlCaHcalHBHEMuonProducer);