From 8011580c9c1b09008a99f2750595bbbaa13cfe0a Mon Sep 17 00:00:00 2001 From: Jean-Baptiste Sauvan Date: Thu, 12 Oct 2017 17:05:44 +0200 Subject: [PATCH 1/7] Add sim energy in TC ntuple --- .../HGCalTriggerNtupleHGCTriggerCells.cc | 118 +++++++++++++++++- .../python/hgcalTriggerNtuples_cfi.py | 6 +- 2 files changed, 121 insertions(+), 3 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc index f6c82c895dea1..c4491b33703f3 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc @@ -1,4 +1,8 @@ +#include "SimDataFormats/CaloHit/interface/PCaloHit.h" +#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" +#include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h" +#include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h" #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h" #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" @@ -7,6 +11,7 @@ + class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase { @@ -17,10 +22,15 @@ class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase virtual void fill(const edm::Event& e, const edm::EventSetup& es) override final; private: + void simhits(const edm::Event& e, std::unordered_map& simhits_ee, std::unordered_map& simhits_fh, std::unordered_map& simhits_bh); virtual void clear() override final; edm::EDGetToken trigger_cells_token_; + edm::EDGetToken simhits_ee_token_, simhits_fh_token_, simhits_bh_token_; + bool fill_simenergy_; + edm::ESHandle geometry_; + int tc_n_ ; std::vector tc_id_; @@ -33,6 +43,7 @@ class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase std::vector tc_data_; std::vector tc_mipPt_; std::vector tc_energy_; + std::vector tc_simenergy_; std::vector tc_eta_; std::vector tc_phi_; std::vector tc_pt_; @@ -48,6 +59,7 @@ DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, HGCalTriggerNtupleHGCTriggerCells:: HGCalTriggerNtupleHGCTriggerCells(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf) { + fill_simenergy_ = conf.getParameter("FillSimEnergy"); } void @@ -56,6 +68,13 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& { trigger_cells_token_ = collector.consumes(conf.getParameter("TriggerCells")); + if (fill_simenergy_) + { + simhits_ee_token_ = collector.consumes(conf.getParameter("eeSimHits")); + simhits_fh_token_ = collector.consumes(conf.getParameter("fhSimHits")); + simhits_bh_token_ = collector.consumes(conf.getParameter("bhSimHits")); + } + tree.Branch("tc_n", &tc_n_, "tc_n/I"); tree.Branch("tc_id", &tc_id_); tree.Branch("tc_subdet", &tc_subdet_); @@ -67,6 +86,7 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("tc_data", &tc_data_); tree.Branch("tc_mipPt", &tc_mipPt_); tree.Branch("tc_energy", &tc_energy_); + if(fill_simenergy_) tree.Branch("tc_simenergy", &tc_simenergy_); tree.Branch("tc_eta", &tc_eta_); tree.Branch("tc_phi", &tc_phi_); tree.Branch("tc_pt", &tc_pt_); @@ -85,8 +105,13 @@ fill(const edm::Event& e, const edm::EventSetup& es) const l1t::HGCalTriggerCellBxCollection& trigger_cells = *trigger_cells_h; // retrieve geometry - edm::ESHandle geometry; - es.get().get(geometry); + es.get().get(geometry_); + + // sim hit association + std::unordered_map simhits_ee; + std::unordered_map simhits_fh; + std::unordered_map simhits_bh; + if(fill_simenergy_) simhits(e, simhits_ee, simhits_fh, simhits_bh); clear(); for(auto tc_itr=trigger_cells.begin(0); tc_itr!=trigger_cells.end(0); tc_itr++) @@ -111,11 +136,99 @@ fill(const edm::Event& e, const edm::EventSetup& es) tc_phi_.emplace_back(tc_itr->phi()); tc_pt_.emplace_back(tc_itr->pt()); tc_z_.emplace_back(tc_itr->position().z()); + + if(fill_simenergy_) + { + double energy = 0; + int subdet = id.subdetId(); + // search for simhit for all the cells inside the trigger cell + for(uint32_t c_id : geometry_->getCellsFromTriggerCell(id)) + { + switch(subdet) + { + case ForwardSubdetector::HGCEE: + { + auto itr = simhits_ee.find(c_id); + if(itr!=simhits_ee.end()) energy += itr->second; + break; + } + case ForwardSubdetector::HGCHEF: + { + auto itr = simhits_fh.find(c_id); + if(itr!=simhits_fh.end()) energy += itr->second; + break; + } + case ForwardSubdetector::HGCHEB: + { + auto itr = simhits_bh.find(c_id); + if(itr!=simhits_bh.end()) energy += itr->second; + break; + } + default: + break; + } + } + tc_simenergy_.emplace_back(energy); + } } } } +void +HGCalTriggerNtupleHGCTriggerCells:: +simhits(const edm::Event& e, std::unordered_map& simhits_ee, std::unordered_map& simhits_fh, std::unordered_map& simhits_bh) +{ + edm::Handle ee_simhits_h; + e.getByToken(simhits_ee_token_,ee_simhits_h); + const edm::PCaloHitContainer& ee_simhits = *ee_simhits_h; + edm::Handle fh_simhits_h; + e.getByToken(simhits_fh_token_,fh_simhits_h); + const edm::PCaloHitContainer& fh_simhits = *fh_simhits_h; + edm::Handle bh_simhits_h; + e.getByToken(simhits_bh_token_,bh_simhits_h); + const edm::PCaloHitContainer& bh_simhits = *bh_simhits_h; + + //EE + int layer=0,cell=0, sec=0, subsec=0, zp=0,subdet=0; + ForwardSubdetector mysubdet; + for( const auto& simhit : ee_simhits ) + { + HGCalTestNumbering::unpackHexagonIndex(simhit.id(), subdet, zp, layer, sec, subsec, cell); + mysubdet = (ForwardSubdetector)subdet; + std::pair recoLayerCell = geometry_->eeTopology().dddConstants().simToReco(cell,layer,sec,geometry_->eeTopology().detectorType()); + cell = recoLayerCell.first; + layer = recoLayerCell.second; + if (layer<0 || cell<0) continue; + auto itr_insert = simhits_ee.emplace(HGCalDetId(mysubdet,zp,layer,subsec,sec,cell), 0.); + itr_insert.first->second += simhit.energy(); + } + + // FH + layer=0; cell=0; sec=0; subsec=0; zp=0; subdet=0; + for( const auto& simhit : fh_simhits ) + { + HGCalTestNumbering::unpackHexagonIndex(simhit.id(), subdet, zp, layer, sec, subsec, cell); + mysubdet = (ForwardSubdetector)(subdet); + std::pair recoLayerCell = geometry_->fhTopology().dddConstants().simToReco(cell,layer,sec,geometry_->fhTopology().detectorType()); + cell = recoLayerCell.first; + layer = recoLayerCell.second; + if (layer<0 || cell<0) continue; + auto itr_insert = simhits_fh.emplace(HGCalDetId(mysubdet,zp,layer,subsec,sec,cell), 0.); + itr_insert.first->second += simhit.energy(); + } + + // BH + for( const auto& simhit : bh_simhits ) + { + HcalDetId id = HcalHitRelabeller::relabel(simhit.id(), geometry_->bhTopology().dddConstants()); + if (id.subdetId()!=HcalEndcap) continue; + auto itr_insert = simhits_bh.emplace(id, 0.); + itr_insert.first->second += simhit.energy(); + } +} + + void HGCalTriggerNtupleHGCTriggerCells:: clear() @@ -131,6 +244,7 @@ clear() tc_data_.clear(); tc_mipPt_.clear(); tc_energy_.clear(); + tc_simenergy_.clear(); tc_eta_.clear(); tc_phi_.clear(); tc_pt_.clear(); diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py b/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py index ac0369628eef6..5224555555866 100644 --- a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py @@ -35,7 +35,11 @@ ntuple_triggercells = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleHGCTriggerCells'), - TriggerCells = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:calibratedTriggerCells') + TriggerCells = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:calibratedTriggerCells'), + FillSimEnergy = cms.bool(False), + eeSimHits = cms.InputTag('g4SimHits:HGCHitsEE'), + fhSimHits = cms.InputTag('g4SimHits:HGCHitsHEfront'), + bhSimHits = cms.InputTag('g4SimHits:HcalHits') ) ntuple_clusters = cms.PSet( From 028aaa05e283832c0703ab7ff29d629cc7770221 Mon Sep 17 00:00:00 2001 From: Jean-Baptiste Sauvan Date: Fri, 13 Oct 2017 09:50:41 +0200 Subject: [PATCH 2/7] Add links between trigger cells and clusters --- .../ntuples/HGCalTriggerNtupleHGCClusters.cc | 51 +++++++++++++++---- .../HGCalTriggerNtupleHGCMulticlusters.cc | 22 ++++---- .../HGCalTriggerNtupleHGCTriggerCells.cc | 48 ++++++++++++++++- .../python/hgcalTriggerNtuples_cfi.py | 6 ++- 4 files changed, 103 insertions(+), 24 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc index f0fdd41ddbe66..57408ff6d516e 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc @@ -1,6 +1,7 @@ #include #include "DataFormats/L1THGCal/interface/HGCalCluster.h" +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" @@ -20,8 +21,8 @@ class HGCalTriggerNtupleHGCClusters : public HGCalTriggerNtupleBase private: virtual void clear() override final; - - edm::EDGetToken clusters_token_; + bool filter_clusters_; + edm::EDGetToken clusters_token_, multiclusters_token_; int cl_n_ ; std::vector cl_mipPt_; @@ -30,8 +31,10 @@ class HGCalTriggerNtupleHGCClusters : public HGCalTriggerNtupleBase std::vector cl_eta_; std::vector cl_phi_; std::vector cl_layer_; - std::vector cl_ncells_; - std::vector> cl_cells_; + std::vector cl_cells_n_; + std::vector> cl_cells_id_; + std::vector cl_multicluster_id_; + std::vector cl_multicluster_pt_; }; @@ -43,6 +46,7 @@ DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, HGCalTriggerNtupleHGCClusters:: HGCalTriggerNtupleHGCClusters(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf) { + filter_clusters_ = conf.getParameter("FilterClusters"); } void @@ -50,6 +54,7 @@ HGCalTriggerNtupleHGCClusters:: initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& collector) { clusters_token_ = collector.consumes(conf.getParameter("Clusters")); + multiclusters_token_ = collector.consumes(conf.getParameter("Multiclusters")); tree.Branch("cl_n", &cl_n_, "cl_n/I"); tree.Branch("cl_mipPt", &cl_mipPt_); @@ -58,8 +63,10 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("cl_eta", &cl_eta_); tree.Branch("cl_phi", &cl_phi_); tree.Branch("cl_layer", &cl_layer_); - tree.Branch("cl_ncells", &cl_ncells_); - tree.Branch("cl_cells", &cl_cells_); + tree.Branch("cl_cells_n", &cl_cells_n_); + tree.Branch("cl_cells_id", &cl_cells_id_); + tree.Branch("cl_multicluster_id", &cl_multicluster_id_); + tree.Branch("cl_multicluster_pt", &cl_multicluster_pt_); } void @@ -71,14 +78,32 @@ fill(const edm::Event& e, const edm::EventSetup& es) edm::Handle clusters_h; e.getByToken(clusters_token_, clusters_h); const l1t::HGCalClusterBxCollection& clusters = *clusters_h; + edm::Handle multiclusters_h; + e.getByToken(multiclusters_token_, multiclusters_h); + const l1t::HGCalMulticlusterBxCollection& multiclusters = *multiclusters_h; // retrieve geometry edm::ESHandle geometry; es.get().get(geometry); + // Associate cells to clusters + std::unordered_map cluster2multicluster; + for(auto mcl_itr=multiclusters.begin(0); mcl_itr!=multiclusters.end(0); mcl_itr++) + { + // loop on 2D clusters inside 3D clusters + for(const auto& cl_ptr : mcl_itr->constituents()) + { + cluster2multicluster.emplace(cl_ptr->detId(), mcl_itr); + } + } + clear(); for(auto cl_itr=clusters.begin(0); cl_itr!=clusters.end(0); cl_itr++) { + auto mcl_itr = cluster2multicluster.find(cl_itr->detId()); + uint32_t mcl_id = (mcl_itr!=cluster2multicluster.end() ? mcl_itr->second->detId() : 0); + float mcl_pt = (mcl_itr!=cluster2multicluster.end() ? mcl_itr->second->pt() : 0.); + if(filter_clusters_ && mcl_id==0) continue; cl_n_++; cl_mipPt_.emplace_back(cl_itr->mipPt()); // physical values @@ -87,12 +112,14 @@ fill(const edm::Event& e, const edm::EventSetup& es) cl_eta_.emplace_back(cl_itr->eta()); cl_phi_.emplace_back(cl_itr->phi()); cl_layer_.emplace_back(cl_itr->layer()); - cl_ncells_.emplace_back(cl_itr->constituents().size()); + cl_cells_n_.emplace_back(cl_itr->constituents().size()); // Retrieve indices of trigger cells inside cluster - cl_cells_.emplace_back(cl_itr->constituents().size()); + cl_cells_id_.emplace_back(cl_itr->constituents().size()); std::transform(cl_itr->constituents_begin(), cl_itr->constituents_end(), - cl_cells_.back().begin(), [](const edm::Ptr& tc){return tc.key();} + cl_cells_id_.back().begin(), [](const edm::Ptr& tc){return tc->detId();} ); + cl_multicluster_id_.emplace_back(mcl_id); + cl_multicluster_pt_.emplace_back(mcl_pt); } } @@ -107,8 +134,10 @@ clear() cl_eta_.clear(); cl_phi_.clear(); cl_layer_.clear(); - cl_ncells_.clear(); - cl_cells_.clear(); + cl_cells_n_.clear(); + cl_cells_id_.clear(); + cl_multicluster_id_.clear(); + cl_multicluster_pt_.clear(); } diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc index a0da63fbb63a3..61bab49aba7ea 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc @@ -25,7 +25,9 @@ class HGCalTriggerNtupleHGCMulticlusters : public HGCalTriggerNtupleBase std::vector cl3d_energy_; std::vector cl3d_eta_; std::vector cl3d_phi_; - std::vector cl3d_nclu_; + std::vector cl3d_clusters_n_; + std::vector> cl3d_clusters_id_; + // cluster shower shapes std::vector cl3d_showerlength_; std::vector cl3d_firstlayer_; std::vector cl3d_seetot_; @@ -36,7 +38,6 @@ class HGCalTriggerNtupleHGCMulticlusters : public HGCalTriggerNtupleBase std::vector cl3d_srrtot_; std::vector cl3d_srrmax_; std::vector cl3d_emaxe_; - std::vector> cl3d_clusters_; }; DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, @@ -60,7 +61,8 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("cl3d_energy", &cl3d_energy_); tree.Branch("cl3d_eta", &cl3d_eta_); tree.Branch("cl3d_phi", &cl3d_phi_); - tree.Branch("cl3d_nclu", &cl3d_nclu_); + tree.Branch("cl3d_clusters_n", &cl3d_clusters_n_); + tree.Branch("cl3d_clusters_id", &cl3d_clusters_id_); tree.Branch("cl3d_showerlength", &cl3d_showerlength_); tree.Branch("cl3d_firstlayer", &cl3d_firstlayer_); tree.Branch("cl3d_seetot", &cl3d_seetot_); @@ -70,8 +72,7 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("cl3d_szz", &cl3d_szz_); tree.Branch("cl3d_srrtot", &cl3d_srrtot_); tree.Branch("cl3d_srrmax", &cl3d_srrmax_); - tree.Branch("cl3d_emaxe", &cl3d_emaxe_); - tree.Branch("cl3d_clusters", &cl3d_clusters_); + tree.Branch("cl3d_emaxe", &cl3d_emaxe_); } @@ -98,7 +99,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) cl3d_energy_.emplace_back(cl3d_itr->energy()); cl3d_eta_.emplace_back(cl3d_itr->eta()); cl3d_phi_.emplace_back(cl3d_itr->phi()); - cl3d_nclu_.emplace_back(cl3d_itr->constituents().size()); + cl3d_clusters_n_.emplace_back(cl3d_itr->constituents().size()); cl3d_showerlength_.emplace_back(cl3d_itr->showerLength()); cl3d_firstlayer_.emplace_back(cl3d_itr->firstLayer()); cl3d_seetot_.emplace_back(cl3d_itr->sigmaEtaEtaTot()); @@ -111,9 +112,9 @@ fill(const edm::Event& e, const edm::EventSetup& es) cl3d_emaxe_.emplace_back(cl3d_itr->eMax()/cl3d_itr->energy()); // Retrieve indices of trigger cells inside cluster - cl3d_clusters_.emplace_back(cl3d_itr->constituents().size()); + cl3d_clusters_id_.emplace_back(cl3d_itr->constituents().size()); std::transform(cl3d_itr->constituents_begin(), cl3d_itr->constituents_end(), - cl3d_clusters_.back().begin(), [](const edm::Ptr& cl){return cl.key();} + cl3d_clusters_id_.back().begin(), [](const edm::Ptr& cl){return cl->detId();} ); } } @@ -128,7 +129,8 @@ clear() cl3d_energy_.clear(); cl3d_eta_.clear(); cl3d_phi_.clear(); - cl3d_nclu_.clear(); + cl3d_clusters_n_.clear(); + cl3d_clusters_id_.clear(); cl3d_showerlength_.clear(); cl3d_firstlayer_.clear(); cl3d_seetot_.clear(); @@ -139,8 +141,6 @@ clear() cl3d_srrtot_.clear(); cl3d_srrmax_.clear(); cl3d_emaxe_.clear(); - cl3d_clusters_.clear(); - } diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc index c4491b33703f3..17e544270f988 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc @@ -4,6 +4,7 @@ #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h" #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h" #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h" +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" @@ -26,9 +27,10 @@ class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase virtual void clear() override final; - edm::EDGetToken trigger_cells_token_; + edm::EDGetToken trigger_cells_token_, multiclusters_token_; edm::EDGetToken simhits_ee_token_, simhits_fh_token_, simhits_bh_token_; bool fill_simenergy_; + bool filter_cluster_cells_; edm::ESHandle geometry_; @@ -48,6 +50,9 @@ class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase std::vector tc_phi_; std::vector tc_pt_; std::vector tc_z_; + std::vector tc_cluster_id_; + std::vector tc_multicluster_id_; + std::vector tc_multicluster_pt_; }; @@ -60,6 +65,7 @@ HGCalTriggerNtupleHGCTriggerCells:: HGCalTriggerNtupleHGCTriggerCells(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf) { fill_simenergy_ = conf.getParameter("FillSimEnergy"); + filter_cluster_cells_ = conf.getParameter("FilterClusterCells"); } void @@ -67,6 +73,7 @@ HGCalTriggerNtupleHGCTriggerCells:: initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& collector) { trigger_cells_token_ = collector.consumes(conf.getParameter("TriggerCells")); + multiclusters_token_ = collector.consumes(conf.getParameter("Multiclusters")); if (fill_simenergy_) { @@ -91,6 +98,9 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("tc_phi", &tc_phi_); tree.Branch("tc_pt", &tc_pt_); tree.Branch("tc_z", &tc_z_); + tree.Branch("tc_cluster_id", &tc_cluster_id_); + tree.Branch("tc_multicluster_id", &tc_multicluster_id_); + tree.Branch("tc_multicluster_pt", &tc_multicluster_pt_); } @@ -104,6 +114,11 @@ fill(const edm::Event& e, const edm::EventSetup& es) e.getByToken(trigger_cells_token_, trigger_cells_h); const l1t::HGCalTriggerCellBxCollection& trigger_cells = *trigger_cells_h; + // retrieve clusters + edm::Handle multiclusters_h; + e.getByToken(multiclusters_token_, multiclusters_h); + const l1t::HGCalMulticlusterBxCollection& multiclusters = *multiclusters_h; + // retrieve geometry es.get().get(geometry_); @@ -113,11 +128,35 @@ fill(const edm::Event& e, const edm::EventSetup& es) std::unordered_map simhits_bh; if(fill_simenergy_) simhits(e, simhits_ee, simhits_fh, simhits_bh); + // Associate cells to clusters + std::unordered_map cell2cluster; + std::unordered_map cell2multicluster; + for(auto mcl_itr=multiclusters.begin(0); mcl_itr!=multiclusters.end(0); mcl_itr++) + { + // loop on 2D clusters inside 3D clusters + for(const auto& cl_ptr : mcl_itr->constituents()) + { + // loop on TC inside 2D clusters + for(const auto& tc_ptr : cl_ptr->constituents()) + { + cell2cluster.emplace(tc_ptr->detId(), cl_ptr->detId()); + cell2multicluster.emplace(tc_ptr->detId(), mcl_itr); + } + } + } + clear(); for(auto tc_itr=trigger_cells.begin(0); tc_itr!=trigger_cells.end(0); tc_itr++) { if(tc_itr->hwPt()>0) { + auto cl_itr = cell2cluster.find(tc_itr->detId()); + auto mcl_itr = cell2multicluster.find(tc_itr->detId()); + uint32_t cl_id = (cl_itr!=cell2cluster.end() ? cl_itr->second : 0); + uint32_t mcl_id = (mcl_itr!=cell2multicluster.end() ? mcl_itr->second->detId() : 0); + float mcl_pt = (mcl_itr!=cell2multicluster.end() ? mcl_itr->second->pt() : 0.); + // Filter cells not included in a cluster, if requested + if(filter_cluster_cells_ && mcl_id==0) continue; tc_n_++; // hardware data HGCalDetId id(tc_itr->detId()); @@ -136,6 +175,10 @@ fill(const edm::Event& e, const edm::EventSetup& es) tc_phi_.emplace_back(tc_itr->phi()); tc_pt_.emplace_back(tc_itr->pt()); tc_z_.emplace_back(tc_itr->position().z()); + // Links between TC and clusters + tc_cluster_id_.emplace_back(cl_id); + tc_multicluster_id_.emplace_back(mcl_id); + tc_multicluster_pt_.emplace_back(mcl_pt); if(fill_simenergy_) { @@ -249,6 +292,9 @@ clear() tc_phi_.clear(); tc_pt_.clear(); tc_z_.clear(); + tc_cluster_id_.clear(); + tc_multicluster_id_.clear(); + tc_multicluster_pt_.clear(); } diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py b/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py index 5224555555866..436681e4f86fc 100644 --- a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py @@ -36,7 +36,9 @@ ntuple_triggercells = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleHGCTriggerCells'), TriggerCells = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:calibratedTriggerCells'), + Multiclusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster3D'), FillSimEnergy = cms.bool(False), + FilterClusterCells = cms.bool(True), eeSimHits = cms.InputTag('g4SimHits:HGCHitsEE'), fhSimHits = cms.InputTag('g4SimHits:HGCHitsHEfront'), bhSimHits = cms.InputTag('g4SimHits:HcalHits') @@ -44,7 +46,9 @@ ntuple_clusters = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleHGCClusters'), - Clusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster2D') + Clusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster2D'), + Multiclusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster3D'), + FilterClusters = cms.bool(True) ) ntuple_multicluster = cms.PSet( From 457f149b04b49fbdc4541f05144c967b3e04affc Mon Sep 17 00:00:00 2001 From: Jean-Baptiste Sauvan Date: Mon, 30 Oct 2017 13:50:14 +0100 Subject: [PATCH 3/7] Add variables in ntuples Conflicts: L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc --- .../plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc | 8 ++++++++ .../plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc | 4 ++++ .../plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc | 8 ++++++++ 3 files changed, 20 insertions(+) diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc index 57408ff6d516e..bbde5c0a9b625 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc @@ -26,11 +26,13 @@ class HGCalTriggerNtupleHGCClusters : public HGCalTriggerNtupleBase int cl_n_ ; std::vector cl_mipPt_; + std::vector cl_id_; std::vector cl_pt_; std::vector cl_energy_; std::vector cl_eta_; std::vector cl_phi_; std::vector cl_layer_; + std::vector cl_subdet_; std::vector cl_cells_n_; std::vector> cl_cells_id_; std::vector cl_multicluster_id_; @@ -58,11 +60,13 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("cl_n", &cl_n_, "cl_n/I"); tree.Branch("cl_mipPt", &cl_mipPt_); + tree.Branch("cl_id", &cl_id_); tree.Branch("cl_pt", &cl_pt_); tree.Branch("cl_energy", &cl_energy_); tree.Branch("cl_eta", &cl_eta_); tree.Branch("cl_phi", &cl_phi_); tree.Branch("cl_layer", &cl_layer_); + tree.Branch("cl_subdet", &cl_subdet_); tree.Branch("cl_cells_n", &cl_cells_n_); tree.Branch("cl_cells_id", &cl_cells_id_); tree.Branch("cl_multicluster_id", &cl_multicluster_id_); @@ -111,7 +115,9 @@ fill(const edm::Event& e, const edm::EventSetup& es) cl_energy_.emplace_back(cl_itr->energy()); cl_eta_.emplace_back(cl_itr->eta()); cl_phi_.emplace_back(cl_itr->phi()); + cl_id_.emplace_back(cl_itr->detId()); cl_layer_.emplace_back(cl_itr->layer()); + cl_subdet_.emplace_back(cl_itr->subdetId()); cl_cells_n_.emplace_back(cl_itr->constituents().size()); // Retrieve indices of trigger cells inside cluster cl_cells_id_.emplace_back(cl_itr->constituents().size()); @@ -129,11 +135,13 @@ HGCalTriggerNtupleHGCClusters:: clear() { cl_n_ = 0; + cl_id_.clear(); cl_pt_.clear(); cl_energy_.clear(); cl_eta_.clear(); cl_phi_.clear(); cl_layer_.clear(); + cl_subdet_.clear(); cl_cells_n_.clear(); cl_cells_id_.clear(); cl_multicluster_id_.clear(); diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc index 61bab49aba7ea..0b95f3265f8f1 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc @@ -21,6 +21,7 @@ class HGCalTriggerNtupleHGCMulticlusters : public HGCalTriggerNtupleBase edm::EDGetToken multiclusters_token_; int cl3d_n_ ; + std::vector cl3d_id_; std::vector cl3d_pt_; std::vector cl3d_energy_; std::vector cl3d_eta_; @@ -57,6 +58,7 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& multiclusters_token_ = collector.consumes(conf.getParameter("Multiclusters")); tree.Branch("cl3d_n", &cl3d_n_, "cl3d_n/I"); + tree.Branch("cl3d_id", &cl3d_id_); tree.Branch("cl3d_pt", &cl3d_pt_); tree.Branch("cl3d_energy", &cl3d_energy_); tree.Branch("cl3d_eta", &cl3d_eta_); @@ -94,6 +96,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) for(auto cl3d_itr=multiclusters.begin(0); cl3d_itr!=multiclusters.end(0); cl3d_itr++) { cl3d_n_++; + cl3d_id_.emplace_back(cl3d_itr->detId()); // physical values cl3d_pt_.emplace_back(cl3d_itr->pt()); cl3d_energy_.emplace_back(cl3d_itr->energy()); @@ -125,6 +128,7 @@ HGCalTriggerNtupleHGCMulticlusters:: clear() { cl3d_n_ = 0; + cl3d_id_.clear(); cl3d_pt_.clear(); cl3d_energy_.clear(); cl3d_eta_.clear(); diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc index 17e544270f988..47e1a332870d4 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc @@ -49,6 +49,8 @@ class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase std::vector tc_eta_; std::vector tc_phi_; std::vector tc_pt_; + std::vector tc_x_; + std::vector tc_y_; std::vector tc_z_; std::vector tc_cluster_id_; std::vector tc_multicluster_id_; @@ -97,6 +99,8 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("tc_eta", &tc_eta_); tree.Branch("tc_phi", &tc_phi_); tree.Branch("tc_pt", &tc_pt_); + tree.Branch("tc_x", &tc_x_); + tree.Branch("tc_y", &tc_y_); tree.Branch("tc_z", &tc_z_); tree.Branch("tc_cluster_id", &tc_cluster_id_); tree.Branch("tc_multicluster_id", &tc_multicluster_id_); @@ -174,6 +178,8 @@ fill(const edm::Event& e, const edm::EventSetup& es) tc_eta_.emplace_back(tc_itr->eta()); tc_phi_.emplace_back(tc_itr->phi()); tc_pt_.emplace_back(tc_itr->pt()); + tc_x_.emplace_back(tc_itr->position().x()); + tc_y_.emplace_back(tc_itr->position().y()); tc_z_.emplace_back(tc_itr->position().z()); // Links between TC and clusters tc_cluster_id_.emplace_back(cl_id); @@ -291,6 +297,8 @@ clear() tc_eta_.clear(); tc_phi_.clear(); tc_pt_.clear(); + tc_x_.clear(); + tc_y_.clear(); tc_z_.clear(); tc_cluster_id_.clear(); tc_multicluster_id_.clear(); From 770e9c53f9581dca66046dd1f9aab0826ad3b794 Mon Sep 17 00:00:00 2001 From: Jean-Baptiste Sauvan Date: Mon, 30 Oct 2017 14:07:44 +0100 Subject: [PATCH 4/7] Change parameter names --- .../plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc | 6 +++--- .../plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc | 8 ++++---- L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py | 8 ++++---- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc index bbde5c0a9b625..6f8a279d0ec35 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc @@ -21,7 +21,7 @@ class HGCalTriggerNtupleHGCClusters : public HGCalTriggerNtupleBase private: virtual void clear() override final; - bool filter_clusters_; + bool filter_clusters_in_multiclusters_; edm::EDGetToken clusters_token_, multiclusters_token_; int cl_n_ ; @@ -48,7 +48,7 @@ DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, HGCalTriggerNtupleHGCClusters:: HGCalTriggerNtupleHGCClusters(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf) { - filter_clusters_ = conf.getParameter("FilterClusters"); + filter_clusters_in_multiclusters_ = conf.getParameter("FilterClustersInMulticlusters"); } void @@ -107,7 +107,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) auto mcl_itr = cluster2multicluster.find(cl_itr->detId()); uint32_t mcl_id = (mcl_itr!=cluster2multicluster.end() ? mcl_itr->second->detId() : 0); float mcl_pt = (mcl_itr!=cluster2multicluster.end() ? mcl_itr->second->pt() : 0.); - if(filter_clusters_ && mcl_id==0) continue; + if(filter_clusters_in_multiclusters_ && mcl_id==0) continue; cl_n_++; cl_mipPt_.emplace_back(cl_itr->mipPt()); // physical values diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc index 47e1a332870d4..1804f4a303038 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc @@ -30,7 +30,7 @@ class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase edm::EDGetToken trigger_cells_token_, multiclusters_token_; edm::EDGetToken simhits_ee_token_, simhits_fh_token_, simhits_bh_token_; bool fill_simenergy_; - bool filter_cluster_cells_; + bool filter_cells_in_multiclusters_; edm::ESHandle geometry_; @@ -67,7 +67,7 @@ HGCalTriggerNtupleHGCTriggerCells:: HGCalTriggerNtupleHGCTriggerCells(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf) { fill_simenergy_ = conf.getParameter("FillSimEnergy"); - filter_cluster_cells_ = conf.getParameter("FilterClusterCells"); + filter_cells_in_multiclusters_ = conf.getParameter("FilterCellsInMulticlusters"); } void @@ -159,8 +159,8 @@ fill(const edm::Event& e, const edm::EventSetup& es) uint32_t cl_id = (cl_itr!=cell2cluster.end() ? cl_itr->second : 0); uint32_t mcl_id = (mcl_itr!=cell2multicluster.end() ? mcl_itr->second->detId() : 0); float mcl_pt = (mcl_itr!=cell2multicluster.end() ? mcl_itr->second->pt() : 0.); - // Filter cells not included in a cluster, if requested - if(filter_cluster_cells_ && mcl_id==0) continue; + // Filter cells not included in a multicluster, if requested + if(filter_cells_in_multiclusters_ && mcl_id==0) continue; tc_n_++; // hardware data HGCalDetId id(tc_itr->detId()); diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py b/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py index 436681e4f86fc..78bbbeaf71587 100644 --- a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py @@ -37,18 +37,18 @@ NtupleName = cms.string('HGCalTriggerNtupleHGCTriggerCells'), TriggerCells = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:calibratedTriggerCells'), Multiclusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster3D'), - FillSimEnergy = cms.bool(False), - FilterClusterCells = cms.bool(True), eeSimHits = cms.InputTag('g4SimHits:HGCHitsEE'), fhSimHits = cms.InputTag('g4SimHits:HGCHitsHEfront'), - bhSimHits = cms.InputTag('g4SimHits:HcalHits') + bhSimHits = cms.InputTag('g4SimHits:HcalHits'), + FillSimEnergy = cms.bool(False), + FilterCellsInMulticlusters = cms.bool(True) ) ntuple_clusters = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleHGCClusters'), Clusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster2D'), Multiclusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster3D'), - FilterClusters = cms.bool(True) + FilterClustersInMulticlusters = cms.bool(True) ) ntuple_multicluster = cms.PSet( From e45003866714d45f296c346e739adbde8455a080 Mon Sep 17 00:00:00 2001 From: Jean-Baptiste Sauvan Date: Mon, 30 Oct 2017 14:25:41 +0100 Subject: [PATCH 5/7] Fix ntuple variable clear --- .../L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc index 6f8a279d0ec35..e8e5201586b77 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc @@ -136,6 +136,7 @@ clear() { cl_n_ = 0; cl_id_.clear(); + cl_mipPt_.clear(); cl_pt_.clear(); cl_energy_.clear(); cl_eta_.clear(); From b33ebce12ac8f965e44bd3b8cc525dcdedb792b2 Mon Sep 17 00:00:00 2001 From: Jean-Baptiste Sauvan Date: Mon, 30 Oct 2017 21:24:22 +0100 Subject: [PATCH 6/7] Add new shower shape variables for e/g id --- .../L1THGCal/interface/HGCalClusterT.h | 9 ++ DataFormats/L1THGCal/src/classes_def.xml | 6 +- .../be_algorithms/HGCalShowerShape.h | 4 + .../HGCalTriggerNtupleHGCMulticlusters.cc | 12 +++ .../be_algorithms/HGCalMulticlusteringImpl.cc | 4 +- .../src/be_algorithms/HGCalShowerShape.cc | 90 ++++++++++++++++++- 6 files changed, 120 insertions(+), 5 deletions(-) diff --git a/DataFormats/L1THGCal/interface/HGCalClusterT.h b/DataFormats/L1THGCal/interface/HGCalClusterT.h index 93d949b527c87..7bd04e79f4540 100644 --- a/DataFormats/L1THGCal/interface/HGCalClusterT.h +++ b/DataFormats/L1THGCal/interface/HGCalClusterT.h @@ -145,7 +145,9 @@ namespace l1t //shower shape int showerLength() const { return showerLength_; } + int coreShowerLength() const { return coreShowerLength_; } int firstLayer() const { return firstLayer_; } + int maxLayer() const { return maxLayer_; } float eMax() const { return eMax_; } float sigmaEtaEtaMax() const { return sigmaEtaEtaMax_; } float sigmaPhiPhiMax() const { return sigmaPhiPhiMax_; } @@ -154,9 +156,12 @@ namespace l1t float sigmaZZ() const { return sigmaZZ_; } float sigmaRRTot() const { return sigmaRRTot_; } float sigmaRRMax() const { return sigmaRRMax_; } + float sigmaRRMean() const { return sigmaRRMean_; } void set_showerLength(int showerLength) { showerLength_ = showerLength;} + void set_coreShowerLength(int coreShowerLength) { coreShowerLength_ = coreShowerLength;} void set_firstLayer(int firstLayer) { firstLayer_ = firstLayer;} + void set_maxLayer(int maxLayer) { maxLayer_ = maxLayer;} void set_eMax(float eMax) { eMax_ = eMax;} void set_sigmaEtaEtaMax(float sigmaEtaEtaMax) { sigmaEtaEtaMax_ = sigmaEtaEtaMax;} void set_sigmaEtaEtaTot(float sigmaEtaEtaTot) { sigmaEtaEtaTot_ = sigmaEtaEtaTot;} @@ -164,6 +169,7 @@ namespace l1t void set_sigmaPhiPhiTot(float sigmaPhiPhiTot) { sigmaPhiPhiTot_ = sigmaPhiPhiTot;} void set_sigmaRRMax(float sigmaRRMax) { sigmaRRMax_ = sigmaRRMax;} void set_sigmaRRTot(float sigmaRRTot) { sigmaRRTot_ = sigmaRRTot;} + void set_sigmaRRMean(float sigmaRRMean) { sigmaRRMean_ = sigmaRRMean;} void set_sigmaZZ(float sigmaZZ) { sigmaZZ_ = sigmaZZ;} /* operators */ @@ -186,7 +192,9 @@ namespace l1t //shower shape int showerLength_; + int coreShowerLength_; int firstLayer_; + int maxLayer_; float eMax_; float sigmaEtaEtaMax_; float sigmaPhiPhiMax_; @@ -194,6 +202,7 @@ namespace l1t float sigmaEtaEtaTot_; float sigmaPhiPhiTot_; float sigmaRRTot_; + float sigmaRRMean_; float sigmaZZ_; ClusterShapes shapes_; diff --git a/DataFormats/L1THGCal/src/classes_def.xml b/DataFormats/L1THGCal/src/classes_def.xml index 8933b357f2327..95ce8287cdf98 100644 --- a/DataFormats/L1THGCal/src/classes_def.xml +++ b/DataFormats/L1THGCal/src/classes_def.xml @@ -20,7 +20,8 @@ - + + @@ -30,7 +31,8 @@ - + + diff --git a/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h index 5b2129507505b..ef2f77a515d7f 100644 --- a/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h +++ b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h @@ -16,7 +16,10 @@ class HGCalShowerShape{ int firstLayer(const l1t::HGCalMulticluster& c3d) const; int lastLayer(const l1t::HGCalMulticluster& c3d) const; + int maxLayer(const l1t::HGCalMulticluster& c3d) const; int showerLength(const l1t::HGCalMulticluster& c3d) const {return lastLayer(c3d)-firstLayer(c3d)+1; }//in number of layers + // Maximum number of consecutive layers in the cluster + int coreShowerLength(const l1t::HGCalMulticluster& c3d) const; float eMax(const l1t::HGCalMulticluster& c3d) const; @@ -33,6 +36,7 @@ class HGCalShowerShape{ float sigmaRRTot(const l1t::HGCalMulticluster& c3d) const; float sigmaRRTot(const l1t::HGCalCluster& c2d) const; float sigmaRRMax(const l1t::HGCalMulticluster& c3d) const; + float sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius=5.) const; private: diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc index 0b95f3265f8f1..524da02bfa3bb 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc @@ -30,7 +30,9 @@ class HGCalTriggerNtupleHGCMulticlusters : public HGCalTriggerNtupleBase std::vector> cl3d_clusters_id_; // cluster shower shapes std::vector cl3d_showerlength_; + std::vector cl3d_coreshowerlength_; std::vector cl3d_firstlayer_; + std::vector cl3d_maxlayer_; std::vector cl3d_seetot_; std::vector cl3d_seemax_; std::vector cl3d_spptot_; @@ -38,6 +40,7 @@ class HGCalTriggerNtupleHGCMulticlusters : public HGCalTriggerNtupleBase std::vector cl3d_szz_; std::vector cl3d_srrtot_; std::vector cl3d_srrmax_; + std::vector cl3d_srrmean_; std::vector cl3d_emaxe_; }; @@ -66,7 +69,9 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("cl3d_clusters_n", &cl3d_clusters_n_); tree.Branch("cl3d_clusters_id", &cl3d_clusters_id_); tree.Branch("cl3d_showerlength", &cl3d_showerlength_); + tree.Branch("cl3d_coreshowerlength", &cl3d_coreshowerlength_); tree.Branch("cl3d_firstlayer", &cl3d_firstlayer_); + tree.Branch("cl3d_maxlayer", &cl3d_maxlayer_); tree.Branch("cl3d_seetot", &cl3d_seetot_); tree.Branch("cl3d_seemax", &cl3d_seemax_); tree.Branch("cl3d_spptot", &cl3d_spptot_); @@ -74,6 +79,7 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("cl3d_szz", &cl3d_szz_); tree.Branch("cl3d_srrtot", &cl3d_srrtot_); tree.Branch("cl3d_srrmax", &cl3d_srrmax_); + tree.Branch("cl3d_srrmean", &cl3d_srrmean_); tree.Branch("cl3d_emaxe", &cl3d_emaxe_); } @@ -104,7 +110,9 @@ fill(const edm::Event& e, const edm::EventSetup& es) cl3d_phi_.emplace_back(cl3d_itr->phi()); cl3d_clusters_n_.emplace_back(cl3d_itr->constituents().size()); cl3d_showerlength_.emplace_back(cl3d_itr->showerLength()); + cl3d_coreshowerlength_.emplace_back(cl3d_itr->coreShowerLength()); cl3d_firstlayer_.emplace_back(cl3d_itr->firstLayer()); + cl3d_maxlayer_.emplace_back(cl3d_itr->maxLayer()); cl3d_seetot_.emplace_back(cl3d_itr->sigmaEtaEtaTot()); cl3d_seemax_.emplace_back(cl3d_itr->sigmaEtaEtaMax()); cl3d_spptot_.emplace_back(cl3d_itr->sigmaPhiPhiTot()); @@ -112,6 +120,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) cl3d_szz_.emplace_back(cl3d_itr->sigmaZZ()); cl3d_srrtot_.emplace_back(cl3d_itr->sigmaRRTot()); cl3d_srrmax_.emplace_back(cl3d_itr->sigmaRRMax()); + cl3d_srrmean_.emplace_back(cl3d_itr->sigmaRRMean()); cl3d_emaxe_.emplace_back(cl3d_itr->eMax()/cl3d_itr->energy()); // Retrieve indices of trigger cells inside cluster @@ -136,7 +145,9 @@ clear() cl3d_clusters_n_.clear(); cl3d_clusters_id_.clear(); cl3d_showerlength_.clear(); + cl3d_coreshowerlength_.clear(); cl3d_firstlayer_.clear(); + cl3d_maxlayer_.clear(); cl3d_seetot_.clear(); cl3d_seemax_.clear(); cl3d_spptot_.clear(); @@ -144,6 +155,7 @@ clear() cl3d_szz_.clear(); cl3d_srrtot_.clear(); cl3d_srrmax_.clear(); + cl3d_srrmean_.clear(); cl3d_emaxe_.clear(); } diff --git a/L1Trigger/L1THGCal/src/be_algorithms/HGCalMulticlusteringImpl.cc b/L1Trigger/L1THGCal/src/be_algorithms/HGCalMulticlusteringImpl.cc index aeecab5c1dba9..0b098b3f833d9 100644 --- a/L1Trigger/L1THGCal/src/be_algorithms/HGCalMulticlusteringImpl.cc +++ b/L1Trigger/L1THGCal/src/be_algorithms/HGCalMulticlusteringImpl.cc @@ -145,7 +145,9 @@ void HGCalMulticlusteringImpl::clusterizeDR( const edm::PtrVector& clustersPtrs = c3d.constituents(); + std::unordered_map layers_pt; + float max_pt = 0.; + int max_layer = 0; + for(const auto& cluster_ptr : clustersPtrs){ + int layer = HGC_layer(cluster_ptr->subdetId(),cluster_ptr->layer()); + auto itr_insert = layers_pt.emplace(layer, 0.); + itr_insert.first->second += cluster_ptr->pt(); + if(itr_insert.first->second>max_pt){ + max_pt = itr_insert.first->second; + max_layer = layer; + } + } + return max_layer; +} + int HGCalShowerShape::lastLayer(const l1t::HGCalMulticluster& c3d) const { @@ -126,8 +144,25 @@ int HGCalShowerShape::lastLayer(const l1t::HGCalMulticluster& c3d) const { } - - +int HGCalShowerShape::coreShowerLength(const l1t::HGCalMulticluster& c3d) const +{ + const edm::PtrVector& clustersPtrs = c3d.constituents(); + std::vector layers(kLayersEE_+kLayersFH_+kLayersBH_); + for(const auto& cluster_ptr : clustersPtrs) + { + int layer = HGC_layer(cluster_ptr->subdetId(), cluster_ptr->layer()); + layers[layer-1] = true; + } + int length = 0; + int maxlength = 0; + for(bool layer : layers) + { + if(layer) length++; + else length = 0; + if(length>maxlength) maxlength = length; + } + return maxlength; +} float HGCalShowerShape::sigmaEtaEtaTot(const l1t::HGCalMulticluster& c3d) const { @@ -343,6 +378,57 @@ float HGCalShowerShape::sigmaRRMax(const l1t::HGCalMulticluster& c3d) const { } +float HGCalShowerShape::sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius) const { + + const edm::PtrVector& clustersPtrs = c3d.constituents(); + // group trigger cells by layer + std::unordered_map> > layers_tcs; + for(const auto& clu : clustersPtrs){ + int layer = HGC_layer(clu->subdetId(),clu->layer()); + const edm::PtrVector& triggerCells = clu->constituents(); + for(const auto& tc : triggerCells){ + layers_tcs[layer].emplace_back(tc); + } + } + + // Select trigger cells within X cm of the max TC in the layer + std::unordered_map > > tc_layers_energy_r; + for(const auto& layer_tcs : layers_tcs){ + int layer = layer_tcs.first; + edm::Ptr max_tc = layer_tcs.second.front(); + for(const auto& tc : layer_tcs.second){ + if(tc->energy()>max_tc->energy()) max_tc = tc; + } + for(const auto& tc : layer_tcs.second){ + double dx = tc->position().x() - max_tc->position().x(); + double dy = tc->position().y() - max_tc->position().y(); + double distance_to_max = std::sqrt(dx*dx+dy*dy); + if(distance_to_maxposition().x()*tc->position().x() + tc->position().y()*tc->position().y())/std::abs(tc->position().z()); + tc_layers_energy_r[layer].emplace_back( std::make_pair(tc->energy(), r)); + } + } + } + + // Compute srr layer by layer + std::vector> layers_energy_srr2; + for(const auto& layer_energy_r : tc_layers_energy_r){ + const auto& energy_r = layer_energy_r.second; + float r_mean_layer = meanX(energy_r); + float srr = sigmaXX(energy_r,r_mean_layer); + double energy_sum = 0.; + for(const auto& energy_r : energy_r){ + energy_sum += energy_r.first; + } + layers_energy_srr2.emplace_back(std::make_pair(energy_sum, srr*srr)); + } + // Combine all layer srr + float srr2_mean = meanX(layers_energy_srr2); + return std::sqrt(srr2_mean); +} + + + float HGCalShowerShape::eMax(const l1t::HGCalMulticluster& c3d) const { std::unordered_map layer_energy; From 454ac835ef8b68d6668ba89703e5100c66849281 Mon Sep 17 00:00:00 2001 From: Jean-Baptiste Sauvan Date: Tue, 31 Oct 2017 10:04:52 +0100 Subject: [PATCH 7/7] Fix ntuple lines --- .../plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc | 4 ++-- .../plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc index e8e5201586b77..2a0282ceed5f5 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc @@ -25,8 +25,8 @@ class HGCalTriggerNtupleHGCClusters : public HGCalTriggerNtupleBase edm::EDGetToken clusters_token_, multiclusters_token_; int cl_n_ ; - std::vector cl_mipPt_; std::vector cl_id_; + std::vector cl_mipPt_; std::vector cl_pt_; std::vector cl_energy_; std::vector cl_eta_; @@ -59,8 +59,8 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& multiclusters_token_ = collector.consumes(conf.getParameter("Multiclusters")); tree.Branch("cl_n", &cl_n_, "cl_n/I"); - tree.Branch("cl_mipPt", &cl_mipPt_); tree.Branch("cl_id", &cl_id_); + tree.Branch("cl_mipPt", &cl_mipPt_); tree.Branch("cl_pt", &cl_pt_); tree.Branch("cl_energy", &cl_energy_); tree.Branch("cl_eta", &cl_eta_); diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc index 1804f4a303038..e560a31d33405 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc @@ -44,11 +44,11 @@ class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase std::vector tc_cell_; std::vector tc_data_; std::vector tc_mipPt_; + std::vector tc_pt_; std::vector tc_energy_; std::vector tc_simenergy_; std::vector tc_eta_; std::vector tc_phi_; - std::vector tc_pt_; std::vector tc_x_; std::vector tc_y_; std::vector tc_z_; @@ -93,12 +93,12 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("tc_wafertype", &tc_wafertype_); tree.Branch("tc_cell", &tc_cell_); tree.Branch("tc_data", &tc_data_); + tree.Branch("tc_pt", &tc_pt_); tree.Branch("tc_mipPt", &tc_mipPt_); tree.Branch("tc_energy", &tc_energy_); if(fill_simenergy_) tree.Branch("tc_simenergy", &tc_simenergy_); tree.Branch("tc_eta", &tc_eta_); tree.Branch("tc_phi", &tc_phi_); - tree.Branch("tc_pt", &tc_pt_); tree.Branch("tc_x", &tc_x_); tree.Branch("tc_y", &tc_y_); tree.Branch("tc_z", &tc_z_); @@ -174,10 +174,10 @@ fill(const edm::Event& e, const edm::EventSetup& es) tc_data_.emplace_back(tc_itr->hwPt()); tc_mipPt_.emplace_back(tc_itr->mipPt()); // physical values + tc_pt_.emplace_back(tc_itr->pt()); tc_energy_.emplace_back(tc_itr->energy()); tc_eta_.emplace_back(tc_itr->eta()); tc_phi_.emplace_back(tc_itr->phi()); - tc_pt_.emplace_back(tc_itr->pt()); tc_x_.emplace_back(tc_itr->position().x()); tc_y_.emplace_back(tc_itr->position().y()); tc_z_.emplace_back(tc_itr->position().z()); @@ -292,11 +292,11 @@ clear() tc_cell_.clear(); tc_data_.clear(); tc_mipPt_.clear(); + tc_pt_.clear(); tc_energy_.clear(); tc_simenergy_.clear(); tc_eta_.clear(); tc_phi_.clear(); - tc_pt_.clear(); tc_x_.clear(); tc_y_.clear(); tc_z_.clear();