From 2f38e0e93a47296a0c921d57f15047999ac1910a Mon Sep 17 00:00:00 2001 From: Sunanda Date: Fri, 21 May 2021 16:05:46 +0200 Subject: [PATCH 1/3] Add a new analyzer for studies of HCAL using isotrack and muons --- .../plugins/HcalIsoTrkAnalyzer.cc | 62 +-- .../test/HcalHBHEMuonSimAnalyzer.cc | 476 ++++++++++++++++++ 2 files changed, 507 insertions(+), 31 deletions(-) create mode 100644 Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc diff --git a/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc b/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc index bc32ef3278ee6..1f9feda8a8eba 100644 --- a/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc +++ b/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc @@ -515,47 +515,47 @@ 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 >& finalDecisions = l1GtUtils_->decisionsFinal(); - for (const auto& decision : finalDecisions) { - if (decision.first.find(l1TrigName_) != std::string::npos) { - t_L1Bit = decision.second; - break; + edm::Handle triggerResults; + if (!ignoreTrigger_) { + //L1 + l1GtUtils_->retrieveL1(iEvent, iSetup, tok_alg_); + const std::vector >& 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 triggerResults; - iEvent.getByToken(tok_trigRes_, triggerResults); - if (triggerResults.isValid()) { - const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults); - const std::vector& 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& 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 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..77576b69e5a09 --- /dev/null +++ b/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc @@ -0,0 +1,476 @@ +#include +#include +#include + +#include + +// 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 { + +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 actHB_, actHE_; + + double tMinE_, tMaxE_, tMinH_, tMaxH_; + edm::Service fs_; + edm::EDGetTokenT tok_SimTk_; + edm::EDGetTokenT tok_SimVtx_; + edm::EDGetTokenT tok_caloEB_, tok_caloEE_; + edm::EDGetTokenT tok_caloHH_; + edm::ESGetToken tok_HRNDC_; + edm::ESGetToken tok_hdcons_; + edm::ESGetToken tok_geom_; + edm::ESGetToken tok_magField_; + edm::ESGetToken tok_topo_; + edm::ESGetToken tok_htopo_; + + static const int depthMax_=7; + const int idMuon_=13; + TTree *tree_; + unsigned int runNumber_, eventNumber_, lumiNumber_, bxNumber_; + std::vector ptGlob_, etaGlob_, phiGlob_, pMuon_; + std::vector ecal3x3Energy_, hcal1x1Energy_; + std::vector ecalDetId_, hcalDetId_, hcalHot_; + std::vector matchedId_; + std::vector hcalDepthEnergy_[depthMax_]; + std::vector hcalDepthActiveLength_[depthMax_]; + std::vector hcalDepthEnergyHot_[depthMax_]; + std::vector hcalDepthActiveLengthHot_[depthMax_]; + std::vector hcalActiveLength_, hcalActiveLengthHot_; +}; + +HcalHBHEMuonSimAnalyzer::HcalHBHEMuonSimAnalyzer(const edm::ParameterSet& iConfig) { + + usesResource(TFileService::kSharedResource); + + //now do what ever initialization is needed + g4Label_ = iConfig.getParameter("ModuleLabel"); + ebLabel_ = iConfig.getParameter("EBCollection"); + eeLabel_ = iConfig.getParameter("EECollection"); + hcLabel_ = iConfig.getParameter("HCCollection"); + verbosity_ = iConfig.getUntrackedParameter("Verbosity",0); + maxDepth_ = iConfig.getUntrackedParameter("MaxDepth",4); + etaMax_ = iConfig.getUntrackedParameter("EtaMax", 3.0); + tMinE_ = iConfig.getUntrackedParameter("TimeMinCutECAL", -500.); + tMaxE_ = iConfig.getUntrackedParameter("TimeMaxCutECAL", 500.); + tMinH_ = iConfig.getUntrackedParameter("TimeMinCutHCAL", -500.); + tMaxH_ = iConfig.getUntrackedParameter("TimeMaxCutHCAL", 500.); + + tok_SimTk_ = consumes(edm::InputTag(g4Label_)); + tok_SimVtx_ = consumes(edm::InputTag(g4Label_)); + tok_caloEB_ = consumes(edm::InputTag(g4Label_,ebLabel_)); + tok_caloEE_ = consumes(edm::InputTag(g4Label_,eeLabel_)); + tok_caloHH_ = consumes(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(); + tok_hdcons_ = esConsumes(); + tok_geom_ = esConsumes(); + tok_magField_ = esConsumes(); + tok_topo_ = esConsumes(); + tok_htopo_ = esConsumes(); +} + +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 SimTk; + iEvent.getByToken(tok_SimTk_,SimTk); + edm::SimTrackContainer::const_iterator simTrkItr; + edm::Handle SimVtx; + iEvent.getByToken(tok_SimVtx_,SimVtx); + + //get Handles to PCaloHitContainers of eb/ee/hbhe + edm::Handle pcaloeb; + iEvent.getByToken(tok_caloEB_, pcaloeb); + edm::Handle pcaloee; + iEvent.getByToken(tok_caloEE_, pcaloee); + edm::Handle pcalohh; + iEvent.getByToken(tok_caloHH_, pcalohh); + std::vector calohh; + bool testN(false); + for (unsigned int k=1; ksize(); ++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 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 > 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= 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 > ehdepth; + spr::energyHCALCell(hotCell, calohh, ehdepth, maxDepth_, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_, debug); + for (unsigned int i=0; i= 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 0) tree_->Fill(); +} + +void HcalHBHEMuonSimAnalyzer::beginJob() { + + tree_ = fs_->make("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; kBranch(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("ModuleLabel","g4SimHits"); + desc.add("EBCollection","EcalHitsEB"); + desc.add("EECollection","EcalHitsEE"); + desc.add("HCCollection","HcalHits"); + desc.addUntracked("Verbosity",0); + desc.addUntracked("MaxDepth",4); + desc.addUntracked("EtaMax",3.0); + desc.addUntracked("TimeMinCutECAL",-500.); + desc.addUntracked("TimeMaxCutECAL",500.); + desc.addUntracked("TimeMinCutHCAL",-500.); + desc.addUntracked("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 Date: Fri, 21 May 2021 16:17:54 +0200 Subject: [PATCH 2/3] Code check --- .../plugins/HcalIsoTrkAnalyzer.cc | 35 +- .../test/HcalHBHEMuonSimAnalyzer.cc | 498 +++++++++--------- 2 files changed, 259 insertions(+), 274 deletions(-) diff --git a/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc b/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc index 1f9feda8a8eba..47c4de3418bc1 100644 --- a/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc +++ b/Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc @@ -522,12 +522,13 @@ void HcalIsoTrkAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const const std::vector >& finalDecisions = l1GtUtils_->decisionsFinal(); for (const auto& decision : finalDecisions) { if (decision.first.find(l1TrigName_) != std::string::npos) { - t_L1Bit = decision.second; - break; + 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 @@ -536,24 +537,26 @@ void HcalIsoTrkAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults); const std::vector& 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; + 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 } diff --git a/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc b/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc index 77576b69e5a09..f87e908dbc204 100644 --- a/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc +++ b/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc @@ -24,8 +24,8 @@ #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 "Calibration/IsolatedParticles/interface/eECALMatrix.h" +#include "Calibration/IsolatedParticles/interface/eHCALMatrix.h" #include "Geometry/CaloGeometry/interface/CaloGeometry.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" @@ -40,8 +40,7 @@ #define EDM_ML_DEBUG -class HcalHBHEMuonSimAnalyzer : public edm::one::EDAnalyzer { - +class HcalHBHEMuonSimAnalyzer : public edm::one::EDAnalyzer { public: explicit HcalHBHEMuonSimAnalyzer(const edm::ParameterSet&); ~HcalHBHEMuonSimAnalyzer(); @@ -55,74 +54,73 @@ class HcalHBHEMuonSimAnalyzer : public edm::one::EDAnalyzer actHB_, actHE_; - double tMinE_, tMaxE_, tMinH_, tMaxH_; - edm::Service fs_; - edm::EDGetTokenT tok_SimTk_; - edm::EDGetTokenT tok_SimVtx_; - edm::EDGetTokenT tok_caloEB_, tok_caloEE_; - edm::EDGetTokenT tok_caloHH_; + double tMinE_, tMaxE_, tMinH_, tMaxH_; + edm::Service fs_; + edm::EDGetTokenT tok_SimTk_; + edm::EDGetTokenT tok_SimVtx_; + edm::EDGetTokenT tok_caloEB_, tok_caloEE_; + edm::EDGetTokenT tok_caloHH_; edm::ESGetToken tok_HRNDC_; edm::ESGetToken tok_hdcons_; - edm::ESGetToken tok_geom_; - edm::ESGetToken tok_magField_; - edm::ESGetToken tok_topo_; - edm::ESGetToken tok_htopo_; - - static const int depthMax_=7; - const int idMuon_=13; - TTree *tree_; - unsigned int runNumber_, eventNumber_, lumiNumber_, bxNumber_; - std::vector ptGlob_, etaGlob_, phiGlob_, pMuon_; - std::vector ecal3x3Energy_, hcal1x1Energy_; - std::vector ecalDetId_, hcalDetId_, hcalHot_; - std::vector matchedId_; - std::vector hcalDepthEnergy_[depthMax_]; - std::vector hcalDepthActiveLength_[depthMax_]; - std::vector hcalDepthEnergyHot_[depthMax_]; - std::vector hcalDepthActiveLengthHot_[depthMax_]; - std::vector hcalActiveLength_, hcalActiveLengthHot_; + edm::ESGetToken tok_geom_; + edm::ESGetToken tok_magField_; + edm::ESGetToken tok_topo_; + edm::ESGetToken tok_htopo_; + + static const int depthMax_ = 7; + const int idMuon_ = 13; + TTree* tree_; + unsigned int runNumber_, eventNumber_, lumiNumber_, bxNumber_; + std::vector ptGlob_, etaGlob_, phiGlob_, pMuon_; + std::vector ecal3x3Energy_, hcal1x1Energy_; + std::vector ecalDetId_, hcalDetId_, hcalHot_; + std::vector matchedId_; + std::vector hcalDepthEnergy_[depthMax_]; + std::vector hcalDepthActiveLength_[depthMax_]; + std::vector hcalDepthEnergyHot_[depthMax_]; + std::vector hcalDepthActiveLengthHot_[depthMax_]; + std::vector hcalActiveLength_, hcalActiveLengthHot_; }; HcalHBHEMuonSimAnalyzer::HcalHBHEMuonSimAnalyzer(const edm::ParameterSet& iConfig) { - usesResource(TFileService::kSharedResource); //now do what ever initialization is needed - g4Label_ = iConfig.getParameter("ModuleLabel"); - ebLabel_ = iConfig.getParameter("EBCollection"); - eeLabel_ = iConfig.getParameter("EECollection"); - hcLabel_ = iConfig.getParameter("HCCollection"); - verbosity_ = iConfig.getUntrackedParameter("Verbosity",0); - maxDepth_ = iConfig.getUntrackedParameter("MaxDepth",4); - etaMax_ = iConfig.getUntrackedParameter("EtaMax", 3.0); - tMinE_ = iConfig.getUntrackedParameter("TimeMinCutECAL", -500.); - tMaxE_ = iConfig.getUntrackedParameter("TimeMaxCutECAL", 500.); - tMinH_ = iConfig.getUntrackedParameter("TimeMinCutHCAL", -500.); - tMaxH_ = iConfig.getUntrackedParameter("TimeMaxCutHCAL", 500.); - - tok_SimTk_ = consumes(edm::InputTag(g4Label_)); + g4Label_ = iConfig.getParameter("ModuleLabel"); + ebLabel_ = iConfig.getParameter("EBCollection"); + eeLabel_ = iConfig.getParameter("EECollection"); + hcLabel_ = iConfig.getParameter("HCCollection"); + verbosity_ = iConfig.getUntrackedParameter("Verbosity", 0); + maxDepth_ = iConfig.getUntrackedParameter("MaxDepth", 4); + etaMax_ = iConfig.getUntrackedParameter("EtaMax", 3.0); + tMinE_ = iConfig.getUntrackedParameter("TimeMinCutECAL", -500.); + tMaxE_ = iConfig.getUntrackedParameter("TimeMaxCutECAL", 500.); + tMinH_ = iConfig.getUntrackedParameter("TimeMinCutHCAL", -500.); + tMaxH_ = iConfig.getUntrackedParameter("TimeMaxCutHCAL", 500.); + + tok_SimTk_ = consumes(edm::InputTag(g4Label_)); tok_SimVtx_ = consumes(edm::InputTag(g4Label_)); - tok_caloEB_ = consumes(edm::InputTag(g4Label_,ebLabel_)); - tok_caloEE_ = consumes(edm::InputTag(g4Label_,eeLabel_)); - tok_caloHH_ = consumes(edm::InputTag(g4Label_,hcLabel_)); - if (maxDepth_ > depthMax_) maxDepth_ = depthMax_; - else if (maxDepth_ < 1) maxDepth_ = 4; + tok_caloEB_ = consumes(edm::InputTag(g4Label_, ebLabel_)); + tok_caloEE_ = consumes(edm::InputTag(g4Label_, eeLabel_)); + tok_caloHH_ = consumes(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; + 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(); tok_hdcons_ = esConsumes(); @@ -132,15 +130,13 @@ HcalHBHEMuonSimAnalyzer::HcalHBHEMuonSimAnalyzer(const edm::ParameterSet& iConfi tok_htopo_ = esConsumes(); } -HcalHBHEMuonSimAnalyzer::~HcalHBHEMuonSimAnalyzer() { } - -void HcalHBHEMuonSimAnalyzer::analyze(const edm::Event& iEvent, - const edm::EventSetup& iSetup) { +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); + 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 @@ -148,17 +144,17 @@ void HcalHBHEMuonSimAnalyzer::analyze(const edm::Event& iEvent, const HcalDDDRecConstants* hcons = &iSetup.getData(tok_hdcons_); - runNumber_ = iEvent.id().run(); + runNumber_ = iEvent.id().run(); eventNumber_ = iEvent.id().event(); - lumiNumber_ = iEvent.id().luminosityBlock(); - bxNumber_ = iEvent.bunchCrossing(); - + lumiNumber_ = iEvent.id().luminosityBlock(); + bxNumber_ = iEvent.bunchCrossing(); + //get Handles to SimTracks and SimHits edm::Handle SimTk; - iEvent.getByToken(tok_SimTk_,SimTk); + iEvent.getByToken(tok_SimTk_, SimTk); edm::SimTrackContainer::const_iterator simTrkItr; edm::Handle SimVtx; - iEvent.getByToken(tok_SimVtx_,SimVtx); + iEvent.getByToken(tok_SimVtx_, SimVtx); //get Handles to PCaloHitContainers of eb/ee/hbhe edm::Handle pcaloeb; @@ -169,39 +165,40 @@ void HcalHBHEMuonSimAnalyzer::analyze(const edm::Event& iEvent, iEvent.getByToken(tok_caloHH_, pcalohh); std::vector calohh; bool testN(false); - for (unsigned int k=1; ksize(); ++k) { + 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;} + 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) { + for (edm::PCaloHitContainer::const_iterator itr = pcalohh->begin(); itr != pcalohh->end(); ++itr) { PCaloHit hit(*itr); - DetId newid = HcalHitRelabeller::relabel(hit.id(),hcons); + 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; + 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()); + 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 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++) { + 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_)) { + (std::abs(simTrkItr->momentum().eta()) < etaMax_)) { unsigned int thisTrk = simTrkItr->trackId(); spr::propagatedTrackDirection trkD = spr::propagateCALO(thisTrk, SimTk, SimVtx, geo, bField, debug); @@ -209,220 +206,207 @@ void HcalHBHEMuonSimAnalyzer::analyze(const edm::Event& iEvent, double eHcalDepth[depthMax_], eHcalDepthHot[depthMax_]; double activeL[depthMax_], activeHotL[depthMax_]; unsigned int isHot(0); - bool tmpmatch(false); - for (int i=0; i 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; + 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_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()); + // 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; + 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); + // 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; + 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 > 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= 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; + HcalSubdetector subdet = HcalDetId(closestCell).subdet(); + int ieta = HcalDetId(closestCell).ieta(); + int iphi = HcalDetId(closestCell).iphi(); + bool hbhe = (std::abs(ieta) == 16); + std::vector > 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; + if ((verbosity_ % 10) > 0) + std::cout << hcid0 << " E " << ehdepth[i].first << " L " << actL << std::endl; #endif - } + } - HcalDetId hotCell; + HcalDetId hotCell; #ifdef EDM_ML_DEBUG - double h3x3 = + double h3x3 = #endif - spr::eHCALmatrix(geo,theHBHETopology, closestCell, calohh, 1,1, hotCell, debug); - isHot = matchId(closestCell,hotCell); + 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; + 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 > ehdepth; - spr::energyHCALCell(hotCell, calohh, ehdepth, maxDepth_, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_, debug); - for (unsigned int i=0; i= 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; + if (hotCell != HcalDetId()) { + subdet = HcalDetId(hotCell).subdet(); + ieta = HcalDetId(hotCell).ieta(); + iphi = HcalDetId(hotCell).iphi(); + hbhe = (std::abs(ieta) == 16); + std::vector > 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; + 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 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(); + if (hcalHot_.size() > 0) + tree_->Fill(); } void HcalHBHEMuonSimAnalyzer::beginJob() { - tree_ = fs_->make("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_); + 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; kBranch(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)); + 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_); +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("ModuleLabel","g4SimHits"); - desc.add("EBCollection","EcalHitsEB"); - desc.add("EECollection","EcalHitsEE"); - desc.add("HCCollection","HcalHits"); - desc.addUntracked("Verbosity",0); - desc.addUntracked("MaxDepth",4); - desc.addUntracked("EtaMax",3.0); - desc.addUntracked("TimeMinCutECAL",-500.); - desc.addUntracked("TimeMaxCutECAL",500.); - desc.addUntracked("TimeMinCutHCAL",-500.); - desc.addUntracked("TimeMaxCutHCAL",500.); - descriptions.add("hcalHBHEMuonSim",desc); + desc.add("ModuleLabel", "g4SimHits"); + desc.add("EBCollection", "EcalHitsEB"); + desc.add("EECollection", "EcalHitsEE"); + desc.add("HCCollection", "HcalHits"); + desc.addUntracked("Verbosity", 0); + desc.addUntracked("MaxDepth", 4); + desc.addUntracked("EtaMax", 3.0); + desc.addUntracked("TimeMinCutECAL", -500.); + desc.addUntracked("TimeMaxCutECAL", 500.); + desc.addUntracked("TimeMinCutHCAL", -500.); + desc.addUntracked("TimeMaxCutHCAL", 500.); + descriptions.add("hcalHBHEMuonSim", desc); } void HcalHBHEMuonSimAnalyzer::clearVectors() { - ///clearing vectots - runNumber_ = -99999; + runNumber_ = -99999; eventNumber_ = -99999; - lumiNumber_ = -99999; - bxNumber_ = -99999; - + lumiNumber_ = -99999; + bxNumber_ = -99999; + ptGlob_.clear(); - etaGlob_.clear(); - phiGlob_.clear(); + etaGlob_.clear(); + phiGlob_.clear(); pMuon_.clear(); matchedId_.clear(); ecal3x3Energy_.clear(); @@ -432,7 +416,7 @@ void HcalHBHEMuonSimAnalyzer::clearVectors() { hcalHot_.clear(); hcalActiveLength_.clear(); hcalActiveLengthHot_.clear(); - for (int k=0; k Date: Mon, 24 May 2021 15:09:12 +0200 Subject: [PATCH 3/3] Correct for Clang error --- Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc b/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc index f87e908dbc204..684111e5f2484 100644 --- a/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc +++ b/Calibration/HcalCalibAlgos/test/HcalHBHEMuonSimAnalyzer.cc @@ -38,7 +38,7 @@ #include "MagneticField/Engine/interface/MagneticField.h" #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" -#define EDM_ML_DEBUG +//#define EDM_ML_DEBUG class HcalHBHEMuonSimAnalyzer : public edm::one::EDAnalyzer { public: @@ -152,7 +152,6 @@ void HcalHBHEMuonSimAnalyzer::analyze(const edm::Event& iEvent, const edm::Event //get Handles to SimTracks and SimHits edm::Handle SimTk; iEvent.getByToken(tok_SimTk_, SimTk); - edm::SimTrackContainer::const_iterator simTrkItr; edm::Handle SimVtx; iEvent.getByToken(tok_SimVtx_, SimVtx);