diff --git a/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc b/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc
index bc32ef3278ee6..47c4de3418bc1 100644
--- a/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc
+++ b/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc
@@ -515,47 +515,50 @@ void HcalIsoTrkAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const
   t_L1Bit = true;
   t_TrigPass = false;
 
-  //L1
-  l1GtUtils_->retrieveL1(iEvent, iSetup, tok_alg_);
-  const std::vector<std::pair<std::string, bool> >& finalDecisions = l1GtUtils_->decisionsFinal();
-  for (const auto& decision : finalDecisions) {
-    if (decision.first.find(l1TrigName_) != std::string::npos) {
-      t_L1Bit = decision.second;
-      break;
+  edm::Handle<edm::TriggerResults> triggerResults;
+  if (!ignoreTrigger_) {
+    //L1
+    l1GtUtils_->retrieveL1(iEvent, iSetup, tok_alg_);
+    const std::vector<std::pair<std::string, bool> >& finalDecisions = l1GtUtils_->decisionsFinal();
+    for (const auto& decision : finalDecisions) {
+      if (decision.first.find(l1TrigName_) != std::string::npos) {
+        t_L1Bit = decision.second;
+        break;
+      }
     }
-  }
 #ifdef EDM_ML_DEBUG
-  edm::LogVerbatim("HcalIsoTrack") << "Trigger Information for " << l1TrigName_ << " is " << t_L1Bit
-                                   << " from a list of " << finalDecisions.size() << " decisions";
+    edm::LogVerbatim("HcalIsoTrack") << "Trigger Information for " << l1TrigName_ << " is " << t_L1Bit
+                                     << " from a list of " << finalDecisions.size() << " decisions";
 #endif
 
-  //HLT
-  edm::Handle<edm::TriggerResults> triggerResults;
-  iEvent.getByToken(tok_trigRes_, triggerResults);
-  if (triggerResults.isValid()) {
-    const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
-    const std::vector<std::string>& names = triggerNames.triggerNames();
-    if (!trigNames_.empty()) {
-      for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
-        int hlt = triggerResults->accept(iHLT);
-        for (unsigned int i = 0; i < trigNames_.size(); ++i) {
-          if (names[iHLT].find(trigNames_[i]) != std::string::npos) {
-            t_trgbits->at(i) = (hlt > 0);
-            t_hltbits->at(i) = (hlt > 0);
-            if (hlt > 0)
-              t_TrigPass = true;
+    //HLT
+    iEvent.getByToken(tok_trigRes_, triggerResults);
+    if (triggerResults.isValid()) {
+      const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
+      const std::vector<std::string>& names = triggerNames.triggerNames();
+      if (!trigNames_.empty()) {
+        for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
+          int hlt = triggerResults->accept(iHLT);
+          for (unsigned int i = 0; i < trigNames_.size(); ++i) {
+            if (names[iHLT].find(trigNames_[i]) != std::string::npos) {
+              t_trgbits->at(i) = (hlt > 0);
+              t_hltbits->at(i) = (hlt > 0);
+              if (hlt > 0)
+                t_TrigPass = true;
 #ifdef EDM_ML_DEBUG
-            edm::LogVerbatim("HcalIsoTrack")
-                << "This trigger " << names[iHLT] << " Flag " << hlt << ":" << t_trgbits->at(i);
+              edm::LogVerbatim("HcalIsoTrack")
+                  << "This trigger " << names[iHLT] << " Flag " << hlt << ":" << t_trgbits->at(i);
 #endif
+            }
           }
         }
       }
     }
-  }
 #ifdef EDM_ML_DEBUG
-  edm::LogVerbatim("HcalIsoTrack") << "HLT Information shows " << t_TrigPass << ":" << trigNames_.empty() << ":" << okC;
+    edm::LogVerbatim("HcalIsoTrack") << "HLT Information shows " << t_TrigPass << ":" << trigNames_.empty() << ":"
+                                     << okC;
 #endif
+  }
 
   std::array<int, 3> ntksave{{0, 0, 0}};
   if (ignoreTrigger_ || useL1Trigger_) {
diff --git a/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc b/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc
new file mode 100644
index 0000000000000..684111e5f2484
--- /dev/null
+++ b/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc
@@ -0,0 +1,457 @@
+#include <memory>
+#include <iostream>
+#include <vector>
+
+#include <TTree.h>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "FWCore/Common/interface/TriggerNames.h"
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+
+#include "DataFormats/EcalDetId/interface/EBDetId.h"
+#include "DataFormats/EcalDetId/interface/EEDetId.h"
+#include "DataFormats/HcalDetId/interface/HcalDetId.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
+#include "SimDataFormats/Track/interface/SimTrack.h"
+#include "SimDataFormats/Track/interface/SimTrackContainer.h"
+#include "SimDataFormats/Vertex/interface/SimVertex.h"
+#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
+
+#include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
+#include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
+#include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
+
+#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
+#include "Geometry/Records/interface/CaloGeometryRecord.h"
+#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
+#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
+#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
+#include "Geometry/CaloTopology/interface/HcalTopology.h"
+#include "Geometry/CaloTopology/interface/CaloTopology.h"
+#include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
+#include "MagneticField/Engine/interface/MagneticField.h"
+#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
+
+//#define EDM_ML_DEBUG
+
+class HcalHBHEMuonSimAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
+public:
+  explicit HcalHBHEMuonSimAnalyzer(const edm::ParameterSet&);
+  ~HcalHBHEMuonSimAnalyzer();
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  virtual void beginJob() override;
+  virtual void analyze(edm::Event const&, edm::EventSetup const&) override;
+  virtual void beginRun(edm::Run const&, edm::EventSetup const&) override;
+  virtual void endRun(edm::Run const&, edm::EventSetup const&) override {}
+  virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
+  virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
+  void clearVectors();
+  unsigned int matchId(const HcalDetId&, const HcalDetId&);
+  double activeLength(const DetId&);
+
+  std::string g4Label_, ebLabel_, eeLabel_;
+  std::string hcLabel_;
+  int verbosity_, maxDepth_;
+  double etaMax_;
+  std::vector<HcalDDDRecConstants::HcalActiveLength> actHB_, actHE_;
+
+  double tMinE_, tMaxE_, tMinH_, tMaxH_;
+  edm::Service<TFileService> fs_;
+  edm::EDGetTokenT<edm::SimTrackContainer> tok_SimTk_;
+  edm::EDGetTokenT<edm::SimVertexContainer> tok_SimVtx_;
+  edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloEB_, tok_caloEE_;
+  edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloHH_;
+  edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_HRNDC_;
+  edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_hdcons_;
+  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
+  edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
+  edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_topo_;
+  edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
+
+  static const int depthMax_ = 7;
+  const int idMuon_ = 13;
+  TTree* tree_;
+  unsigned int runNumber_, eventNumber_, lumiNumber_, bxNumber_;
+  std::vector<double> ptGlob_, etaGlob_, phiGlob_, pMuon_;
+  std::vector<double> ecal3x3Energy_, hcal1x1Energy_;
+  std::vector<unsigned int> ecalDetId_, hcalDetId_, hcalHot_;
+  std::vector<double> matchedId_;
+  std::vector<double> hcalDepthEnergy_[depthMax_];
+  std::vector<double> hcalDepthActiveLength_[depthMax_];
+  std::vector<double> hcalDepthEnergyHot_[depthMax_];
+  std::vector<double> hcalDepthActiveLengthHot_[depthMax_];
+  std::vector<double> hcalActiveLength_, hcalActiveLengthHot_;
+};
+
+HcalHBHEMuonSimAnalyzer::HcalHBHEMuonSimAnalyzer(const edm::ParameterSet& iConfig) {
+  usesResource(TFileService::kSharedResource);
+
+  //now do what ever initialization is needed
+  g4Label_ = iConfig.getParameter<std::string>("ModuleLabel");
+  ebLabel_ = iConfig.getParameter<std::string>("EBCollection");
+  eeLabel_ = iConfig.getParameter<std::string>("EECollection");
+  hcLabel_ = iConfig.getParameter<std::string>("HCCollection");
+  verbosity_ = iConfig.getUntrackedParameter<int>("Verbosity", 0);
+  maxDepth_ = iConfig.getUntrackedParameter<int>("MaxDepth", 4);
+  etaMax_ = iConfig.getUntrackedParameter<double>("EtaMax", 3.0);
+  tMinE_ = iConfig.getUntrackedParameter<double>("TimeMinCutECAL", -500.);
+  tMaxE_ = iConfig.getUntrackedParameter<double>("TimeMaxCutECAL", 500.);
+  tMinH_ = iConfig.getUntrackedParameter<double>("TimeMinCutHCAL", -500.);
+  tMaxH_ = iConfig.getUntrackedParameter<double>("TimeMaxCutHCAL", 500.);
+
+  tok_SimTk_ = consumes<edm::SimTrackContainer>(edm::InputTag(g4Label_));
+  tok_SimVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag(g4Label_));
+  tok_caloEB_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, ebLabel_));
+  tok_caloEE_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, eeLabel_));
+  tok_caloHH_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hcLabel_));
+  if (maxDepth_ > depthMax_)
+    maxDepth_ = depthMax_;
+  else if (maxDepth_ < 1)
+    maxDepth_ = 4;
+#ifdef EDM_ML_DEBUG
+  std::cout << "Labels: " << g4Label_ << ":" << ebLabel_ << ":" << eeLabel_ << ":" << hcLabel_ << "\nVerbosity "
+            << verbosity_ << " MaxDepth " << maxDepth_ << " Maximum Eta " << etaMax_ << " tMin|tMax " << tMinE_ << ":"
+            << tMaxE_ << ":" << tMinH_ << ":" << tMaxH_ << std::endl;
+#endif
+  tok_HRNDC_ = esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>();
+  tok_hdcons_ = esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>();
+  tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
+  tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
+  tok_topo_ = esConsumes<CaloTopology, CaloTopologyRecord>();
+  tok_htopo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
+}
+
+HcalHBHEMuonSimAnalyzer::~HcalHBHEMuonSimAnalyzer() {}
+
+void HcalHBHEMuonSimAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  clearVectors();
+  bool debug(false);
+#ifdef EDM_ML_DEBUG
+  debug = ((verbosity_ / 10) > 0);
+#endif
+  // depthHE is the first depth index for HE for |ieta| = 16
+  // It used to be 3 for all runs preceding 2017 and 4 beyond that
+  int depthHE = (maxDepth_ <= 6) ? 3 : 4;
+
+  const HcalDDDRecConstants* hcons = &iSetup.getData(tok_hdcons_);
+
+  runNumber_ = iEvent.id().run();
+  eventNumber_ = iEvent.id().event();
+  lumiNumber_ = iEvent.id().luminosityBlock();
+  bxNumber_ = iEvent.bunchCrossing();
+
+  //get Handles to SimTracks and SimHits
+  edm::Handle<edm::SimTrackContainer> SimTk;
+  iEvent.getByToken(tok_SimTk_, SimTk);
+  edm::Handle<edm::SimVertexContainer> SimVtx;
+  iEvent.getByToken(tok_SimVtx_, SimVtx);
+
+  //get Handles to PCaloHitContainers of eb/ee/hbhe
+  edm::Handle<edm::PCaloHitContainer> pcaloeb;
+  iEvent.getByToken(tok_caloEB_, pcaloeb);
+  edm::Handle<edm::PCaloHitContainer> pcaloee;
+  iEvent.getByToken(tok_caloEE_, pcaloee);
+  edm::Handle<edm::PCaloHitContainer> pcalohh;
+  iEvent.getByToken(tok_caloHH_, pcalohh);
+  std::vector<PCaloHit> calohh;
+  bool testN(false);
+  for (unsigned int k = 1; k < pcalohh->size(); ++k) {
+    // if it is a standard DetId bits 28..31 will carry the det #
+    // for HCAL det # is 4 and if there is at least one hit in the collection
+    // have det # which is not 4 this collection is created using TestNumbering
+    int det = ((((*pcalohh)[k].id()) >> 28) & 0xF);
+    if (det != 4) {
+      testN = true;
+      break;
+    }
+  }
+  if (testN) {
+    for (edm::PCaloHitContainer::const_iterator itr = pcalohh->begin(); itr != pcalohh->end(); ++itr) {
+      PCaloHit hit(*itr);
+      DetId newid = HcalHitRelabeller::relabel(hit.id(), hcons);
+#ifdef EDM_ML_DEBUG
+      std::cout << "Old ID " << std::hex << hit.id() << std::dec << " New " << HcalDetId(newid) << std::endl;
+#endif
+      hit.setID(newid.rawId());
+      calohh.push_back(hit);
+    }
+  } else {
+    calohh.insert(calohh.end(), pcalohh->begin(), pcalohh->end());
+  }
+
+  // get handles to calogeometry and calotopology
+  const CaloGeometry* geo = &iSetup.getData(tok_geom_);
+  const MagneticField* bField = &iSetup.getData(tok_magField_);
+  const CaloTopology* caloTopology = &iSetup.getData(tok_topo_);
+  const HcalTopology* theHBHETopology = &iSetup.getData(tok_htopo_);
+
+  // Loop over all SimTracks
+  for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
+    if ((std::abs(simTrkItr->type()) == idMuon_) && (simTrkItr->vertIndex() == 0) &&
+        (std::abs(simTrkItr->momentum().eta()) < etaMax_)) {
+      unsigned int thisTrk = simTrkItr->trackId();
+      spr::propagatedTrackDirection trkD = spr::propagateCALO(thisTrk, SimTk, SimVtx, geo, bField, debug);
+
+      double eEcal(0), eHcal(0), activeLengthTot(0), activeLengthHotTot(0);
+      double eHcalDepth[depthMax_], eHcalDepthHot[depthMax_];
+      double activeL[depthMax_], activeHotL[depthMax_];
+      unsigned int isHot(0);
+      bool tmpmatch(false);
+      for (int i = 0; i < depthMax_; ++i)
+        eHcalDepth[i] = eHcalDepthHot[i] = activeL[i] = activeHotL[i] = -10000;
+
+#ifdef EDM_ML_DEBUG
+      if ((verbosity_ % 10) > 0)
+        std::cout << "Track Type " << simTrkItr->type() << " Vertex " << simTrkItr->vertIndex() << " Charge "
+                  << simTrkItr->charge() << " Momentum " << simTrkItr->momentum().P() << ":"
+                  << simTrkItr->momentum().eta() << ":" << simTrkItr->momentum().phi() << " ECAL|HCAL " << trkD.okECAL
+                  << ":" << trkD.okHCAL << " Point " << trkD.pointECAL << ":" << trkD.pointHCAL << " Direction "
+                  << trkD.directionECAL.eta() << ":" << trkD.directionECAL.phi() << " | " << trkD.directionHCAL.eta()
+                  << ":" << trkD.directionHCAL.phi() << std::endl;
+#endif
+      bool propageback(false);
+      spr::propagatedTrackDirection trkD_back = spr::propagateHCALBack(thisTrk, SimTk, SimVtx, geo, bField, debug);
+      HcalDetId closestCell_back;
+      if (trkD_back.okHCAL) {
+        closestCell_back = (HcalDetId)(trkD_back.detIdHCAL);
+        propageback = true;
+      }
+      if (trkD.okHCAL) {
+        // Muon properties
+        spr::trackAtOrigin tkvx = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
+        ptGlob_.push_back(tkvx.momentum.perp());
+        etaGlob_.push_back(tkvx.momentum.eta());
+        phiGlob_.push_back(tkvx.momentum.phi());
+        pMuon_.push_back(tkvx.momentum.mag());
+#ifdef EDM_ML_DEBUG
+        if ((verbosity_ % 10) > 0)
+          std::cout << "Track at vertex " << tkvx.ok << " position " << tkvx.position << " Momentum "
+                    << tkvx.momentum.mag() << ":" << tkvx.momentum.eta() << ":" << tkvx.momentum.phi() << " Charge "
+                    << tkvx.charge << std::endl;
+#endif
+
+        // Energy in ECAL
+        DetId isoCell;
+        if (trkD.okECAL) {
+          isoCell = trkD.detIdECAL;
+          eEcal = spr::eECALmatrix(
+              isoCell, pcaloeb, pcaloee, geo, caloTopology, 1, 1, -100.0, -100.0, tMinE_, tMaxE_, debug);
+        }
+
+        // Energy in  Hcal
+        const DetId closestCell(trkD.detIdHCAL);
+        if ((propageback) && (HcalDetId(closestCell).ieta() == HcalDetId(closestCell_back).ieta()) &&
+            (HcalDetId(closestCell).iphi() == HcalDetId(closestCell_back).iphi()))
+          tmpmatch = true;
+
+        eHcal = spr::eHCALmatrix(
+            theHBHETopology, closestCell, calohh, 0, 0, false, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_, debug);
+#ifdef EDM_ML_DEBUG
+        if ((verbosity_ % 10) > 0)
+          std::cout << "eEcal " << trkD.okECAL << ":" << eEcal << " eHcal " << eHcal << std::endl;
+#endif
+
+        HcalSubdetector subdet = HcalDetId(closestCell).subdet();
+        int ieta = HcalDetId(closestCell).ieta();
+        int iphi = HcalDetId(closestCell).iphi();
+        bool hbhe = (std::abs(ieta) == 16);
+        std::vector<std::pair<double, int> > ehdepth;
+        spr::energyHCALCell(
+            (HcalDetId)closestCell, calohh, ehdepth, maxDepth_, -100.0, -100.0, -100.0, -100.0, -500.0, 500.0, debug);
+        for (unsigned int i = 0; i < ehdepth.size(); ++i) {
+          eHcalDepth[ehdepth[i].second - 1] = ehdepth[i].first;
+          HcalSubdetector subdet0 = (hbhe) ? ((ehdepth[i].second >= depthHE) ? HcalEndcap : HcalBarrel) : subdet;
+          HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
+          double actL = activeLength(DetId(hcid0));
+          activeL[ehdepth[i].second - 1] = actL;
+          activeLengthTot += actL;
+#ifdef EDM_ML_DEBUG
+          if ((verbosity_ % 10) > 0)
+            std::cout << hcid0 << " E " << ehdepth[i].first << " L " << actL << std::endl;
+#endif
+        }
+
+        HcalDetId hotCell;
+#ifdef EDM_ML_DEBUG
+        double h3x3 =
+#endif
+            spr::eHCALmatrix(geo, theHBHETopology, closestCell, calohh, 1, 1, hotCell, debug);
+        isHot = matchId(closestCell, hotCell);
+#ifdef EDM_ML_DEBUG
+        if ((verbosity_ % 10) > 0)
+          std::cout << "hcal 3X3  < " << h3x3 << ">"
+                    << " ClosestCell <" << (HcalDetId)(closestCell) << "> hotCell id < " << hotCell << "> isHot"
+                    << isHot << std::endl;
+#endif
+
+        if (hotCell != HcalDetId()) {
+          subdet = HcalDetId(hotCell).subdet();
+          ieta = HcalDetId(hotCell).ieta();
+          iphi = HcalDetId(hotCell).iphi();
+          hbhe = (std::abs(ieta) == 16);
+          std::vector<std::pair<double, int> > ehdepth;
+          spr::energyHCALCell(
+              hotCell, calohh, ehdepth, maxDepth_, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_, debug);
+          for (unsigned int i = 0; i < ehdepth.size(); ++i) {
+            eHcalDepthHot[ehdepth[i].second - 1] = ehdepth[i].first;
+            HcalSubdetector subdet0 = (hbhe) ? ((ehdepth[i].second >= depthHE) ? HcalEndcap : HcalBarrel) : subdet;
+            HcalDetId hcid0(subdet0, ieta, iphi, ehdepth[i].second);
+            double actL = activeLength(DetId(hcid0));
+            activeHotL[ehdepth[i].second - 1] = actL;
+            activeLengthHotTot += actL;
+#ifdef EDM_ML_DEBUG
+            if ((verbosity_ % 10) > 0)
+              std::cout << hcid0 << " E " << ehdepth[i].first << " L " << actL << std::endl;
+#endif
+          }
+        }
+#ifdef EDM_ML_DEBUG
+        if ((verbosity_ % 10) > 0) {
+          for (int k = 0; k < depthMax_; ++k)
+            std::cout << "Depth " << k << " E " << eHcalDepth[k] << ":" << eHcalDepthHot[k] << std::endl;
+        }
+#endif
+        matchedId_.push_back(tmpmatch);
+        ecal3x3Energy_.push_back(eEcal);
+        ecalDetId_.push_back(isoCell.rawId());
+        hcal1x1Energy_.push_back(eHcal);
+        hcalDetId_.push_back(closestCell.rawId());
+        for (int k = 0; k < depthMax_; ++k) {
+          hcalDepthEnergy_[k].push_back(eHcalDepth[k]);
+          hcalDepthActiveLength_[k].push_back(activeL[k]);
+          hcalDepthEnergyHot_[k].push_back(eHcalDepthHot[k]);
+          hcalDepthActiveLengthHot_[k].push_back(activeHotL[k]);
+        }
+        hcalHot_.push_back(isHot);
+        hcalActiveLengthHot_.push_back(activeLengthHotTot);
+      }
+    }
+  }
+  if (hcalHot_.size() > 0)
+    tree_->Fill();
+}
+
+void HcalHBHEMuonSimAnalyzer::beginJob() {
+  tree_ = fs_->make<TTree>("TREE", "TREE");
+  tree_->Branch("Run_No", &runNumber_);
+  tree_->Branch("Event_No", &eventNumber_);
+  tree_->Branch("LumiNumber", &lumiNumber_);
+  tree_->Branch("BXNumber", &bxNumber_);
+  tree_->Branch("pt_of_muon", &ptGlob_);
+  tree_->Branch("eta_of_muon", &etaGlob_);
+  tree_->Branch("phi_of_muon", &phiGlob_);
+  tree_->Branch("p_of_muon", &pMuon_);
+  tree_->Branch("matchedId", &matchedId_);
+
+  tree_->Branch("ecal_3x3", &ecal3x3Energy_);
+  tree_->Branch("ecal_detID", &ecalDetId_);
+  tree_->Branch("hcal_1x1", &hcal1x1Energy_);
+  tree_->Branch("hcal_detID", &hcalDetId_);
+  tree_->Branch("hcal_cellHot", &hcalHot_);
+  tree_->Branch("activeLength", &hcalActiveLength_);
+  tree_->Branch("activeLengthHot", &hcalActiveLengthHot_);
+  char name[100];
+  for (int k = 0; k < maxDepth_; ++k) {
+    sprintf(name, "hcal_edepth%d", (k + 1));
+    tree_->Branch(name, &hcalDepthEnergy_[k]);
+    sprintf(name, "hcal_activeL%d", (k + 1));
+    tree_->Branch(name, &hcalDepthActiveLength_[k]);
+    sprintf(name, "hcal_edepthHot%d", (k + 1));
+    tree_->Branch(name, &hcalDepthEnergyHot_[k]);
+    sprintf(name, "hcal_activeHotL%d", (k + 1));
+    tree_->Branch(name, &hcalDepthActiveLength_[k]);
+  }
+}
+
+void HcalHBHEMuonSimAnalyzer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
+  const HcalDDDRecConstants& hdc = iSetup.getData(tok_HRNDC_);
+  actHB_.clear();
+  actHE_.clear();
+  actHB_ = hdc.getThickActive(0);
+  actHE_ = hdc.getThickActive(1);
+}
+
+void HcalHBHEMuonSimAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<std::string>("ModuleLabel", "g4SimHits");
+  desc.add<std::string>("EBCollection", "EcalHitsEB");
+  desc.add<std::string>("EECollection", "EcalHitsEE");
+  desc.add<std::string>("HCCollection", "HcalHits");
+  desc.addUntracked<int>("Verbosity", 0);
+  desc.addUntracked<int>("MaxDepth", 4);
+  desc.addUntracked<double>("EtaMax", 3.0);
+  desc.addUntracked<double>("TimeMinCutECAL", -500.);
+  desc.addUntracked<double>("TimeMaxCutECAL", 500.);
+  desc.addUntracked<double>("TimeMinCutHCAL", -500.);
+  desc.addUntracked<double>("TimeMaxCutHCAL", 500.);
+  descriptions.add("hcalHBHEMuonSim", desc);
+}
+
+void HcalHBHEMuonSimAnalyzer::clearVectors() {
+  ///clearing vectots
+  runNumber_ = -99999;
+  eventNumber_ = -99999;
+  lumiNumber_ = -99999;
+  bxNumber_ = -99999;
+
+  ptGlob_.clear();
+  etaGlob_.clear();
+  phiGlob_.clear();
+  pMuon_.clear();
+  matchedId_.clear();
+  ecal3x3Energy_.clear();
+  ecalDetId_.clear();
+  hcal1x1Energy_.clear();
+  hcalDetId_.clear();
+  hcalHot_.clear();
+  hcalActiveLength_.clear();
+  hcalActiveLengthHot_.clear();
+  for (int k = 0; k < depthMax_; ++k) {
+    hcalDepthEnergy_[k].clear();
+    hcalDepthActiveLength_[k].clear();
+    hcalDepthEnergyHot_[k].clear();
+    hcalDepthActiveLengthHot_[k].clear();
+  }
+}
+
+unsigned int HcalHBHEMuonSimAnalyzer::matchId(const HcalDetId& id1, const HcalDetId& id2) {
+  HcalDetId kd1(id1.subdet(), id1.ieta(), id1.iphi(), 1);
+  HcalDetId kd2(id2.subdet(), id2.ieta(), id2.iphi(), 1);
+  unsigned int match = ((kd1 == kd2) ? 1 : 0);
+  return match;
+}
+
+double HcalHBHEMuonSimAnalyzer::activeLength(const DetId& id_) {
+  HcalDetId id(id_);
+  int ieta = id.ietaAbs();
+  int depth = id.depth();
+  double lx(0);
+  if (id.subdet() == HcalBarrel) {
+    for (unsigned int i = 0; i < actHB_.size(); ++i) {
+      if (ieta == actHB_[i].ieta && depth == actHB_[i].depth) {
+        lx = actHB_[i].thick;
+        break;
+      }
+    }
+  } else {
+    for (unsigned int i = 0; i < actHE_.size(); ++i) {
+      if (ieta == actHE_[i].ieta && depth == actHE_[i].depth) {
+        lx = actHE_[i].thick;
+        break;
+      }
+    }
+  }
+  return lx;
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(HcalHBHEMuonSimAnalyzer);