From fd56c214360e5a8af00e453b34599ff3e64c9c99 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 13 Sep 2016 10:53:18 +0200 Subject: [PATCH 01/21] Adding Consumes Collector and Events to the BackEnd FWK --- .../HGCalTriggerBackendAlgorithmBase.h | 24 +++++++++++++++---- .../interface/HGCalTriggerBackendProcessor.h | 6 +++-- .../fe_codecs/HGCalBestChoiceCodecImpl.h | 1 + .../plugins/HGCalTriggerDigiFEReproducer.cc | 4 ++-- .../plugins/HGCalTriggerDigiProducer.cc | 4 ++-- .../be_algorithms/FullModuleSumAlgo.cc | 12 ++++++---- .../be_algorithms/RandomClusterAlgo.cc | 12 ++++++---- .../be_algorithms/SingleCellClusterAlgo.cc | 12 ++++++---- .../src/HGCalTriggerBackendProcessor.cc | 9 +++---- 9 files changed, 58 insertions(+), 26 deletions(-) diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h index 42661e672b34d..2137a673f5425 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h @@ -29,9 +29,17 @@ class HGCalTriggerBackendAlgorithmBase { public: + /* HGCalTriggerBackendAlgorithmBase(const edm::ParameterSet& conf) : name_(conf.getParameter("AlgorithmName")) {} + */ + + // Allow HGCalTriggerBackend to be passed a consume collector + HGCalTriggerBackendAlgorithmBase(const edm::ParameterSet& conf, edm::ConsumesCollector &cc) : + name_(conf.getParameter("AlgorithmName")) + {} + virtual ~HGCalTriggerBackendAlgorithmBase() {} const std::string& name() const { return name_; } @@ -40,7 +48,9 @@ class HGCalTriggerBackendAlgorithmBase { virtual void setProduces(edm::EDProducer& prod) const = 0; virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom) = 0; + const std::unique_ptr& geom, + const edm::Event &e + ) = 0; virtual void putInEvent(edm::Event& evt) = 0; @@ -57,9 +67,15 @@ namespace HGCalTriggerBackend { template class Algorithm : public HGCalTriggerBackendAlgorithmBase { public: + /* Algorithm(const edm::ParameterSet& conf) : - HGCalTriggerBackendAlgorithmBase(conf), - codec_(conf.getParameterSet("FECodec")){ } + HGCalTriggerBackendAlgorithmBase(conf), + codec_(conf.getParameterSet("FECodec")){ } + */ + + Algorithm(const edm::ParameterSet& conf, edm::ConsumesCollector &cc ) : + HGCalTriggerBackendAlgorithmBase(conf, cc), + codec_(conf.getParameterSet("FECodec")){ } protected: FECODEC codec_; @@ -67,6 +83,6 @@ namespace HGCalTriggerBackend { } #include "FWCore/PluginManager/interface/PluginFactory.h" -typedef edmplugin::PluginFactory< HGCalTriggerBackendAlgorithmBase* (const edm::ParameterSet&) > HGCalTriggerBackendAlgorithmFactory; +typedef edmplugin::PluginFactory< HGCalTriggerBackendAlgorithmBase* (const edm::ParameterSet&,edm::ConsumesCollector & ) > HGCalTriggerBackendAlgorithmFactory; #endif diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h index e93f8067e50d9..7a9c5d2816d4e 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h @@ -31,12 +31,14 @@ class HGCalTriggerBackendProcessor { public: typedef std::unique_ptr algo_ptr; - HGCalTriggerBackendProcessor(const edm::ParameterSet& conf); + HGCalTriggerBackendProcessor(const edm::ParameterSet& conf, edm::ConsumesCollector&&cc); void setProduces(edm::EDProducer& prod) const; void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom); + const std::unique_ptr& geom, + const edm::Event&e + ); void putInEvent(edm::Event& evt); diff --git a/L1Trigger/L1THGCal/interface/fe_codecs/HGCalBestChoiceCodecImpl.h b/L1Trigger/L1THGCal/interface/fe_codecs/HGCalBestChoiceCodecImpl.h index 21b9206e7e4f5..e2cf1c4f93bea 100644 --- a/L1Trigger/L1THGCal/interface/fe_codecs/HGCalBestChoiceCodecImpl.h +++ b/L1Trigger/L1THGCal/interface/fe_codecs/HGCalBestChoiceCodecImpl.h @@ -21,6 +21,7 @@ struct HGCalBestChoiceDataPayload { payload.fill(0); } + }; diff --git a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc index 7932985197839..4362f30e0f4ea 100644 --- a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc +++ b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc @@ -40,7 +40,7 @@ DEFINE_FWK_MODULE(HGCalTriggerDigiFEReproducer); /*****************************************************************/ HGCalTriggerDigiFEReproducer::HGCalTriggerDigiFEReproducer(const edm::ParameterSet& conf): inputdigi_(consumes(conf.getParameter("feDigis"))), - backEndProcessor_(conf.getParameterSet("BEConfiguration")) + backEndProcessor_(conf.getParameterSet("BEConfiguration"), consumesCollector()) /*****************************************************************/ { //setup geometry configuration @@ -113,7 +113,7 @@ void HGCalTriggerDigiFEReproducer::produce(edm::Event& e, const edm::EventSetup& auto fe_digis_coll = *fe_digis_handle; //now we run the emulation of the back-end processor - backEndProcessor_.run(fe_digis_coll,triggerGeometry_); + backEndProcessor_.run(fe_digis_coll,triggerGeometry_,e); backEndProcessor_.putInEvent(e); backEndProcessor_.reset(); } diff --git a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc index 821bac2ca3a15..f5f3311be6130 100644 --- a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc +++ b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc @@ -40,7 +40,7 @@ HGCalTriggerDigiProducer(const edm::ParameterSet& conf): inputee_(consumes(conf.getParameter("eeDigis"))), inputfh_(consumes(conf.getParameter("fhDigis"))), //inputbh_(consumes(conf.getParameter("bhDigis"))), - backEndProcessor_(conf.getParameterSet("BEConfiguration")) { + backEndProcessor_(conf.getParameterSet("BEConfiguration"),consumesCollector()) { //setup geometry configuration const edm::ParameterSet& geometryConfig = @@ -147,7 +147,7 @@ void HGCalTriggerDigiProducer::produce(edm::Event& e, const edm::EventSetup& es) auto fe_digis_coll = *fe_digis_handle; //now we run the emulation of the back-end processor - backEndProcessor_.run(fe_digis_coll,triggerGeometry_); + backEndProcessor_.run(fe_digis_coll,triggerGeometry_,e); backEndProcessor_.putInEvent(e); backEndProcessor_.reset(); } diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc index d627217bf3e64..70b80309e64f6 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc @@ -10,8 +10,8 @@ class FullModuleSumAlgo : public Algorithm { public: - FullModuleSumAlgo(const edm::ParameterSet& conf): - Algorithm(conf), + FullModuleSumAlgo(const edm::ParameterSet& conf, edm::ConsumesCollector& cc): + Algorithm(conf,cc), cluster_product_( new l1t::HGCalClusterBxCollection ){} virtual void setProduces(edm::EDProducer& prod) const override final @@ -20,7 +20,9 @@ class FullModuleSumAlgo : public Algorithm } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom) override final; + const std::unique_ptr& geom, + const edm::Event&evt + ) override final; virtual void putInEvent(edm::Event& evt) override final { @@ -39,7 +41,9 @@ class FullModuleSumAlgo : public Algorithm /*****************************************************************/ void FullModuleSumAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom) + const std::unique_ptr& geom, + const edm::Event&evt + ) /*****************************************************************/ { for( const auto& digi : coll ) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc index dcceb29bfb208..cb026c03b00cb 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc @@ -8,8 +8,8 @@ using namespace HGCalTriggerBackend; class RandomClusterAlgo : public Algorithm { public: - RandomClusterAlgo(const edm::ParameterSet& conf): - Algorithm(conf), + RandomClusterAlgo(const edm::ParameterSet& conf,edm::ConsumesCollector &cc): + Algorithm(conf,cc), cluster_product_( new l1t::HGCalClusterBxCollection ){ } @@ -18,7 +18,9 @@ class RandomClusterAlgo : public Algorithm { } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom) override final; + const std::unique_ptr& geom, + const edm::Event&evt + ) override final; virtual void putInEvent(edm::Event& evt) override final { evt.put(std::move(cluster_product_),name()); @@ -34,7 +36,9 @@ class RandomClusterAlgo : public Algorithm { }; void RandomClusterAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom) { + const std::unique_ptr& geom, + const edm::Event&evt + ) { for( const auto& digi : coll ) { HGCal64BitRandomCodec::data_type my_data; digi.decode(codec_,my_data); diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc index c0b4042c52a8c..a3c1023a25d71 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc @@ -10,8 +10,8 @@ class SingleCellClusterAlgo : public Algorithm { public: - SingleCellClusterAlgo(const edm::ParameterSet& conf): - Algorithm(conf), + SingleCellClusterAlgo(const edm::ParameterSet& conf,edm::ConsumesCollector&cc): + Algorithm(conf,cc), cluster_product_( new l1t::HGCalClusterBxCollection ){} virtual void setProduces(edm::EDProducer& prod) const override final @@ -20,7 +20,9 @@ class SingleCellClusterAlgo : public Algorithm } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom) override final; + const std::unique_ptr& geom, + const edm::Event&evt + ) override final; virtual void putInEvent(edm::Event& evt) override final { @@ -39,7 +41,9 @@ class SingleCellClusterAlgo : public Algorithm /*****************************************************************/ void SingleCellClusterAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom) + const std::unique_ptr& geom, + const edm::Event&evt + ) /*****************************************************************/ { for( const auto& digi : coll ) diff --git a/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc b/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc index edfdcf7b6865a..afe0feb10d893 100644 --- a/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc +++ b/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc @@ -1,14 +1,14 @@ #include "L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h" HGCalTriggerBackendProcessor:: -HGCalTriggerBackendProcessor(const edm::ParameterSet& conf) { +HGCalTriggerBackendProcessor(const edm::ParameterSet& conf, edm::ConsumesCollector&& cc) { const std::vector be_confs = conf.getParameterSetVector("algorithms"); for( const auto& algo_cfg : be_confs ) { const std::string& algo_name = algo_cfg.getParameter("AlgorithmName"); HGCalTriggerBackendAlgorithmBase* algo = - HGCalTriggerBackendAlgorithmFactory::get()->create(algo_name,algo_cfg); + HGCalTriggerBackendAlgorithmFactory::get()->create(algo_name,algo_cfg,cc); algorithms_.emplace_back(algo); } } @@ -21,9 +21,10 @@ void HGCalTriggerBackendProcessor::setProduces(edm::EDProducer& prod) const { void HGCalTriggerBackendProcessor:: run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom) { + const std::unique_ptr& geom, + const edm::Event &e) { for( auto& algo : algorithms_ ) { - algo->run(coll,geom); + algo->run(coll,geom,e); } } From ee51c97c06f9014e668f0e03c2a3aeb4952fb2e6 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 13 Sep 2016 10:53:37 +0200 Subject: [PATCH 02/21] incomplete version of Simclustering --- .../plugins/be_algorithms/HGCalSimCluster.cc | 135 ++++++++++++++++++ 1 file changed, 135 insertions(+) create mode 100644 L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc new file mode 100644 index 0000000000000..9587d515a0eca --- /dev/null +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -0,0 +1,135 @@ +// HGCal Trigger +#include "L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h" +#include "L1Trigger/L1THGCal/interface/fe_codecs/HGCalBestChoiceCodec.h" +#include "DataFormats/ForwardDetId/interface/HGCTriggerDetId.h" + +// HGCalClusters and detId +#include "DataFormats/L1THGCal/interface/HGCalCluster.h" +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" + +// PF Cluster definition +#include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h" + +// Consumes +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/Framework/interface/EDConsumerBase.h" + + + +/* Original Author: Andrea Carlo Marini + * Original Dat: 23 Aug 2016 + * + * This backend algorithm is supposed to use the sim cluster information to handle an + * optimal clustering algorithm, benchmark the performances of the enconding ... + */ + + +namespace HGCalTriggerBackend{ + template + class HGCalTriggerSimCluster : public Algorithm + { + private: + std::unique_ptr cluster_product; + // token + edm::EDGetTokenT< reco::PFClusterCollection > sim_token; + // handle + edm::Handle< reco::PFClusterCollection > sim_handle; + + protected: + using Algorithm::codec_; + + public: + // Constructor + // + using Algorithm::Algorithm; + using Algorithm::name; + using Algorithm::run; + using Algorithm::putInEvent; + using Algorithm::setProduces; + using Algorithm::reset; + + //Consumes tokens + HGCalTriggerSimCluster(const edm::ParameterSet& conf,edm::ConsumesCollector&cc) : Algorithm(conf,cc) { + // I need to consumes the PF Cluster Collection with the sim clustering, TODO: make it configurable (?) + sim_token = cc.consumes< reco::PFClusterCollection >(edm::InputTag("particleFlowClusterHGCal")); + } + + // setProduces + virtual void setProduces(edm::EDProducer& prod) const override final + { + prod.produces(name()); + } + + // putInEvent + virtual void putInEvent(edm::Event& evt) override final + { + evt.put(std::move(cluster_product),name()); + } + + //reset + virtual void reset() override final + { + cluster_product.reset( new l1t::HGCalClusterBxCollection ); + } + + // run, actual algorithm + virtual void run( const l1t::HGCFETriggerDigiCollection & coll, + const std::unique_ptr&geom, + const edm::Event&evt + ) + { + //1. construct a cluster container that hosts the cluster per truth-particle + std::unordered_map > cluster_container;// PID-> bx,cluster + evt.getByToken(sim_token,sim_handle); + //2. run on the digits, + for( const auto& digi : coll ) + { + DATA data; + const HGCTriggerDetId& tcellId = digi.getDetId(); + digi.decode(codec_,data); + //2.A get the trigger-cell information energy/id + //uint32_t digiEnergy = data.payload; i + auto digiEnergy= data.energy(); // if it is implemented an energy() method, etherwise it will not compile + // XXX this depends on the DATA type + //2.B get the HGCAL-base-cell associated to it / geometry + //const auto& tc=geom->triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& + //for(const auto& cell : tc.components() ) // HGcell -- unsigned + for(const auto& cell : geom->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned + { + HGCalDetId cellId(cell); + //2.C get the particleId and energy fractions + //2.D add to the corresponding cluster + } + } + + //3. Push the clusters in the cluster_product + uint32_t clusterEnergyHw=0; + uint32_t clusterEtaHw = 0 ;//tcellId(); + //const GlobalPoint& tcellPosition = geom->getTriggerCellPosition( tcellId()); + + // construct it from *both* physical and integer values + l1t::HGCalCluster cluster( reco::LeafCandidate::LorentzVector(), + clusterEnergyHw, clusterEtaHw, 0); + + cluster_product->push_back(0,cluster); // bx,cluster + + } + + + }; + +}// namespace + + +// define plugins, template needs to be spelled out here, in order to allow the compiler to compile, and the factory to be populated +typedef HGCalTriggerBackend::HGCalTriggerSimCluster HGCalTriggerSimClusterBestChoice; +DEFINE_EDM_PLUGIN(HGCalTriggerBackendAlgorithmFactory, HGCalTriggerSimClusterBestChoice,"HGCalTriggerSimClusterBestChoice"); +//DEFINE_EDM_PLUGIN(HGCalTriggerBackendAlgorithmFactory, HGCalTriggerSimCluster,"HGCalTriggerSimCluster"); + +// Local Variables: +// mode:c++ +// indent-tabs-mode:nil +// tab-width:4 +// c-basic-offset:4 +// End: +// vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 From e8717dabaea5f5333ce59632d38e92d58f3e9131 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Thu, 22 Sep 2016 10:29:39 +0200 Subject: [PATCH 03/21] needs still some work --- .../plugins/be_algorithms/HGCalSimCluster.cc | 62 +++++++++++++++---- 1 file changed, 51 insertions(+), 11 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 9587d515a0eca..296af04701a7e 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -14,6 +14,8 @@ #include "FWCore/Framework/interface/ConsumesCollector.h" #include "FWCore/Framework/interface/EDConsumerBase.h" +// Print out something before crashing or throwing exceptions +#include /* Original Author: Andrea Carlo Marini @@ -24,7 +26,9 @@ */ + namespace HGCalTriggerBackend{ + template class HGCalTriggerSimCluster : public Algorithm { @@ -81,24 +85,60 @@ namespace HGCalTriggerBackend{ //1. construct a cluster container that hosts the cluster per truth-particle std::unordered_map > cluster_container;// PID-> bx,cluster evt.getByToken(sim_token,sim_handle); + + if (not sim_handle.isValid()) { std::cout<<"[HGCalTriggerSimCluster]::[run]::[ERROR] PFCluster collection for HGC sim clustering not available"< [ (pid, fraction), ... + std::unordered_map > > simclusters; + for (auto& cluster : *sim_handle) + { + auto& pid= cluster.particleId(); // not pdgId + auto& hf = cluster.hits_and_fractions(); + for (const auto & p : hf ) + { + simclusters[p.first].push_back( std::pair( pid,p.second) ) ; + } + } + + //2. run on the digits, for( const auto& digi : coll ) { DATA data; - const HGCTriggerDetId& tcellId = digi.getDetId(); digi.decode(codec_,data); //2.A get the trigger-cell information energy/id - //uint32_t digiEnergy = data.payload; i - auto digiEnergy= data.energy(); // if it is implemented an energy() method, etherwise it will not compile - // XXX this depends on the DATA type - //2.B get the HGCAL-base-cell associated to it / geometry - //const auto& tc=geom->triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& - //for(const auto& cell : tc.components() ) // HGcell -- unsigned - for(const auto& cell : geom->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned + const HGCTriggerDetId& moduleId = digi.getDetId(); // this is a module Det Id + + const auto& trcells=geom->getOrderedTriggerCellsFromModule( moduleId() ); + + // check + if ( trcells.size() != data.payload.size() ) + { + std::cout<<"[HGCalTriggerSimCluster]::[run]::[ERROR] Mapping outside of assumptions."<triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& + //for(const auto& cell : tc.components() ) // HGcell -- unsigned + for(const auto& cell : geom->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned + { + HGCalDetId cellId(cell); + //2.C get the particleId and energy fractions + //2.D add to the corresponding cluster + } } } From d9e36bd6ff6dc618748819b2a9c5be4e22f2f641 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Mon, 26 Sep 2016 13:42:41 +0200 Subject: [PATCH 04/21] compilation checkpoint --- .../plugins/be_algorithms/HGCalSimCluster.cc | 49 ++++++++++--------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 296af04701a7e..968b0a220b8a6 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -9,6 +9,8 @@ // PF Cluster definition #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h" +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" + // Consumes #include "FWCore/Framework/interface/ConsumesCollector.h" @@ -35,9 +37,11 @@ namespace HGCalTriggerBackend{ private: std::unique_ptr cluster_product; // token - edm::EDGetTokenT< reco::PFClusterCollection > sim_token; + //edm::EDGetTokenT< reco::PFClusterCollection > sim_token; + edm::EDGetTokenT< std::vector > sim_token; // handle - edm::Handle< reco::PFClusterCollection > sim_handle; + //edm::Handle< reco::PFClusterCollection > sim_handle; + edm::Handle< std::vector > sim_handle; protected: using Algorithm::codec_; @@ -55,7 +59,7 @@ namespace HGCalTriggerBackend{ //Consumes tokens HGCalTriggerSimCluster(const edm::ParameterSet& conf,edm::ConsumesCollector&cc) : Algorithm(conf,cc) { // I need to consumes the PF Cluster Collection with the sim clustering, TODO: make it configurable (?) - sim_token = cc.consumes< reco::PFClusterCollection >(edm::InputTag("particleFlowClusterHGCal")); + sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("particleFlowClusterHGCal")); } // setProduces @@ -93,8 +97,8 @@ namespace HGCalTriggerBackend{ std::unordered_map > > simclusters; for (auto& cluster : *sim_handle) { - auto& pid= cluster.particleId(); // not pdgId - auto& hf = cluster.hits_and_fractions(); + auto pid= cluster.particleId(); // not pdgId + const auto& hf = cluster.hits_and_fractions(); for (const auto & p : hf ) { simclusters[p.first].push_back( std::pair( pid,p.second) ) ; @@ -122,24 +126,25 @@ namespace HGCalTriggerBackend{ // there is a loss of generality here, due to the restriction imposed by the data formats // it will work if inside a module there is a data.payload with an ordered list of all the energies // one may think to add on top of it a wrapper if this stop to be the case for some of the data classes - for( const auto valIt=data.payload.begin(), const auto tcIt = trcells.begin(); - valIt != data.payload.end() and tcIt != trcells.end(); - valIt++, tcIt++ - ) - { - const HGCTriggerDetId tcellId(*tcIt); - //uint32_t digiEnergy = data.payload; i - auto digiEnergy=*valIt; // if it is implemented an energy() method, etherwise it will not compile - //2.B get the HGCAL-base-cell associated to it / geometry - //const auto& tc=geom->triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& - //for(const auto& cell : tc.components() ) // HGcell -- unsigned - for(const auto& cell : geom->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned + { //scope + auto valIt=data.payload.begin(); + auto tcIt = trcells.begin(); + for( ; valIt != data.payload.end() and tcIt != trcells.end(); valIt++, tcIt++) { - HGCalDetId cellId(cell); - //2.C get the particleId and energy fractions - //2.D add to the corresponding cluster - } - } + const HGCTriggerDetId tcellId(*tcIt); + //uint32_t digiEnergy = data.payload; i + auto digiEnergy=*valIt; // if it is implemented an energy() method, etherwise it will not compile + //2.B get the HGCAL-base-cell associated to it / geometry + //const auto& tc=geom->triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& + //for(const auto& cell : tc.components() ) // HGcell -- unsigned + for(const auto& cell : geom->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned + { + HGCalDetId cellId(cell); + //2.C get the particleId and energy fractions + //2.D add to the corresponding cluster + } + }// end of for loop + } //end of for-scope } //3. Push the clusters in the cluster_product From c879713c3e502b5beb8cce8f2e09c8a843e52ec6 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 27 Sep 2016 14:27:34 +0200 Subject: [PATCH 05/21] commit --- L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 968b0a220b8a6..cf8ff4c9383b1 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -37,10 +37,8 @@ namespace HGCalTriggerBackend{ private: std::unique_ptr cluster_product; // token - //edm::EDGetTokenT< reco::PFClusterCollection > sim_token; edm::EDGetTokenT< std::vector > sim_token; // handle - //edm::Handle< reco::PFClusterCollection > sim_handle; edm::Handle< std::vector > sim_handle; protected: @@ -59,7 +57,9 @@ namespace HGCalTriggerBackend{ //Consumes tokens HGCalTriggerSimCluster(const edm::ParameterSet& conf,edm::ConsumesCollector&cc) : Algorithm(conf,cc) { // I need to consumes the PF Cluster Collection with the sim clustering, TODO: make it configurable (?) - sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("particleFlowClusterHGCal")); + // vector "mix" "MergedCaloTruth" "HLT" + // pf clusters cannot be safely cast to SimCluster + sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("MergedCaloTruth")); } // setProduces From b9344b37c5d1d8d1c0c7e9ac12f723f62b3cdbaa Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Wed, 5 Oct 2016 11:43:56 +0200 Subject: [PATCH 06/21] mostly implemented. miss eta/phi info --- .../plugins/be_algorithms/HGCalSimCluster.cc | 55 ++++++++++++++++--- 1 file changed, 48 insertions(+), 7 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index cf8ff4c9383b1..7f42ed20ad31f 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -16,6 +16,9 @@ #include "FWCore/Framework/interface/ConsumesCollector.h" #include "FWCore/Framework/interface/EDConsumerBase.h" +// +#include "DataFormats/Math/interface/LorentzVector.h" + // Print out something before crashing or throwing exceptions #include @@ -40,6 +43,33 @@ namespace HGCalTriggerBackend{ edm::EDGetTokenT< std::vector > sim_token; // handle edm::Handle< std::vector > sim_handle; + // add to cluster + void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float energy,float eta, float phi) + { + auto iterator = cluster_container.find (pid); + if (iterator == cluster_container.end()) + { + //create an empty cluster + cluster_container[pid] = std::pair(0,l1t::HGCalCluster()); + iterator = cluster_container.find (pid); + iterator -> second . second . setPdgId(pdgid); + } + // p4 += p4' + math::PtEtaPhiMLorentzVectorD p4; + p4.SetPt ( iterator -> second . second . pt() ) ; + p4.SetEta( iterator -> second . second . eta() ) ; + p4.SetPhi( iterator -> second . second . phi() ) ; + p4.SetM ( iterator -> second . second . mass() ) ; + math::PtEtaPhiMLorentzVectorD pp4; + float t = std::exp (- eta); + pp4.SetPt ( energy * (1-t*t)/(1+t*t) ) ; + pp4.SetEta( eta ) ; + pp4.SetPhi( phi ) ; + pp4.SetM ( 0 ) ; + p4 += pp4; + iterator -> second . second . setP4(p4); + return ; + } protected: using Algorithm::codec_; @@ -141,22 +171,33 @@ namespace HGCalTriggerBackend{ { HGCalDetId cellId(cell); //2.C get the particleId and energy fractions - //2.D add to the corresponding cluster + const auto& particles = simclusters[cellId]; // vector pid fractions + for ( const auto& p: particles ) + { + const auto & pid= p.first; + const auto & fraction=p.second; + auto energy = fraction*digiEnergy; + //2.D add to the corresponding cluster + //void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float & energy,float & eta, float &phi) + //addToCluster(cluster_container, pid, 0 energy,ETA/PHI? ) ; + addToCluster(cluster_container, pid, 0,energy,0.,0. ) ; // how do I get eta, phi w/o the hgcal geometry? + } } }// end of for loop } //end of for-scope } //3. Push the clusters in the cluster_product - uint32_t clusterEnergyHw=0; - uint32_t clusterEtaHw = 0 ;//tcellId(); + //uint32_t clusterEnergyHw=0; + //uint32_t clusterEtaHw = 0 ;//tcellId(); //const GlobalPoint& tcellPosition = geom->getTriggerCellPosition( tcellId()); // construct it from *both* physical and integer values - l1t::HGCalCluster cluster( reco::LeafCandidate::LorentzVector(), - clusterEnergyHw, clusterEtaHw, 0); - - cluster_product->push_back(0,cluster); // bx,cluster + //l1t::HGCalCluster cluster( reco::LeafCandidate::LorentzVector(), + // clusterEnergyHw, clusterEtaHw, 0); + // + for (auto& p : cluster_container) + cluster_product->push_back(p.second.first,p.second.second); // bx,cluster } From 2c3d92613a96832cbafcc73dd163f2fc28920e89 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Mon, 10 Oct 2016 16:43:02 +0200 Subject: [PATCH 07/21] sim clustering --- .../plugins/be_algorithms/HGCalSimCluster.cc | 15 +++++++++++++++ L1Trigger/L1THGCal/test/testHGCalL1T_cfg.py | 3 ++- 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 7f42ed20ad31f..2fe4bc2e517cf 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -30,6 +30,8 @@ * optimal clustering algorithm, benchmark the performances of the enconding ... */ +#define DEBUG + namespace HGCalTriggerBackend{ @@ -116,6 +118,9 @@ namespace HGCalTriggerBackend{ const edm::Event&evt ) { +#ifdef DEBUG + cout<<"[HGCalTriggerSimCluster]::[run] Start"< > cluster_container;// PID-> bx,cluster evt.getByToken(sim_token,sim_handle); @@ -123,6 +128,9 @@ namespace HGCalTriggerBackend{ if (not sim_handle.isValid()) { std::cout<<"[HGCalTriggerSimCluster]::[run]::[ERROR] PFCluster collection for HGC sim clustering not available"< [ (pid, fraction), ... std::unordered_map > > simclusters; for (auto& cluster : *sim_handle) @@ -135,6 +143,9 @@ namespace HGCalTriggerBackend{ } } +#ifdef DEBUG + cout<<"[HGCalTriggerSimCluster]::[run] Run on digis"< Date: Tue, 18 Oct 2016 18:57:28 +0200 Subject: [PATCH 08/21] testing and adapting to pre12 --- .../HGCalTriggerBackendAlgorithmBase.h | 2 +- .../interface/HGCalTriggerBackendProcessor.h | 3 +- .../plugins/HGCalTriggerDigiFEReproducer.cc | 12 ++--- .../plugins/HGCalTriggerDigiProducer.cc | 11 ++--- .../be_algorithms/FullModuleSumAlgo.cc | 6 +-- .../plugins/be_algorithms/HGCalSimCluster.cc | 44 +++++++++++++++---- .../be_algorithms/RandomClusterAlgo.cc | 4 +- .../be_algorithms/SingleCellClusterAlgo.cc | 8 ++-- .../src/HGCalTriggerBackendProcessor.cc | 3 +- 9 files changed, 62 insertions(+), 31 deletions(-) diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h index 0f8ded5b2790d..c50e10aede028 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h @@ -51,7 +51,7 @@ class HGCalTriggerBackendAlgorithmBase { virtual void setProduces(edm::EDProducer& prod) const = 0; virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom, + const edm::ESHandle& geom, const edm::Event &e ) = 0; diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h index 8be39318950c3..eff1dca3b73fe 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h @@ -38,7 +38,8 @@ class HGCalTriggerBackendProcessor { void setProduces(edm::EDProducer& prod) const; void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom, + const edm::ESHandle& geom, + //const edm::ESHandle & geom, const edm::Event&e ); diff --git a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc index 9a1d7d19b9360..7db51345e218d 100644 --- a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc +++ b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc @@ -39,7 +39,7 @@ DEFINE_FWK_MODULE(HGCalTriggerDigiFEReproducer); /*****************************************************************/ HGCalTriggerDigiFEReproducer::HGCalTriggerDigiFEReproducer(const edm::ParameterSet& conf): inputdigi_(consumes(conf.getParameter("feDigis"))), - backEndProcessor_(conf.getParameterSet("BEConfiguration"), consumesCollector()) + backEndProcessor_(new HGCalTriggerBackendProcessor(conf.getParameterSet("BEConfiguration"), consumesCollector())) /*****************************************************************/ { //setup FE codec @@ -50,8 +50,8 @@ HGCalTriggerDigiFEReproducer::HGCalTriggerDigiFEReproducer(const edm::ParameterS codec_->unSetDataPayload(); produces(); - //setup BE processor - backEndProcessor_ = std::make_unique(conf.getParameterSet("BEConfiguration")); + //setup BE processor it's in the constructor + //backEndProcessor_ = std::make_unique(conf.getParameterSet("BEConfiguration")); // register backend processor products backEndProcessor_->setProduces(*this); } @@ -99,7 +99,7 @@ void HGCalTriggerDigiFEReproducer::produce(edm::Event& e, const edm::EventSetup& auto fe_digis_coll = *fe_digis_handle; //now we run the emulation of the back-end processor - backEndProcessor_.run(fe_digis_coll,triggerGeometry_,e); - backEndProcessor_.putInEvent(e); - backEndProcessor_.reset(); + backEndProcessor_->run(fe_digis_coll,triggerGeometry_,e); + backEndProcessor_->putInEvent(e); + backEndProcessor_->reset(); } diff --git a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc index 35c73fb96c186..946af452c3e56 100644 --- a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc +++ b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc @@ -40,7 +40,7 @@ HGCalTriggerDigiProducer(const edm::ParameterSet& conf): inputee_(consumes(conf.getParameter("eeDigis"))), inputfh_(consumes(conf.getParameter("fhDigis"))), //inputbh_(consumes(conf.getParameter("bhDigis"))), - backEndProcessor_(conf.getParameterSet("BEConfiguration"),consumesCollector()) + backEndProcessor_(new HGCalTriggerBackendProcessor(conf.getParameterSet("BEConfiguration"),consumesCollector()) ) { //setup FE codec const edm::ParameterSet& feCodecConfig = @@ -54,7 +54,7 @@ HGCalTriggerDigiProducer(const edm::ParameterSet& conf): produces(); //setup BE processor - backEndProcessor_ = std::make_unique(conf.getParameterSet("BEConfiguration")); + //backEndProcessor_ = std::make_unique(conf.getParameterSet("BEConfiguration")); // register backend processor products backEndProcessor_->setProduces(*this); } @@ -132,7 +132,8 @@ void HGCalTriggerDigiProducer::produce(edm::Event& e, const edm::EventSetup& es) auto fe_digis_coll = *fe_digis_handle; //now we run the emulation of the back-end processor - backEndProcessor_.run(fe_digis_coll,triggerGeometry_,e); - backEndProcessor_.putInEvent(e); - backEndProcessor_.reset(); + backEndProcessor_->reset(); + backEndProcessor_->run(fe_digis_coll,triggerGeometry_,e); + backEndProcessor_->putInEvent(e); + //backEndProcessor_->reset(); } diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc index d47ff04ad4538..2fabca29535eb 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc @@ -11,7 +11,7 @@ class FullModuleSumAlgo : public Algorithm public: FullModuleSumAlgo(const edm::ParameterSet& conf, edm::ConsumesCollector& cc): - Algorithm(conf,cc), + Algorithm(conf,cc), cluster_product_( new l1t::HGCalClusterBxCollection ){} virtual void setProduces(edm::EDProducer& prod) const override final @@ -20,7 +20,7 @@ class FullModuleSumAlgo : public Algorithm } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom, + const edm::ESHandle& geom, const edm::Event&evt ) override final; @@ -41,7 +41,7 @@ class FullModuleSumAlgo : public Algorithm /*****************************************************************/ void FullModuleSumAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom, + const edm::ESHandle& geom, const edm::Event&evt ) /*****************************************************************/ diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 5afc0e23e5b61..e274beece9d4f 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -1,7 +1,7 @@ // HGCal Trigger #include "L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h" //#include "L1Trigger/L1THGCal/interface/fe_codecs/HGCalBestChoiceCodec.h" -#include "L1Trigger/L1THGCal/interface/fe_codecs/HGCalTriggerBestChoiceCodec.h" +#include "L1Trigger/L1THGCal/interface/fe_codecs/HGCalTriggerCellBestChoiceCodec.h" #include "DataFormats/ForwardDetId/interface/HGCTriggerDetId.h" // HGCalClusters and detId @@ -90,9 +90,9 @@ namespace HGCalTriggerBackend{ //Consumes tokens HGCalTriggerSimCluster(const edm::ParameterSet& conf,edm::ConsumesCollector&cc) : Algorithm(conf,cc) { // I need to consumes the PF Cluster Collection with the sim clustering, TODO: make it configurable (?) - // vector "mix" "MergedCaloTruth" "HLT" + // vector "mix" "MergedCaloTruth" "HLT/DIGI" // pf clusters cannot be safely cast to SimCluster - sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("MergedCaloTruth")); + sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","DIGI")); } // setProduces @@ -104,18 +104,31 @@ namespace HGCalTriggerBackend{ // putInEvent virtual void putInEvent(edm::Event& evt) override final { +#ifdef DEBUG + cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] start"<&geom, + const edm::ESHandle&geom, const edm::Event&evt ) { @@ -166,6 +179,8 @@ namespace HGCalTriggerBackend{ const HGCalDetId tcellId(triggercell.detId()); //uint32_t digiEnergy = data.payload; i auto digiEnergy=triggercell.p4().E(); + double eta=triggercell.p4().Eta(); + double phi=triggercell.p4().Phi(); //2.B get the HGCAL-base-cell associated to it / geometry //const auto& tc=geom->triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& //for(const auto& cell : tc.components() ) // HGcell -- unsigned @@ -182,7 +197,7 @@ namespace HGCalTriggerBackend{ //2.D add to the corresponding cluster //void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float & energy,float & eta, float &phi) //addToCluster(cluster_container, pid, 0 energy,ETA/PHI? ) ; - addToCluster(cluster_container, pid, 0,energy,0.,0. ) ; // how do I get eta, phi w/o the hgcal geometry? + addToCluster(cluster_container, pid, 0,energy,eta,phi ) ; // how do I get eta, phi w/o the hgcal geometry? } } } //end of for-loop @@ -202,12 +217,25 @@ namespace HGCalTriggerBackend{ // clusterEnergyHw, clusterEtaHw, 0); // for (auto& p : cluster_container) + { + //std::unordered_map > + #ifdef DEBUG + cout<<"[HGCalTriggerSimCluster]::[run] Cluster: pid="<< p.first<push_back(p.second.first,p.second.second); // bx,cluster + } - } +#ifdef DEBUG + cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] cluster product (!=NULL) "< HGCalTriggerSimClusterBestChoice; //DEFINE_EDM_PLUGIN(HGCalTriggerBackendAlgorithmFactory, HGCalTriggerSimClusterBestChoice,"HGCalTriggerSimClusterBestChoice"); -typedef HGCalTriggerBackend::HGCalTriggerSimCluster HGCalTriggerSimClusterBestChoice; +typedef HGCalTriggerBackend::HGCalTriggerSimCluster HGCalTriggerSimClusterBestChoice; DEFINE_EDM_PLUGIN(HGCalTriggerBackendAlgorithmFactory, HGCalTriggerSimClusterBestChoice,"HGCalTriggerSimClusterBestChoice"); //DEFINE_EDM_PLUGIN(HGCalTriggerBackendAlgorithmFactory, HGCalTriggerSimCluster,"HGCalTriggerSimCluster"); diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc index cb026c03b00cb..7344b7957f782 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc @@ -18,7 +18,7 @@ class RandomClusterAlgo : public Algorithm { } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom, + const edm::ESHandle& geom, const edm::Event&evt ) override final; @@ -36,7 +36,7 @@ class RandomClusterAlgo : public Algorithm { }; void RandomClusterAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom, + const edm::ESHandle& geom, const edm::Event&evt ) { for( const auto& digi : coll ) { diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc index d7638299fb9c6..b2fd041cd651a 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc @@ -10,8 +10,8 @@ class SingleCellClusterAlgo : public Algorithm { public: - SingleCellClusterAlgo(const edm::ParameterSet& conf): - Algorithm(conf), + SingleCellClusterAlgo(const edm::ParameterSet& conf, edm::ConsumesCollector &cc ): + Algorithm(conf,cc), cluster_product_( new l1t::HGCalClusterBxCollection ){} virtual void setProduces(edm::EDProducer& prod) const override final @@ -20,7 +20,7 @@ class SingleCellClusterAlgo : public Algorithm } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom, + const edm::ESHandle& geom, const edm::Event&evt ) override final; @@ -41,7 +41,7 @@ class SingleCellClusterAlgo : public Algorithm /*****************************************************************/ void SingleCellClusterAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom, + const edm::ESHandle& geom, const edm::Event&evt ) /*****************************************************************/ diff --git a/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc b/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc index 7e93fec177e13..eb3318542734b 100644 --- a/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc +++ b/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc @@ -28,7 +28,8 @@ void HGCalTriggerBackendProcessor::setProduces(edm::EDProducer& prod) const { void HGCalTriggerBackendProcessor:: run(const l1t::HGCFETriggerDigiCollection& coll, - const std::unique_ptr& geom, + //const std::unique_ptr& geom, + const edm::ESHandle & geom, const edm::Event &e) { for( auto& algo : algorithms_ ) { algo->run(coll,geom,e); From 16f79e6095815ef930436d09be570b9fc99656c3 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Wed, 2 Nov 2016 13:48:53 +0100 Subject: [PATCH 09/21] requested changes --- .../HGCalTriggerBackendAlgorithmBase.h | 2 +- .../be_algorithms/FullModuleSumAlgo.cc | 2 -- .../plugins/be_algorithms/HGCalSimCluster.cc | 30 +++++++++++-------- .../be_algorithms/RandomClusterAlgo.cc | 2 -- .../be_algorithms/SingleCellClusterAlgo.cc | 2 -- .../src/HGCalTriggerBackendProcessor.cc | 2 +- 6 files changed, 20 insertions(+), 20 deletions(-) diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h index c50e10aede028..ae00c44f35d36 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h @@ -38,6 +38,7 @@ class HGCalTriggerBackendAlgorithmBase { // Allow HGCalTriggerBackend to be passed a consume collector HGCalTriggerBackendAlgorithmBase(const edm::ParameterSet& conf, edm::ConsumesCollector &cc) : + geometry_(nullptr), name_(conf.getParameter("AlgorithmName")) {} @@ -51,7 +52,6 @@ class HGCalTriggerBackendAlgorithmBase { virtual void setProduces(edm::EDProducer& prod) const = 0; virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const edm::ESHandle& geom, const edm::Event &e ) = 0; diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc index 2fabca29535eb..8b9016c30ab1c 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc @@ -20,7 +20,6 @@ class FullModuleSumAlgo : public Algorithm } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const edm::ESHandle& geom, const edm::Event&evt ) override final; @@ -41,7 +40,6 @@ class FullModuleSumAlgo : public Algorithm /*****************************************************************/ void FullModuleSumAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, - const edm::ESHandle& geom, const edm::Event&evt ) /*****************************************************************/ diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index e274beece9d4f..8c3a24c4f8e27 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -49,14 +49,17 @@ namespace HGCalTriggerBackend{ // add to cluster void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float energy,float eta, float phi) { - auto iterator = cluster_container.find (pid); - if (iterator == cluster_container.end()) - { - //create an empty cluster - cluster_container[pid] = std::pair(0,l1t::HGCalCluster()); - iterator = cluster_container.find (pid); - iterator -> second . second . setPdgId(pdgid); - } + + auto pair = cluster_container.emplace(pid, std::pair(0,l1t::HGCalCluster() ) ) ; + auto iterator = pair.first; + //auto iterator = cluster_container.find (pid); + //if (iterator == cluster_container.end()) + // { + // //create an empty cluster + // cluster_container[pid] = std::pair(0,l1t::HGCalCluster()); + // iterator = cluster_container.find (pid); + // iterator -> second . second . setPdgId(pdgid); + // } // p4 += p4' math::PtEtaPhiMLorentzVectorD p4; p4.SetPt ( iterator -> second . second . pt() ) ; @@ -74,6 +77,7 @@ namespace HGCalTriggerBackend{ return ; } + using Algorithm::geometry_; protected: using Algorithm::codec_; @@ -128,7 +132,6 @@ namespace HGCalTriggerBackend{ // run, actual algorithm virtual void run( const l1t::HGCFETriggerDigiCollection & coll, - const edm::ESHandle&geom, const edm::Event&evt ) { @@ -167,7 +170,7 @@ namespace HGCalTriggerBackend{ DATA data; digi.decode(codec_,data); //2.A get the trigger-cell information energy/id - const HGCTriggerDetId& moduleId = digi.getDetId(); // this is a module Det Id + //const HGCTriggerDetId& moduleId = digi.getDetId(); // this is a module Det Id // there is a loss of generality here, due to the restriction imposed by the data formats // it will work if inside a module there is a data.payload with an ordered list of all the energies @@ -184,11 +187,14 @@ namespace HGCalTriggerBackend{ //2.B get the HGCAL-base-cell associated to it / geometry //const auto& tc=geom->triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& //for(const auto& cell : tc.components() ) // HGcell -- unsigned - for(const auto& cell : geom->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned + for(const auto& cell : geometry_->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned { HGCalDetId cellId(cell); //2.C get the particleId and energy fractions - const auto& particles = simclusters[cellId]; // vector pid fractions + const auto & iterator= simclusters.find(cellId); + if (iterator == simclusters.end() ) continue; + //const auto& particles = simclusters[cellId]; // vector pid fractions + const auto & particles = iterator->second; for ( const auto& p: particles ) { const auto & pid= p.first; diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc index 7344b7957f782..497747d70683a 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc @@ -18,7 +18,6 @@ class RandomClusterAlgo : public Algorithm { } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const edm::ESHandle& geom, const edm::Event&evt ) override final; @@ -36,7 +35,6 @@ class RandomClusterAlgo : public Algorithm { }; void RandomClusterAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, - const edm::ESHandle& geom, const edm::Event&evt ) { for( const auto& digi : coll ) { diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc index b2fd041cd651a..919f157dfe262 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc @@ -20,7 +20,6 @@ class SingleCellClusterAlgo : public Algorithm } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const edm::ESHandle& geom, const edm::Event&evt ) override final; @@ -41,7 +40,6 @@ class SingleCellClusterAlgo : public Algorithm /*****************************************************************/ void SingleCellClusterAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, - const edm::ESHandle& geom, const edm::Event&evt ) /*****************************************************************/ diff --git a/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc b/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc index eb3318542734b..a5341ceda4f32 100644 --- a/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc +++ b/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc @@ -32,7 +32,7 @@ run(const l1t::HGCFETriggerDigiCollection& coll, const edm::ESHandle & geom, const edm::Event &e) { for( auto& algo : algorithms_ ) { - algo->run(coll,geom,e); + algo->run(coll,e); } } From a5873a3d1af57ed65e8b84c3583881e1a5fb1a6b Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Wed, 2 Nov 2016 13:59:17 +0100 Subject: [PATCH 10/21] avoid double counting of energy --- L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 8c3a24c4f8e27..34294af63a2c2 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -187,6 +187,7 @@ namespace HGCalTriggerBackend{ //2.B get the HGCAL-base-cell associated to it / geometry //const auto& tc=geom->triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& //for(const auto& cell : tc.components() ) // HGcell -- unsigned + unsigned ncells = geometry_->getCellsFromTriggerCell( tcellId()).size(); for(const auto& cell : geometry_->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned { HGCalDetId cellId(cell); @@ -199,7 +200,7 @@ namespace HGCalTriggerBackend{ { const auto & pid= p.first; const auto & fraction=p.second; - auto energy = fraction*digiEnergy; + auto energy = fraction*digiEnergy/ncells; //2.D add to the corresponding cluster //void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float & energy,float & eta, float &phi) //addToCluster(cluster_container, pid, 0 energy,ETA/PHI? ) ; From ede54e16afbee0218d5316f8899defb2bc7f5946 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Mon, 21 Nov 2016 16:56:53 +0100 Subject: [PATCH 11/21] trying to debug --- .../plugins/HGCalTriggerDigiProducer.cc | 27 ++++++++++ .../plugins/be_algorithms/HGCalSimCluster.cc | 28 +++++----- .../L1THGCal/test/testTriggerSimCluster.py | 53 +++++++++++++++++++ 3 files changed, 96 insertions(+), 12 deletions(-) create mode 100644 L1Trigger/L1THGCal/test/testTriggerSimCluster.py diff --git a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc index 946af452c3e56..5bd8d1505931b 100644 --- a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc +++ b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc @@ -15,6 +15,13 @@ #include #include +#define HGCAL_DEBUG + +#ifdef HGCAL_DEBUG + #include + using namespace std; +#endif + class HGCalTriggerDigiProducer : public edm::EDProducer { public: HGCalTriggerDigiProducer(const edm::ParameterSet&); @@ -42,6 +49,9 @@ HGCalTriggerDigiProducer(const edm::ParameterSet& conf): //inputbh_(consumes(conf.getParameter("bhDigis"))), backEndProcessor_(new HGCalTriggerBackendProcessor(conf.getParameterSet("BEConfiguration"),consumesCollector()) ) { + #ifdef HGCAL_DEBUG + cout<<"[HGCalTriggerDigiProducer]::["<< __FUNCTION__ <<"] constructor: file " <<__FILE__<<":"<<__LINE__<(conf.getParameterSet("BEConfiguration")); // register backend processor products backEndProcessor_->setProduces(*this); + + #ifdef HGCAL_DEBUG + cout<<"[HGCalTriggerDigiProducer]::["<< __FUNCTION__ <<"] end: file " <<__FILE__<<":"<<__LINE__<().get(triggerGeometry_); codec_->setGeometry(triggerGeometry_.product()); backEndProcessor_->setGeometry(triggerGeometry_.product()); + #ifdef HGCAL_DEBUG + cout<<"[HGCalTriggerDigiProducer]::["<< __FUNCTION__ <<"] end" < fe_output( new l1t::HGCFETriggerDigiCollection ); @@ -136,4 +160,7 @@ void HGCalTriggerDigiProducer::produce(edm::Event& e, const edm::EventSetup& es) backEndProcessor_->run(fe_digis_coll,triggerGeometry_,e); backEndProcessor_->putInEvent(e); //backEndProcessor_->reset(); + #ifdef HGCAL_DEBUG + cout<<"[HGCalTriggerDigiProducer]::["<< __FUNCTION__ <<"] end. file " <<__FILE__<<":"<<__LINE__< "mix" "MergedCaloTruth" "HLT/DIGI" // pf clusters cannot be safely cast to SimCluster - sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","DIGI")); + //sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","DIGI")); + sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","")); } // setProduces virtual void setProduces(edm::EDProducer& prod) const override final { +#ifdef HGCAL_DEBUG + cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] start: Setting producer to produce:"<(name()); } // putInEvent virtual void putInEvent(edm::Event& evt) override final { -#ifdef DEBUG +#ifdef HGCAL_DEBUG cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] start"< [ (pid, fraction), ... @@ -160,7 +164,7 @@ namespace HGCalTriggerBackend{ } } -#ifdef DEBUG +#ifdef HGCAL_DEBUG cout<<"[HGCalTriggerSimCluster]::[run] Run on digis"< > - #ifdef DEBUG + #ifdef HGCAL_DEBUG cout<<"[HGCalTriggerSimCluster]::[run] Cluster: pid="<< p.first<push_back(p.second.first,p.second.second); // bx,cluster } -#ifdef DEBUG +#ifdef HGCAL_DEBUG cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] cluster product (!=NULL) "< "hgcalTriggerPrimitiveDigiProducer" "HGCalTriggerSimClusterBestChoice" "DIGI" +#clusters, clusterLabel = Handle("BXVector"), "hgcalTriggerPrimitiveDigiProducer:HGCalTriggerSimClusterBestChoice" +clusters, clusterLabel = Handle("l1t::HGCalClusterBxCollection"), "hgcalTriggerPrimitiveDigiProducer:HGCalTriggerSimClusterBestChoice" + +verbose=True + +for i,event in enumerate(events): + if verbose: + print "\nEvent", i + if i > 10: break + + event.getByLabel(clusterLabel, clusters) + + ## check that handle is valid + if not clusters.isValid(): print "Cluster handle is invalid" + + ## -- + if verbose: print "-> BX=0" + if verbose: print "-> clusters=",clusters,"\n","product=",clusters.product() + + print "what's saved in clusters?" + for bx in clusters.product(): + print "bx=",bx + print "product", clusters.product()[bx] + + ## get 0 bunch crossing vector + for bx in [0]: + + if verbose: print "-> 0 pos",clusters.product()[0] + + vector = clusters.product()[bx] + for idx in range(0,vector.size()): + print "* accessing position ", idx + print " ",vector[idx].pt(),vector[idx].eta() + From c0385472840bdfb6e76949d9aad41970c5739b80 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Wed, 14 Dec 2016 12:20:50 +0100 Subject: [PATCH 12/21] merging conflicts p2 --- .../HGCalTriggerBackendAlgorithmBase.h | 3 +- .../be_algorithms/FullModuleSumAlgo.cc | 18 +++---- .../plugins/be_algorithms/HGCalSimCluster.cc | 51 ++++++++++--------- .../be_algorithms/RandomClusterAlgo.cc | 2 +- .../be_algorithms/SingleCellClusterAlgo.cc | 2 +- 5 files changed, 39 insertions(+), 37 deletions(-) diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h index fa335b75b7c49..df62746c85287 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h @@ -51,7 +51,8 @@ class HGCalTriggerBackendAlgorithmBase { //runs the trigger algorithm, storing internally the results virtual void setProduces(edm::EDProducer& prod) const = 0; - virtual void run(const l1t::HGCFETriggerDigiCollection& coll, edm::EventSetup& es, + virtual void run(const l1t::HGCFETriggerDigiCollection& coll, + const edm::EventSetup& es, const edm::Event &e ) = 0; diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc index 774566f3fd11f..975b5a9fbdce0 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/FullModuleSumAlgo.cc @@ -18,13 +18,11 @@ class FullModuleSumAlgo : public Algorithm using Algorithm::codec_; public: - //FullModuleSumAlgo(const edm::ParameterSet& conf): - //Algorithm(conf), - //cluster_product_( new l1t::HGCalClusterBxCollection ){} FullModuleSumAlgo(const edm::ParameterSet& conf, edm::ConsumesCollector& cc): - Algorithm(conf,cc), - cluster_product_( new l1t::HGCalClusterBxCollection ){} + Algorithm(conf,cc), + cluster_product_( new l1t::HGCalClusterBxCollection ) + {} virtual void setProduces(edm::EDProducer& prod) const override final { @@ -48,14 +46,16 @@ class FullModuleSumAlgo : public Algorithm }; /*****************************************************************/ -void FullModuleSumAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, const edm::EventSetup& es, - const edm::Event&evt - ) +template +void FullModuleSumAlgo::run(const l1t::HGCFETriggerDigiCollection& coll, + const edm::EventSetup& es, + const edm::Event&evt + ) /*****************************************************************/ { for( const auto& digi : coll ) { - HGCalTriggerCellBestChoiceCodec::data_type data; + DATA data; data.reset(); const HGCalDetId& moduleId = digi.getDetId(); digi.decode(codec_, data); diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 757b27d4d15db..6fed4d9029fe9 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -50,31 +50,31 @@ namespace HGCalTriggerBackend{ void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float energy,float eta, float phi) { - auto pair = cluster_container.emplace(pid, std::pair(0,l1t::HGCalCluster() ) ) ; - auto iterator = pair.first; - //auto iterator = cluster_container.find (pid); - //if (iterator == cluster_container.end()) - // { - // //create an empty cluster - // cluster_container[pid] = std::pair(0,l1t::HGCalCluster()); - // iterator = cluster_container.find (pid); - // iterator -> second . second . setPdgId(pdgid); - // } - // p4 += p4' - math::PtEtaPhiMLorentzVectorD p4; - p4.SetPt ( iterator -> second . second . pt() ) ; - p4.SetEta( iterator -> second . second . eta() ) ; - p4.SetPhi( iterator -> second . second . phi() ) ; - p4.SetM ( iterator -> second . second . mass() ) ; - math::PtEtaPhiMLorentzVectorD pp4; - float t = std::exp (- eta); - pp4.SetPt ( energy * (1-t*t)/(1+t*t) ) ; - pp4.SetEta( eta ) ; - pp4.SetPhi( phi ) ; - pp4.SetM ( 0 ) ; - p4 += pp4; - iterator -> second . second . setP4(p4); - return ; + auto pair = cluster_container.emplace(pid, std::pair(0,l1t::HGCalCluster() ) ) ; + auto iterator = pair.first; + //auto iterator = cluster_container.find (pid); + //if (iterator == cluster_container.end()) + // { + // //create an empty cluster + // cluster_container[pid] = std::pair(0,l1t::HGCalCluster()); + // iterator = cluster_container.find (pid); + // iterator -> second . second . setPdgId(pdgid); + // } + // p4 += p4' + math::PtEtaPhiMLorentzVectorD p4; + p4.SetPt ( iterator -> second . second . pt() ) ; + p4.SetEta( iterator -> second . second . eta() ) ; + p4.SetPhi( iterator -> second . second . phi() ) ; + p4.SetM ( iterator -> second . second . mass() ) ; + math::PtEtaPhiMLorentzVectorD pp4; + float t = std::exp (- eta); + pp4.SetPt ( energy * (1-t*t)/(1+t*t) ) ; + pp4.SetEta( eta ) ; + pp4.SetPhi( phi ) ; + pp4.SetM ( 0 ) ; + p4 += pp4; + iterator -> second . second . setP4(p4); + return ; } using Algorithm::geometry_; @@ -136,6 +136,7 @@ namespace HGCalTriggerBackend{ // run, actual algorithm virtual void run( const l1t::HGCFETriggerDigiCollection & coll, + const edm::EventSetup& es, const edm::Event&evt ) { diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc index d83abc77e85d3..a0ace45e7a0cd 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/RandomClusterAlgo.cc @@ -18,7 +18,7 @@ class RandomClusterAlgo : public Algorithm { } virtual void run(const l1t::HGCFETriggerDigiCollection& coll, - const edm::EventSetup& es + const edm::EventSetup& es, const edm::Event&evt ) override final; diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc index 6ee593f6153b8..2826afdc006c7 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/SingleCellClusterAlgo.cc @@ -18,7 +18,7 @@ class SingleCellClusterAlgo : public Algorithm public: SingleCellClusterAlgo(const edm::ParameterSet& conf,edm::ConsumesCollector &cc): - Algorithm(conf,cc), + Algorithm(conf,cc), cluster_product_( new l1t::HGCalTriggerCellBxCollection ), HGCalEESensitive_(conf.getParameter("HGCalEESensitive_tag")), HGCalHESiliconSensitive_(conf.getParameter("HGCalHESiliconSensitive_tag")), From 5c5473afcae6c00ca425cf831d5a6246899ec8c0 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Wed, 14 Dec 2016 12:37:02 +0100 Subject: [PATCH 13/21] trailing slash in sim clusternig --- .../plugins/be_algorithms/HGCalSimCluster.cc | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 6fed4d9029fe9..fb79011bc24a5 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -41,11 +41,11 @@ namespace HGCalTriggerBackend{ class HGCalTriggerSimCluster : public Algorithm { private: - std::unique_ptr cluster_product; + std::unique_ptr cluster_product_; // token - edm::EDGetTokenT< std::vector > sim_token; + edm::EDGetTokenT< std::vector > sim_token_; // handle - edm::Handle< std::vector > sim_handle; + edm::Handle< std::vector > sim_handle_; // add to cluster void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float energy,float eta, float phi) { @@ -96,8 +96,8 @@ namespace HGCalTriggerBackend{ // I need to consumes the PF Cluster Collection with the sim clustering, TODO: make it configurable (?) // vector "mix" "MergedCaloTruth" "HLT/DIGI" // pf clusters cannot be safely cast to SimCluster - //sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","DIGI")); - sim_token = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","")); + //sim_token_ = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","DIGI")); + sim_token_ = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","")); } // setProduces @@ -114,9 +114,9 @@ namespace HGCalTriggerBackend{ { #ifdef HGCAL_DEBUG cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] start"< > cluster_container;// PID-> bx,cluster - evt.getByToken(sim_token,sim_handle); + evt.getByToken(sim_token_,sim_handle_); - if (not sim_handle.isValid()) { std::cout<<"[HGCalTriggerSimCluster]::[run]::[ERROR] PFCluster collection for HGC sim clustering not available"< [ (pid, fraction), ... std::unordered_map > > simclusters; - for (auto& cluster : *sim_handle) + for (auto& cluster : *sim_handle_) { auto pid= cluster.particleId(); // not pdgId const auto& hf = cluster.hits_and_fractions(); @@ -219,7 +219,7 @@ namespace HGCalTriggerBackend{ cout<<"[HGCalTriggerSimCluster]::[run] Push clusters in cluster products"<getTriggerCellPosition( tcellId()); @@ -237,11 +237,11 @@ namespace HGCalTriggerBackend{ cout<<"[HGCalTriggerSimCluster]::[run] Cluster: ----------- l1t eta"<< p.second.second.eta() <push_back(p.second.first,p.second.second); // bx,cluster + cluster_product_->push_back(p.second.first,p.second.second); // bx,cluster } #ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] cluster product (!=NULL) "< Date: Wed, 14 Dec 2016 12:42:40 +0100 Subject: [PATCH 14/21] moving to cmsExceptions --- L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index fb79011bc24a5..d6cd2522aadf2 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -147,7 +147,9 @@ namespace HGCalTriggerBackend{ std::unordered_map > cluster_container;// PID-> bx,cluster evt.getByToken(sim_token_,sim_handle_); - if (not sim_handle_.isValid()) { std::cout<<"[HGCalTriggerSimCluster]::[run]::[ERROR] PFCluster collection for HGC sim clustering not available"< Date: Mon, 19 Dec 2016 13:13:51 +0100 Subject: [PATCH 15/21] adding energy calibration --- .../plugins/be_algorithms/HGCalSimCluster.cc | 49 +++++++++++++++++-- .../L1THGCal/test/testTriggerSimCluster.py | 46 ++++++++++++----- 2 files changed, 77 insertions(+), 18 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index d6cd2522aadf2..6e9dd83661c86 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -20,8 +20,12 @@ // #include "DataFormats/Math/interface/LorentzVector.h" +// Energy calibration +#include "L1Trigger/L1THGCal/interface/be_algorithms/HGCalTriggerCellCalibration.h" + // Print out something before crashing or throwing exceptions #include +#include /* Original Author: Andrea Carlo Marini @@ -42,10 +46,21 @@ namespace HGCalTriggerBackend{ { private: std::unique_ptr cluster_product_; - // token - edm::EDGetTokenT< std::vector > sim_token_; // handle edm::Handle< std::vector > sim_handle_; + // calibration handle + edm::ESHandle hgceeTopoHandle_; + edm::ESHandle hgchefTopoHandle_; + + // variables that needs to be init, in the right order + // energy calibration + std::string HGCalEESensitive_; + std::string HGCalHESiliconSensitive_; + HGCalTriggerCellCalibration calibration_; + + // token + edm::EDGetTokenT< std::vector > sim_token_; + // add to cluster void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float energy,float eta, float phi) { @@ -92,7 +107,12 @@ namespace HGCalTriggerBackend{ using Algorithm::reset; //Consumes tokens - HGCalTriggerSimCluster(const edm::ParameterSet& conf,edm::ConsumesCollector&cc) : Algorithm(conf,cc) { + HGCalTriggerSimCluster(const edm::ParameterSet& conf,edm::ConsumesCollector&cc) : + Algorithm(conf,cc), + HGCalEESensitive_(conf.getParameter("HGCalEESensitive_tag")), + HGCalHESiliconSensitive_(conf.getParameter("HGCalHESiliconSensitive_tag")), + calibration_(conf.getParameterSet("calib_parameters")) + { // I need to consumes the PF Cluster Collection with the sim clustering, TODO: make it configurable (?) // vector "mix" "MergedCaloTruth" "HLT/DIGI" // pf clusters cannot be safely cast to SimCluster @@ -150,6 +170,9 @@ namespace HGCalTriggerBackend{ if (not sim_handle_.isValid()){ throw cms::Exception("ContentError")<<"[HGCalTriggerSimCluster]::[run]::[ERROR] PFCluster collection for HGC sim clustering not available"; } + // calibration + es.get().get(HGCalEESensitive_, hgceeTopoHandle_); + es.get().get(HGCalHESiliconSensitive_, hgchefTopoHandle_); // 1.5. pre-process the sim cluster to have easy accessible information #ifdef HGCAL_DEBUG @@ -187,8 +210,22 @@ namespace HGCalTriggerBackend{ if(triggercell.hwPt()<=0) continue; const HGCalDetId tcellId(triggercell.detId()); + // calbration + int subdet = tcellId.subdetId(); + int cellThickness = 0; + + if( subdet == HGCEE ){ + cellThickness = (hgceeTopoHandle_)->dddConstants().waferTypeL((unsigned int)tcellId.wafer() ); + }else if( subdet == HGCHEF ){ + cellThickness = (hgchefTopoHandle_)->dddConstants().waferTypeL((unsigned int)tcellId.wafer() ); + }else if( subdet == HGCHEB ){ + edm::LogWarning("DataNotFound") << "ATTENTION: the BH trgCells are not yet implemented !! "; + } + l1t::HGCalTriggerCell calibratedtriggercell(triggercell); + calibration_.calibrate(calibratedtriggercell, cellThickness); //uint32_t digiEnergy = data.payload; i - auto digiEnergy=triggercell.p4().E(); + //auto digiEnergy=triggercell.p4().E(); + auto calibratedDigiEnergy=calibratedtriggercell.p4().E(); double eta=triggercell.p4().Eta(); double phi=triggercell.p4().Phi(); //2.B get the HGCAL-base-cell associated to it / geometry @@ -207,7 +244,9 @@ namespace HGCalTriggerBackend{ { const auto & pid= p.first; const auto & fraction=p.second; - auto energy = fraction*digiEnergy/ncells; + //auto energy = fraction*digiEnergy/ncells; + auto energy = fraction * calibratedDigiEnergy/ncells; + //2.D add to the corresponding cluster //void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float & energy,float & eta, float &phi) //addToCluster(cluster_container, pid, 0 energy,ETA/PHI? ) ; diff --git a/L1Trigger/L1THGCal/test/testTriggerSimCluster.py b/L1Trigger/L1THGCal/test/testTriggerSimCluster.py index f2cf1fc3fe3e2..8883d7868ccd2 100644 --- a/L1Trigger/L1THGCal/test/testTriggerSimCluster.py +++ b/L1Trigger/L1THGCal/test/testTriggerSimCluster.py @@ -1,9 +1,9 @@ # import ROOT in batch mode import sys oldargv = sys.argv[:] -sys.argv = [ '-b-' ] +#sys.argv = [ '-b-' ] import ROOT -ROOT.gROOT.SetBatch(True) +#ROOT.gROOT.SetBatch(True) sys.argv = oldargv # load FWLite C++ libraries @@ -22,10 +22,15 @@ verbose=True -for i,event in enumerate(events): +hEta=ROOT.TH1D("eta","eta",1000,-20,20) +hPhi=ROOT.TH1D("phi","phi",1000,-3.1416,3.1416) +hPt = ROOT.TH1D("pt","pt",1000,0,1000) +hE = ROOT.TH1D("E","E",1000,0,1000) + +for iev,event in enumerate(events): if verbose: - print "\nEvent", i - if i > 10: break + print "\nEvent %d: run %6d, lumi %4d, event %12d" % (iev,event.eventAuxiliary().run(), event.eventAuxiliary().luminosityBlock(),event.eventAuxiliary().event()) + #if iev > 10: break event.getByLabel(clusterLabel, clusters) @@ -39,15 +44,30 @@ print "what's saved in clusters?" for bx in clusters.product(): print "bx=",bx - print "product", clusters.product()[bx] ## get 0 bunch crossing vector - for bx in [0]: - + for bx,vector in enumerate(clusters.product()): if verbose: print "-> 0 pos",clusters.product()[0] + if vector == None: + print " cluster product is none" + continue + print " pt=",vector.pt(),"eta=",vector.eta(),"phi=",vector.phi() + hEta.Fill(vector.eta() ) + if vector.eta() <8: + hE.Fill(vector.energy() ) + hPt.Fill(vector.pt() ) + hPhi.Fill(vector.phi() ) + + +c=ROOT.TCanvas() +c.Divide(2,2) +c.cd(1) +hPt.Draw("HIST") +c.cd(2) +hEta.Draw("HIST") +c.cd(3) +hPhi.Draw("HIST") +c.cd(4) +hE.Draw("HIST") +raw_input("ok?") - vector = clusters.product()[bx] - for idx in range(0,vector.size()): - print "* accessing position ", idx - print " ",vector[idx].pt(),vector[idx].eta() - From 8f61174b46383a109441be1a36cec454c260efee Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Mon, 19 Dec 2016 14:38:15 +0100 Subject: [PATCH 16/21] removing calling to the geometry --- L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h | 1 - L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc | 2 +- L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc | 2 +- L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc | 5 +++-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h index ff6018f0c0039..f8898c1a5dab0 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h @@ -39,7 +39,6 @@ class HGCalTriggerBackendProcessor { void run(const l1t::HGCFETriggerDigiCollection& coll, const edm::EventSetup& es, - const edm::ESHandle& geom, const edm::Event&e ); diff --git a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc index 1267ece1729eb..4835a57ccd5f9 100644 --- a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc +++ b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc @@ -99,7 +99,7 @@ void HGCalTriggerDigiFEReproducer::produce(edm::Event& e, const edm::EventSetup& auto fe_digis_coll = *fe_digis_handle; //now we run the emulation of the back-end processor - backEndProcessor_->run(fe_digis_coll,es,triggerGeometry_,e); + backEndProcessor_->run(fe_digis_coll,es,e); backEndProcessor_->putInEvent(e); backEndProcessor_->reset(); } diff --git a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc index 6da5cf225ec9b..b23ba4671247c 100644 --- a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc +++ b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc @@ -154,7 +154,7 @@ void HGCalTriggerDigiProducer::produce(edm::Event& e, const edm::EventSetup& es) //now we run the emulation of the back-end processor backEndProcessor_->reset(); - backEndProcessor_->run(fe_digis_coll,es,triggerGeometry_,e); + backEndProcessor_->run(fe_digis_coll,es,e); backEndProcessor_->putInEvent(e); //backEndProcessor_->reset(); #ifdef HGCAL_DEBUG diff --git a/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc b/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc index 1e27ecfa45dbe..3ca18b3165bc9 100644 --- a/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc +++ b/L1Trigger/L1THGCal/src/HGCalTriggerBackendProcessor.cc @@ -26,8 +26,9 @@ void HGCalTriggerBackendProcessor::setProduces(edm::EDProducer& prod) const { } } -void HGCalTriggerBackendProcessor::run(const l1t::HGCFETriggerDigiCollection& coll, const edm::EventSetup& es, - const edm::ESHandle & geom,const edm::Event &e +void HGCalTriggerBackendProcessor::run(const l1t::HGCFETriggerDigiCollection& coll, + const edm::EventSetup& es, + const edm::Event &e ) { for( auto& algo : algorithms_ ) { algo->run(coll, es,e); From d586b56f75bdbc4ed12f4fec6fa77f844bb49746 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Sun, 25 Dec 2016 18:48:15 +0100 Subject: [PATCH 17/21] recompute normalization. just to be safe for the time being. --- .../plugins/be_algorithms/HGCalSimCluster.cc | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 6e9dd83661c86..ef2b3876a30ea 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -232,6 +232,25 @@ namespace HGCalTriggerBackend{ //const auto& tc=geom->triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& //for(const auto& cell : tc.components() ) // HGcell -- unsigned unsigned ncells = geometry_->getCellsFromTriggerCell( tcellId()).size(); + + // normalization loop -- this loop is fast ~4 elem + double norm=0.0; + for(const auto& cell : geometry_->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned + { + HGCalDetId cellId(cell); + //2.C get the particleId and energy fractions + const auto & iterator= simclusters.find(cellId); + if (iterator == simclusters.end() ) continue; + //const auto& particles = simclusters[cellId]; // vector pid fractions + const auto & particles = iterator->second; + for ( const auto& p: particles ) + { + const auto & pid= p.first; + const auto & fraction=p.second; + norm += fraction; + } + } + // for(const auto& cell : geometry_->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned { HGCalDetId cellId(cell); @@ -245,7 +264,7 @@ namespace HGCalTriggerBackend{ const auto & pid= p.first; const auto & fraction=p.second; //auto energy = fraction*digiEnergy/ncells; - auto energy = fraction * calibratedDigiEnergy/ncells; + auto energy = fraction * calibratedDigiEnergy/norm; //2.D add to the corresponding cluster //void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float & energy,float & eta, float &phi) From d1b7f28928bbb1fb2447cc560eced7a9965d39d4 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Wed, 25 Jan 2017 09:49:39 +0100 Subject: [PATCH 18/21] cluster shapes --- .../L1THGCal/interface/ClusterShapes.h | 85 ++++++++++ DataFormats/L1THGCal/interface/HGCalCluster.h | 6 +- DataFormats/L1THGCal/src/ClusterShapes.cc | 146 ++++++++++++++++++ DataFormats/L1THGCal/src/HGCalCluster.cc | 1 + DataFormats/L1THGCal/src/classes.h | 4 + DataFormats/L1THGCal/src/classes_def.xml | 2 + .../plugins/be_algorithms/HGCalSimCluster.cc | 45 ++++++ 7 files changed, 288 insertions(+), 1 deletion(-) create mode 100644 DataFormats/L1THGCal/interface/ClusterShapes.h create mode 100644 DataFormats/L1THGCal/src/ClusterShapes.cc diff --git a/DataFormats/L1THGCal/interface/ClusterShapes.h b/DataFormats/L1THGCal/interface/ClusterShapes.h new file mode 100644 index 0000000000000..7fb33f284f292 --- /dev/null +++ b/DataFormats/L1THGCal/interface/ClusterShapes.h @@ -0,0 +1,85 @@ +#ifndef CLUSTER_SHAPES_H +#define CLUSTER_SHAPES_H +#include +#include + +namespace l1t{ + // this class is design to contain and compute + // efficiently the cluster shapes + // running only once on the cluster members. + class ClusterShapes{ + private: + float sum_e = 0.0; + float sum_e2 = 0.0; + float sum_logE = 0.0; + int n=0.0; + + float emax = 0.0; + + float sum_w =0.0; // just here for clarity + float sum_eta = 0.0; + float sum_r = 0.0; + // i will discriminate using the rms in -pi,pi or in 0,pi + float sum_phi_0 = 0.0; // computed in -pi,pi + float sum_phi_1 = 0.0; // computed in 0, 2pi + + float sum_eta2=0.0; + float sum_r2 = 0.0; + float sum_phi2_0=0.0; //computed in -pi,pi + float sum_phi2_1=0.0; //computed in 0,2pi + + // off diagonal element of the tensor + float sum_eta_r =0.0; + float sum_r_phi_0 = 0.0; + float sum_r_phi_1 = 0.0; + float sum_eta_phi_0 = 0.0; + float sum_eta_phi_1 = 0.0; + + // caching of informations + mutable bool isPhi0 = true; + mutable bool modified_ = false; // check wheneever i need + public: + ClusterShapes(){} + ClusterShapes(float e, float eta, float phi, float r) { Init(e,eta,phi,r);} + ~ClusterShapes(){} + ClusterShapes(const ClusterShapes&x)= default; + //init an empty cluster + void Init(float e, float eta, float phi, float r=0.); + inline void Add(float e, float eta, float phi, float r=0.0) + { (*this) += ClusterShapes(e,eta,phi,r);} + + + // --- this is what I want out: + float SigmaEtaEta() const ; + float SigmaPhiPhi() const ; + float SigmaRR() const ; + // ---- + float Phi() const ; + float R() const ; + float Eta() const ; + inline int N()const {return n;} + // -- + float SigmaEtaR()const ; + float SigmaEtaPhi() const ; + float SigmaRPhi() const ; + // -- + float LogEoverE() const { return sum_logE/sum_e;} + float eD() const { return std::sqrt(sum_e2)/sum_e;} + + ClusterShapes operator+(const ClusterShapes &); + void operator+=(const ClusterShapes &); + ClusterShapes& operator=(const ClusterShapes &) = default; + }; + +}; // end namespace + +#endif + +// ----------------------------------- +// Local Variables: +// mode:c++ +// indent-tabs-mode:nil +// tab-width:4 +// c-basic-offset:4 +// End: +// vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 diff --git a/DataFormats/L1THGCal/interface/HGCalCluster.h b/DataFormats/L1THGCal/interface/HGCalCluster.h index 74cad31fe9e30..fa26b0a251c00 100644 --- a/DataFormats/L1THGCal/interface/HGCalCluster.h +++ b/DataFormats/L1THGCal/interface/HGCalCluster.h @@ -3,9 +3,11 @@ #include "DataFormats/L1Trigger/interface/L1Candidate.h" #include "DataFormats/L1Trigger/interface/BXVector.h" +#include "DataFormats/L1THGCal/interface/ClusterShapes.h" namespace l1t { - + + class HGCalCluster : public L1Candidate { public: HGCalCluster(){} @@ -36,6 +38,8 @@ namespace l1t { uint32_t hOverE() const {return hOverE_;} + ClusterShapes shapes ; + bool operator<(const HGCalCluster& cl) const; bool operator>(const HGCalCluster& cl) const {return cl<*this;}; bool operator<=(const HGCalCluster& cl) const {return !(cl>*this);}; diff --git a/DataFormats/L1THGCal/src/ClusterShapes.cc b/DataFormats/L1THGCal/src/ClusterShapes.cc new file mode 100644 index 0000000000000..e72f86012e8d5 --- /dev/null +++ b/DataFormats/L1THGCal/src/ClusterShapes.cc @@ -0,0 +1,146 @@ +#include "DataFormats/L1THGCal/interface/ClusterShapes.h" +#include + +using namespace l1t; + +ClusterShapes ClusterShapes::operator+(const ClusterShapes& x) +{ + ClusterShapes cs(*this); // copy constructor + cs += x; + return cs; +} + + +void ClusterShapes::operator +=(const ClusterShapes &x){ + + sum_e += x.sum_e; + sum_e2 += x.sum_e2; + sum_logE += x.sum_logE; + n += x.n; + + sum_w += x.sum_w; + + emax = (emax> x.emax) ? emax: x.emax; + + // mid-point + sum_eta += x.sum_eta; + sum_phi_0 += x.sum_phi_0; // + sum_phi_1 += x.sum_phi_1; // + sum_r += x.sum_r; + + // square + sum_eta2 += x.sum_eta2; + sum_phi2_0 += x.sum_phi2_0; + sum_phi2_1 += x.sum_phi2_1; + sum_r2 += x.sum_r2; + + // off diagonal + sum_eta_r += x.sum_eta_r ; + sum_r_phi_0 += x.sum_r_phi_0 ; + sum_r_phi_1 += x.sum_r_phi_1 ; + sum_eta_phi_0 += x.sum_eta_phi_0; + sum_eta_phi_1 += x.sum_eta_phi_1; + +} + + +// -------------- CLUSTER SHAPES --------------- +void ClusterShapes::Init(float e ,float eta, float phi, float r){ + if (e<=0 ) return; + sum_e = e; + sum_e2 = e*e; + sum_logE = std::log(e); + + float w = e; + + n=1; + + sum_w = w; + + sum_phi_0 = w *( phi ); + sum_phi_1 = w* (phi + M_PI); + sum_r = w * r; + sum_eta = w * eta; + + //-- + sum_r2 += w * (r*r); + sum_eta2 += w * (eta*eta); + sum_phi2_0 += w * (phi*phi); + sum_phi2_1 += w * (phi+M_PI)*(phi+M_PI); + + // -- off diagonal + sum_eta_r += w * (r*eta); + sum_r_phi_0 += w* (r *phi); + sum_r_phi_1 += w* r *(phi + M_PI); + sum_eta_phi_0 += w* (eta *phi); + sum_eta_phi_1 += w* eta * (phi+M_PI); + +} +// ------ +float ClusterShapes::Eta()const { return sum_eta/sum_w;} +float ClusterShapes::R() const { return sum_r/sum_w;} + +float ClusterShapes::SigmaEtaEta()const {return sum_eta2/sum_w - Eta()*Eta();} + +float ClusterShapes::SigmaRR()const { return sum_r2/sum_w - R() *R();} + + +float ClusterShapes::SigmaPhiPhi()const { + float phi_0 = (sum_phi_0 / sum_w); + float phi_1 = (sum_phi_1 / sum_w); + float spp_0 = sum_phi2_0 / sum_w - phi_0*phi_0; + float spp_1 = sum_phi2_1 / sum_w - phi_1*phi_1; + + if (spp_0 < spp_1 ) + { + float phi = phi_0; + isPhi0=true; + while (phi < - M_PI) phi += 2*M_PI; + while (phi > M_PI) phi -= 2*M_PI; + return spp_0; + } + else + { + float phi = phi_1 ; + isPhi0=false; + while (phi < - M_PI) phi += 2*M_PI; + while (phi > M_PI) phi -= 2*M_PI; + return spp_1; + } +} + +float ClusterShapes::Phi()const { + SigmaPhiPhi(); //update phi + if (isPhi0) return (sum_phi_0 / sum_w); + else return (sum_phi_1 / sum_w); +} + + +// off - diagonal +float ClusterShapes::SigmaEtaR() const { return -(sum_eta_r / sum_w - Eta() *R()) ;} + +float ClusterShapes::SigmaEtaPhi()const { + SigmaPhiPhi() ; // decide which phi use, update phi + + if (isPhi0) + return -(sum_eta_phi_0 /sum_w - Eta()*(sum_phi_0 / sum_w)); + else + return -(sum_eta_phi_1 / sum_w - Eta()*(sum_phi_1 / sum_w)); +} + +float ClusterShapes::SigmaRPhi()const { + SigmaPhiPhi() ; // decide which phi use, update phi + if (isPhi0) + return -(sum_r_phi_0 / sum_w - R() *(sum_phi_0 / sum_w)); + else + return -(sum_r_phi_1 / sum_w - R() * (sum_phi_1 / sum_w)); +} + +// ----------------------------------- +// Local Variables: +// mode:c++ +// indent-tabs-mode:nil +// tab-width:4 +// c-basic-offset:4 +// End: +// vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 diff --git a/DataFormats/L1THGCal/src/HGCalCluster.cc b/DataFormats/L1THGCal/src/HGCalCluster.cc index 2b90c36fe3a93..610f61752d3bd 100644 --- a/DataFormats/L1THGCal/src/HGCalCluster.cc +++ b/DataFormats/L1THGCal/src/HGCalCluster.cc @@ -31,3 +31,4 @@ bool HGCalCluster::operator<(const HGCalCluster& cl) const } return res; } + diff --git a/DataFormats/L1THGCal/src/classes.h b/DataFormats/L1THGCal/src/classes.h index 594796b421334..1a895906fd1d5 100644 --- a/DataFormats/L1THGCal/src/classes.h +++ b/DataFormats/L1THGCal/src/classes.h @@ -5,6 +5,8 @@ #include "DataFormats/L1THGCal/interface/HGCalTower.h" #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h" +#include "DataFormats/L1THGCal/interface/ClusterShapes.h" + namespace DataFormats { namespace L1THGCal { l1t::HGCFETriggerDigi hgcfetd; @@ -27,5 +29,7 @@ namespace DataFormats { edm::Wrapper w_hgcalTowerBxColl; edm::Wrapper w_hgcalClusterBxColl; edm::Wrapper w_hgcalTriggerCellBxColl; + + l1t::ClusterShapes clusterShapes; } } diff --git a/DataFormats/L1THGCal/src/classes_def.xml b/DataFormats/L1THGCal/src/classes_def.xml index 56a186957b6f6..f9019f19397bd 100644 --- a/DataFormats/L1THGCal/src/classes_def.xml +++ b/DataFormats/L1THGCal/src/classes_def.xml @@ -29,4 +29,6 @@ + + diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index ef2b3876a30ea..9bbc75977b97b 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -6,6 +6,7 @@ // HGCalClusters and detId #include "DataFormats/L1THGCal/interface/HGCalCluster.h" +#include "DataFormats/L1THGCal/interface/ClusterShapes.h" #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" // PF Cluster definition @@ -23,6 +24,7 @@ // Energy calibration #include "L1Trigger/L1THGCal/interface/be_algorithms/HGCalTriggerCellCalibration.h" + // Print out something before crashing or throwing exceptions #include #include @@ -60,7 +62,15 @@ namespace HGCalTriggerBackend{ // token edm::EDGetTokenT< std::vector > sim_token_; + // Digis + edm::EDGetToken inputee_, inputfh_, inputbh_; + // add to cluster shapes, once per trigger cell + void addToClusterShapes(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float energy,float eta, float phi, float r=0.0){ + auto pair = cluster_container.emplace(pid, std::pair(0,l1t::HGCalCluster() ) ) ; + auto iterator = pair.first; + iterator -> second . second . shapes.Add( energy,eta,phi,r); // last is r, for 3d clusters + } // add to cluster void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float energy,float eta, float phi) { @@ -89,6 +99,7 @@ namespace HGCalTriggerBackend{ pp4.SetM ( 0 ) ; p4 += pp4; iterator -> second . second . setP4(p4); + //iterator -> second . second . shapes.Add( energy,eta,phi,r); // last is r, for 3d clusters return ; } @@ -118,6 +129,10 @@ namespace HGCalTriggerBackend{ // pf clusters cannot be safely cast to SimCluster //sim_token_ = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","DIGI")); sim_token_ = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","")); + // -- inputee_ = cc.consumes(conf.getParameter("eeDigis")); + //inputee_ = cc.consumes(edm::InputTag("mix","HGCDigisEE","")); + //inputfh_ = cc.consumes(edm::InputTag("mix","HGCDigisHEfront","")); + // -- inputbh_ = cc.consumes(conf.getParameter("bhDigis")); } // setProduces @@ -163,6 +178,26 @@ namespace HGCalTriggerBackend{ #ifdef HGCAL_DEBUG cout<<"[HGCalTriggerSimCluster]::[run] Start"< energy + /* + edm::Handle ee_digis_h; + edm::Handle fh_digis_h, bh_digis_h; + + evt.getByToken(inputee_,ee_digis_h); + evt.getByToken(inputfh_,fh_digis_h); + //evt.getByToken(inputbh_,bh_digis_h); + + if (not ee_digis_h.isValid()){ + throw cms::Exception("ContentError")<<"[HGCalTriggerSimCluster]::[run]::[ERROR] EE Digis from HGC not available"; + } + if (not fh_digis_h.isValid()){ + throw cms::Exception("ContentError")<<"[HGCalTriggerSimCluster]::[run]::[ERROR] FH Digis from HGC not available"; + } +#ifdef HCAL_DEBUG + std::cout <<"[HGCalTriggerSimCluster]::[run]::[INFO] BH digis not loaded"< > cluster_container;// PID-> bx,cluster evt.getByToken(sim_token_,sim_handle_); @@ -228,6 +263,7 @@ namespace HGCalTriggerBackend{ auto calibratedDigiEnergy=calibratedtriggercell.p4().E(); double eta=triggercell.p4().Eta(); double phi=triggercell.p4().Phi(); + double z = triggercell.position().z(); // may be useful for cluster shapes //2.B get the HGCAL-base-cell associated to it / geometry //const auto& tc=geom->triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& //for(const auto& cell : tc.components() ) // HGcell -- unsigned @@ -235,6 +271,8 @@ namespace HGCalTriggerBackend{ // normalization loop -- this loop is fast ~4 elem double norm=0.0; + map energy_for_cluster_shapes; + for(const auto& cell : geometry_->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned { HGCalDetId cellId(cell); @@ -248,6 +286,7 @@ namespace HGCalTriggerBackend{ const auto & pid= p.first; const auto & fraction=p.second; norm += fraction; + energy_for_cluster_shapes[pid] += calibratedDigiEnergy *fraction; // norm will be done later } } // @@ -272,6 +311,12 @@ namespace HGCalTriggerBackend{ addToCluster(cluster_container, pid, 0,energy,eta,phi ) ; // how do I get eta, phi w/o the hgcal geometry? } } + for(const auto& iterator :energy_for_cluster_shapes) + { + double energy = iterator.second / norm; + unsigned pid = iterator.first; + addToClusterShapes(cluster_container, pid, 0,energy,eta,phi,z ) ;// only one for trigger cell + } } //end of for-loop } From 5c9810651b50e6281a22dd7259c95456b66dbeea Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Thu, 2 Feb 2017 09:37:17 +0100 Subject: [PATCH 19/21] weighted average using simhit energy --- .../plugins/be_algorithms/HGCalSimCluster.cc | 74 ++++++++++++++----- 1 file changed, 56 insertions(+), 18 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index 9bbc75977b97b..f3fe67ee6bdf1 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -24,6 +24,13 @@ // Energy calibration #include "L1Trigger/L1THGCal/interface/be_algorithms/HGCalTriggerCellCalibration.h" +#include "DataFormats/HGCDigi/interface/HGCDigiCollections.h" + +#include "SimDataFormats/CaloHit/interface/PCaloHit.h" +#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" +#include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h" + + // Print out something before crashing or throwing exceptions #include @@ -31,7 +38,7 @@ /* Original Author: Andrea Carlo Marini - * Original Dat: 23 Aug 2016 + * Original Date: 23 Aug 2016 * * This backend algorithm is supposed to use the sim cluster information to handle an * optimal clustering algorithm, benchmark the performances of the enconding ... @@ -129,10 +136,9 @@ namespace HGCalTriggerBackend{ // pf clusters cannot be safely cast to SimCluster //sim_token_ = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","DIGI")); sim_token_ = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","")); - // -- inputee_ = cc.consumes(conf.getParameter("eeDigis")); - //inputee_ = cc.consumes(edm::InputTag("mix","HGCDigisEE","")); - //inputfh_ = cc.consumes(edm::InputTag("mix","HGCDigisHEfront","")); - // -- inputbh_ = cc.consumes(conf.getParameter("bhDigis")); + inputee_ = cc.consumes(edm::InputTag("g4SimHits","HGCHitsEE","")); + inputfh_ = cc.consumes(edm::InputTag("g4SimHits","HGCHitsHEfront","")); + // inputbh_ = cc.consumes(conf.getParameter("g4SimHits:HGCHitsHEback")); } // setProduces @@ -179,24 +185,45 @@ namespace HGCalTriggerBackend{ cout<<"[HGCalTriggerSimCluster]::[run] Start"< energy - /* - edm::Handle ee_digis_h; - edm::Handle fh_digis_h, bh_digis_h; + + edm::Handle ee_simhits_h; + evt.getByToken(inputee_,ee_simhits_h); + + edm::Handle fh_simhits_h; + evt.getByToken(inputfh_,fh_simhits_h); + + edm::Handle bh_simhits_h; + //evt.getByToken(inputbh_,bh_simhits_h); + - evt.getByToken(inputee_,ee_digis_h); - evt.getByToken(inputfh_,fh_digis_h); - //evt.getByToken(inputbh_,bh_digis_h); - if (not ee_digis_h.isValid()){ + if (not ee_simhits_h.isValid()){ throw cms::Exception("ContentError")<<"[HGCalTriggerSimCluster]::[run]::[ERROR] EE Digis from HGC not available"; } - if (not fh_digis_h.isValid()){ + if (not fh_simhits_h.isValid()){ throw cms::Exception("ContentError")<<"[HGCalTriggerSimCluster]::[run]::[ERROR] FH Digis from HGC not available"; } #ifdef HCAL_DEBUG - std::cout <<"[HGCalTriggerSimCluster]::[run]::[INFO] BH digis not loaded"< hgc_simhit_energy; + + if (ee_simhits_h.isValid()) + { + for (const auto& simhit : *ee_simhits_h) + hgc_simhit_energy[simhit.id()] = simhit.energy(); + } + if (fh_simhits_h.isValid()) + { + for (const auto& simhit : *fh_simhits_h) + hgc_simhit_energy[simhit.id()] = simhit.energy(); + } + if (bh_simhits_h.isValid()) + { + for (const auto& simhit : *bh_simhits_h) + hgc_simhit_energy[simhit.id()] = simhit.energy(); + } + //1. construct a cluster container that hosts the cluster per truth-particle std::unordered_map > cluster_container;// PID-> bx,cluster @@ -276,6 +303,12 @@ namespace HGCalTriggerBackend{ for(const auto& cell : geometry_->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned { HGCalDetId cellId(cell); + + //2.C0 find energy of the hgc cell + double hgc_energy=1.0; // 1 -> average if not found / bh + const auto &it = hgc_simhit_energy.find(cell); + if (it != hgc_simhit_energy.end()) { hgc_energy = it->second; } + //2.C get the particleId and energy fractions const auto & iterator= simclusters.find(cellId); if (iterator == simclusters.end() ) continue; @@ -285,14 +318,18 @@ namespace HGCalTriggerBackend{ { const auto & pid= p.first; const auto & fraction=p.second; - norm += fraction; - energy_for_cluster_shapes[pid] += calibratedDigiEnergy *fraction; // norm will be done later + norm += fraction * hgc_energy; + energy_for_cluster_shapes[pid] += calibratedDigiEnergy *fraction *hgc_energy; // norm will be done later, with the above function } } // for(const auto& cell : geometry_->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned { HGCalDetId cellId(cell); + // + double hgc_energy=1.0; // 1 -> average if not found / bh + const auto &it = hgc_simhit_energy.find(cell); + if (it != hgc_simhit_energy.end()) { hgc_energy = it->second; } //2.C get the particleId and energy fractions const auto & iterator= simclusters.find(cellId); if (iterator == simclusters.end() ) continue; @@ -303,7 +340,8 @@ namespace HGCalTriggerBackend{ const auto & pid= p.first; const auto & fraction=p.second; //auto energy = fraction*digiEnergy/ncells; - auto energy = fraction * calibratedDigiEnergy/norm; + //auto energy = fraction * calibratedDigiEnergy/norm; + auto energy = fraction * hgc_energy* calibratedDigiEnergy/norm; //2.D add to the corresponding cluster //void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float & energy,float & eta, float &phi) From e20b8e20797eafa8f6a0af69c51fa817cc986734 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Thu, 9 Feb 2017 10:52:57 +0100 Subject: [PATCH 20/21] bug fixes --- L1Trigger/L1THGCal/plugins/BuildFile.xml | 1 + .../plugins/be_algorithms/HGCalSimCluster.cc | 81 ++++++++++++++-- .../L1THGCal/test/testTriggerSimCluster.py | 97 +++++++++++++++++-- 3 files changed, 163 insertions(+), 16 deletions(-) diff --git a/L1Trigger/L1THGCal/plugins/BuildFile.xml b/L1Trigger/L1THGCal/plugins/BuildFile.xml index e834524c0a915..10f090cd2222e 100644 --- a/L1Trigger/L1THGCal/plugins/BuildFile.xml +++ b/L1Trigger/L1THGCal/plugins/BuildFile.xml @@ -21,6 +21,7 @@ + diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc index f3fe67ee6bdf1..10155d3edcdc9 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCalSimCluster.cc @@ -100,11 +100,16 @@ namespace HGCalTriggerBackend{ p4.SetM ( iterator -> second . second . mass() ) ; math::PtEtaPhiMLorentzVectorD pp4; float t = std::exp (- eta); - pp4.SetPt ( energy * (1-t*t)/(1+t*t) ) ; + //pp4.SetPt ( energy * (1-t*t)/(1+t*t) ) ; + pp4.SetPt ( energy * (2*t)/(1+t*t) ) ; pp4.SetEta( eta ) ; pp4.SetPhi( phi ) ; pp4.SetM ( 0 ) ; + //cout<<" - Adding to Cluster "< second . second . setP4(p4); //iterator -> second . second . shapes.Add( energy,eta,phi,r); // last is r, for 3d clusters return ; @@ -208,21 +213,54 @@ namespace HGCalTriggerBackend{ #endif std::unordered_map hgc_simhit_energy; + edm::ESHandle topo_ee, topo_fh; + es.get().get("HGCalEESensitive",topo_ee); + es.get().get("HGCalHESiliconSensitive",topo_fh); + + //cout <<"--- SIM HITS ----" < recoLayerCell = topo_ee->dddConstants().simToReco(cell,layer,sec,topo_ee->detectorType()); + cell = recoLayerCell.first; + layer = recoLayerCell.second; + if (layer<0 || cell<0) { + continue; + } + unsigned recoCell = HGCalDetId(mysubdet,zp,layer,subsec,sec,cell); + + //cout <<"* Adding ee sim hit "< recoLayerCell = topo_fh->dddConstants().simToReco(cell,layer,sec,topo_fh->detectorType()); + cell = recoLayerCell.first; + layer = recoLayerCell.second; + if (layer<0 || cell<0) { + continue; + } + unsigned recoCell = HGCalDetId(mysubdet,zp,layer,subsec,sec,cell); + //cout <<"* Adding fh sim hit "<().get(HGCalEESensitive_, hgceeTopoHandle_); es.get().get(HGCalHESiliconSensitive_, hgchefTopoHandle_); + // GET HGCAl Geometry --DEBUG + //edm::ESHandle geom_ee, geom_fh; //DEBUG + //es.get().get("HGCalEESensitive", geom_ee); + //es.get().get("HGCalHESiliconSensitive", geom_fh); // 1.5. pre-process the sim cluster to have easy accessible information #ifdef HGCAL_DEBUG cout<<"[HGCalTriggerSimCluster]::[run] processing sim clusters"< [ (pid, fraction), ... + //cout <<"[HGCalTriggerSimCluster]::[run] ---- SIM HITS"< > > simclusters; for (auto& cluster : *sim_handle_) { auto pid= cluster.particleId(); // not pdgId + //cout <<"* Find PID="<getPosition(p.first); + //cout <<" - HIT: "<< cellpos.eta()<<":"<< cellpos.phi()<( pid,p.second) ) ; } } + cout <<"--------------------------"< average if not found / bh + double hgc_energy=1.0e-10; // 1 -> average if not found / bh const auto &it = hgc_simhit_energy.find(cell); if (it != hgc_simhit_energy.end()) { hgc_energy = it->second; } + //else { cout <<"SIM HIT ENERGY NOT FOUND: id="<< cell< average if not found / bh + double hgc_energy=1.0e-10; // 1 -> average if not found / bh const auto &it = hgc_simhit_energy.find(cell); if (it != hgc_simhit_energy.end()) { hgc_energy = it->second; } + /* + cout <<"considering cellId in layer: "< >& cluster_container, uint64_t pid,int pdgid,float & energy,float & eta, float &phi) diff --git a/L1Trigger/L1THGCal/test/testTriggerSimCluster.py b/L1Trigger/L1THGCal/test/testTriggerSimCluster.py index 8883d7868ccd2..b980554c32230 100644 --- a/L1Trigger/L1THGCal/test/testTriggerSimCluster.py +++ b/L1Trigger/L1THGCal/test/testTriggerSimCluster.py @@ -15,17 +15,45 @@ from DataFormats.FWLite import Handle, Events # open file (you can use 'edmFileUtil -d /store/whatever.root' to get the physical file name) -events = Events("file:step2.root") +events = Events("file:test.root") # BXVector "hgcalTriggerPrimitiveDigiProducer" "HGCalTriggerSimClusterBestChoice" "DIGI" #clusters, clusterLabel = Handle("BXVector"), "hgcalTriggerPrimitiveDigiProducer:HGCalTriggerSimClusterBestChoice" clusters, clusterLabel = Handle("l1t::HGCalClusterBxCollection"), "hgcalTriggerPrimitiveDigiProducer:HGCalTriggerSimClusterBestChoice" +genParticles, genParticlesLabel = Handle("vector"),"genParticles" -verbose=True +verbose=False hEta=ROOT.TH1D("eta","eta",1000,-20,20) hPhi=ROOT.TH1D("phi","phi",1000,-3.1416,3.1416) hPt = ROOT.TH1D("pt","pt",1000,0,1000) hE = ROOT.TH1D("E","E",1000,0,1000) +## cluster shapes +hN = ROOT.TH1D("N","N",50,0,1000) +hSEE = ROOT.TH1D("SEE","SEE",1000,0,0.02) +hSPP = ROOT.TH1D("SPP","SPP",1000,0,0.02) +hSRR = ROOT.TH1D("SRR","SRR",1000,0,100.) +## +hLE = ROOT.TH1D("LogEoverE","Log(e)/E",1000,-0.5,1.) +hED = ROOT.TH1D("ED","ED",1000,0.,1.) +hSEP = ROOT.TH1D("SEP","SEP",1000,-0.02,0.02) + +hMatch={ + "deta":ROOT.TH1D("eta_m","eta",1000,-20,20), + "dphi":ROOT.TH1D("phi_m","phi",1000,-3.1416,3.1416), + "pt" : ROOT.TH1D("pt_m","pt",1000,0,1000), + "ptsum" : ROOT.TH1D("pt_sum_m","pt",1000,0,1000), + } + +import math +def deltaPhi(phi1,phi2): + pi=3.14159265359 + dphi=phi1-phi2 + while dphi > pi: dphi -= 2*pi + while dphi < -pi: dphi += 2*pi + return dphi + +def deltaR(eta1,phi1,eta2,phi2): + return math.sqrt( (eta1-eta2)**2 + deltaPhi(phi1,phi2)**2) for iev,event in enumerate(events): if verbose: @@ -33,30 +61,56 @@ #if iev > 10: break event.getByLabel(clusterLabel, clusters) + event.getByLabel(genParticlesLabel,genParticles) ## check that handle is valid if not clusters.isValid(): print "Cluster handle is invalid" + if not genParticles.isValid(): print "GP handle not valid" ## -- if verbose: print "-> BX=0" if verbose: print "-> clusters=",clusters,"\n","product=",clusters.product() - print "what's saved in clusters?" - for bx in clusters.product(): - print "bx=",bx + #print "what's saved in clusters?" + #for bx in clusters.product(): + # print "bx=",bx + print "--------------------------------" + print "GP SIZE=",genParticles.product().size() + for gp in genParticles.product(): + if gp.status() == 1: + print "-> Found GP",gp.pt(),gp.eta(),gp.phi(),gp.pdgId() + print "--------------------------------" ## get 0 bunch crossing vector + ptsum = [ROOT.TLorentzVector(),ROOT.TLorentzVector() ] for bx,vector in enumerate(clusters.product()): if verbose: print "-> 0 pos",clusters.product()[0] if vector == None: print " cluster product is none" continue - print " pt=",vector.pt(),"eta=",vector.eta(),"phi=",vector.phi() + #print " pt=",vector.pt(),"eta=",vector.eta(),"phi=",vector.phi() hEta.Fill(vector.eta() ) if vector.eta() <8: hE.Fill(vector.energy() ) hPt.Fill(vector.pt() ) hPhi.Fill(vector.phi() ) + ## cluster shapes + hN.Fill(vector.shapes.N()) + hSEE.Fill(vector.shapes.SigmaEtaEta()) + hSPP.Fill(vector.shapes.SigmaPhiPhi()) + hSRR.Fill(vector.shapes.SigmaRR()) + hLE.Fill(vector.shapes.LogEoverE() ) + hED.Fill(vector.shapes.eD()) + hSEP.Fill(vector.shapes.SigmaEtaPhi() ) + for idx,gp in enumerate(genParticles.product()): + if gp.status()==1 and deltaR(vector.eta(),vector.phi(),gp.eta(),gp.phi())<0.1: + hMatch["deta"].Fill(vector.eta()-gp.eta()) + hMatch["dphi"].Fill(deltaPhi(vector.phi(),gp.phi())) + hMatch["pt"].Fill(vector.pt()) + ptsum[idx] += ROOT.TLorentzVector(vector.px(),vector.py(),vector.pz(),vector.energy()) + hMatch["ptsum"].Fill( ptsum[0].Pt() ) + hMatch["ptsum"].Fill( ptsum[1].Pt() ) + c=ROOT.TCanvas() @@ -69,5 +123,36 @@ hPhi.Draw("HIST") c.cd(4) hE.Draw("HIST") + +c2=ROOT.TCanvas("c2","c2") +c2.Divide(2,2) +c2.cd(1) +hN.Draw("HIST") +c2.cd(2) +hSEE.Draw("HIST") +c2.cd(3) +hSPP.Draw("HIST") +c2.cd(4) +hSRR.Draw("HIST") + +#c3=ROOT.TCanvas("c3","c3") +#c3.Divide(2,2) +#c3.cd(1) +#hLE.Draw("HIST") +#c3.cd(2) +#hED.Draw("HIST") +#c3.cd(3) +#hSEP.Draw("HIST") + +c4=ROOT.TCanvas("c4","c4") +c4.Divide(2,2) +c4.cd(1) +hMatch["pt"].Draw("HIST") +c4.cd(2) +hMatch["deta"].Draw("HIST") +c4.cd(3) +hMatch["dphi"].Draw("HIST") +c4.cd(4) +hMatch["ptsum"].Draw("HIST") raw_input("ok?") From ea0a3cdc6b27ec9837d6b36ae0d65730d69cfaab Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Thu, 9 Feb 2017 18:02:55 +0100 Subject: [PATCH 21/21] cleaning couts and adding configuartion --- .../plugins/HGCalTriggerDigiProducer.cc | 24 --- .../plugins/be_algorithms/HGCalSimCluster.cc | 145 +++--------------- ...hgcalTriggerPrimitivesSimClustering_cff.py | 17 ++ 3 files changed, 42 insertions(+), 144 deletions(-) create mode 100644 L1Trigger/L1THGCal/python/hgcalTriggerPrimitivesSimClustering_cff.py diff --git a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc index 33668ca192d5d..9eb54b1e83adf 100644 --- a/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc +++ b/L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc @@ -15,12 +15,6 @@ #include #include -#define HGCAL_DEBUG - -#ifdef HGCAL_DEBUG - #include - using namespace std; -#endif class HGCalTriggerDigiProducer : public edm::EDProducer { public: @@ -49,9 +43,6 @@ HGCalTriggerDigiProducer(const edm::ParameterSet& conf): //inputbh_(consumes(conf.getParameter("bhDigis"))), backEndProcessor_(new HGCalTriggerBackendProcessor(conf.getParameterSet("BEConfiguration"),consumesCollector()) ) { - #ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerDigiProducer]::["<< __FUNCTION__ <<"] constructor: file " <<__FILE__<<":"<<__LINE__<("CodecName"); @@ -65,30 +56,18 @@ HGCalTriggerDigiProducer(const edm::ParameterSet& conf): // register backend processor products backEndProcessor_->setProduces(*this); - #ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerDigiProducer]::["<< __FUNCTION__ <<"] end: file " <<__FILE__<<":"<<__LINE__<().get(triggerGeometry_); codec_->setGeometry(triggerGeometry_.product()); backEndProcessor_->setGeometry(triggerGeometry_.product()); - #ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerDigiProducer]::["<< __FUNCTION__ <<"] end" < fe_output( new l1t::HGCFETriggerDigiCollection ); @@ -157,7 +136,4 @@ void HGCalTriggerDigiProducer::produce(edm::Event& e, const edm::EventSetup& es) backEndProcessor_->run(fe_digis_coll,es,e); backEndProcessor_->putInEvent(e); //backEndProcessor_->reset(); - #ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerDigiProducer]::["<< __FUNCTION__ <<"] end. file " <<__FILE__<<":"<<__LINE__<(0,l1t::HGCalCluster() ) ) ; auto iterator = pair.first; - //auto iterator = cluster_container.find (pid); - //if (iterator == cluster_container.end()) - // { - // //create an empty cluster - // cluster_container[pid] = std::pair(0,l1t::HGCalCluster()); - // iterator = cluster_container.find (pid); - // iterator -> second . second . setPdgId(pdgid); - // } - // p4 += p4' math::PtEtaPhiMLorentzVectorD p4; p4.SetPt ( iterator -> second . second . pt() ) ; p4.SetEta( iterator -> second . second . eta() ) ; @@ -100,16 +92,11 @@ namespace HGCalTriggerBackend{ p4.SetM ( iterator -> second . second . mass() ) ; math::PtEtaPhiMLorentzVectorD pp4; float t = std::exp (- eta); - //pp4.SetPt ( energy * (1-t*t)/(1+t*t) ) ; pp4.SetPt ( energy * (2*t)/(1+t*t) ) ; pp4.SetEta( eta ) ; pp4.SetPhi( phi ) ; pp4.SetM ( 0 ) ; - //cout<<" - Adding to Cluster "< second . second . setP4(p4); //iterator -> second . second . shapes.Add( energy,eta,phi,r); // last is r, for 3d clusters return ; @@ -136,48 +123,29 @@ namespace HGCalTriggerBackend{ HGCalHESiliconSensitive_(conf.getParameter("HGCalHESiliconSensitive_tag")), calibration_(conf.getParameterSet("calib_parameters")) { - // I need to consumes the PF Cluster Collection with the sim clustering, TODO: make it configurable (?) - // vector "mix" "MergedCaloTruth" "HLT/DIGI" - // pf clusters cannot be safely cast to SimCluster - //sim_token_ = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","DIGI")); - sim_token_ = cc.consumes< std::vector< SimCluster > >(edm::InputTag("mix","MergedCaloTruth","")); - inputee_ = cc.consumes(edm::InputTag("g4SimHits","HGCHitsEE","")); - inputfh_ = cc.consumes(edm::InputTag("g4SimHits","HGCHitsHEfront","")); + sim_token_ = cc.consumes< std::vector< SimCluster > >(conf.getParameter("simcollection")); + inputee_ = cc.consumes( conf.getParameter("simhitsee")); + inputfh_ = cc.consumes( conf.getParameter("simhitsfh")); // inputbh_ = cc.consumes(conf.getParameter("g4SimHits:HGCHitsHEback")); + edm::LogWarning("HGCalTriggerSimCluster") <<"WARNING: BH simhits not loaded"; } // setProduces virtual void setProduces(edm::EDProducer& prod) const override final { -#ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] start: Setting producer to produce:"<(name()); } // putInEvent virtual void putInEvent(edm::Event& evt) override final { -#ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] start"< energy edm::Handle ee_simhits_h; @@ -200,7 +165,6 @@ namespace HGCalTriggerBackend{ edm::Handle bh_simhits_h; //evt.getByToken(inputbh_,bh_simhits_h); - if (not ee_simhits_h.isValid()){ throw cms::Exception("ContentError")<<"[HGCalTriggerSimCluster]::[run]::[ERROR] EE Digis from HGC not available"; @@ -208,16 +172,13 @@ namespace HGCalTriggerBackend{ if (not fh_simhits_h.isValid()){ throw cms::Exception("ContentError")<<"[HGCalTriggerSimCluster]::[run]::[ERROR] FH Digis from HGC not available"; } -#ifdef HCAL_DEBUG - std::cout <<"[HGCalTriggerSimCluster]::[run]::[INFO] BH simhits not loaded"< hgc_simhit_energy; edm::ESHandle topo_ee, topo_fh; es.get().get("HGCalEESensitive",topo_ee); es.get().get("HGCalHESiliconSensitive",topo_fh); - //cout <<"--- SIM HITS ----" < > cluster_container;// PID-> bx,cluster evt.getByToken(sim_token_,sim_handle_); if (not sim_handle_.isValid()){ - throw cms::Exception("ContentError")<<"[HGCalTriggerSimCluster]::[run]::[ERROR] PFCluster collection for HGC sim clustering not available"; + throw cms::Exception("ContentError")<<"[HGCalTriggerSimCluster]::[run]::[ERROR] Sim Cluster collection for HGC sim clustering not available"; } // calibration es.get().get(HGCalEESensitive_, hgceeTopoHandle_); es.get().get(HGCalHESiliconSensitive_, hgchefTopoHandle_); - // GET HGCAl Geometry --DEBUG - //edm::ESHandle geom_ee, geom_fh; //DEBUG - //es.get().get("HGCalEESensitive", geom_ee); - //es.get().get("HGCalHESiliconSensitive", geom_fh); // 1.5. pre-process the sim cluster to have easy accessible information -#ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerSimCluster]::[run] processing sim clusters"< [ (pid, fraction), ... - //cout <<"[HGCalTriggerSimCluster]::[run] ---- SIM HITS"< > > simclusters; for (auto& cluster : *sim_handle_) { auto pid= cluster.particleId(); // not pdgId - //cout <<"* Find PID="<getPosition(p.first); - //cout <<" - HIT: "<< cellpos.eta()<<":"<< cellpos.phi()<( pid,p.second) ) ; } } - cout <<"--------------------------"<triggerCells()[ tcellId() ] ;//HGCalTriggerGeometry::TriggerCell& - //for(const auto& cell : tc.components() ) // HGcell -- unsigned - unsigned ncells = geometry_->getCellsFromTriggerCell( tcellId()).size(); - // normalization loop -- this loop is fast ~4 elem + // normalization loop double norm=0.0; map energy_for_cluster_shapes; @@ -353,47 +294,34 @@ namespace HGCalTriggerBackend{ { HGCalDetId cellId(cell); - //2.C0 find energy of the hgc cell - double hgc_energy=1.0e-10; // 1 -> average if not found / bh + //2.C0 find energy of the hgc cell. default is very small value + double hgc_energy=1.0e-10; //average if not found / bh const auto &it = hgc_simhit_energy.find(cell); if (it != hgc_simhit_energy.end()) { hgc_energy = it->second; } - //else { cout <<"SIM HIT ENERGY NOT FOUND: id="<< cell<second; for ( const auto& p: particles ) { const auto & pid= p.first; const auto & fraction=p.second; norm += fraction * hgc_energy; - energy_for_cluster_shapes[pid] += calibratedDigiEnergy *fraction *hgc_energy; // norm will be done later, with the above function + energy_for_cluster_shapes[pid] += calibratedDigiEnergy *fraction *hgc_energy; // norm will be done later, with the above norm } } - // + + // second loop counting the energy for(const auto& cell : geometry_->getCellsFromTriggerCell( tcellId()) ) // HGCcell -- unsigned { HGCalDetId cellId(cell); - // double hgc_energy=1.0e-10; // 1 -> average if not found / bh const auto &it = hgc_simhit_energy.find(cell); if (it != hgc_simhit_energy.end()) { hgc_energy = it->second; } - /* - cout <<"considering cellId in layer: "<second; for ( const auto& p: particles ) { @@ -405,11 +333,12 @@ namespace HGCalTriggerBackend{ //auto energy = fraction * hgc_energy* hgc_energy/norm; //2.D add to the corresponding cluster - //void addToCluster(std::unordered_map >& cluster_container, uint64_t pid,int pdgid,float & energy,float & eta, float &phi) - //addToCluster(cluster_container, pid, 0 energy,ETA/PHI? ) ; - addToCluster(cluster_container, pid, 0,energy,eta,phi ) ; // how do I get eta, phi w/o the hgcal geometry? + //eta/phi are the position of the trigger cell, to consider degradation + addToCluster(cluster_container, pid, 0,energy,eta,phi ) ; } } + + // third loop, cluster shapes to ensure the correct counts for(const auto& iterator :energy_for_cluster_shapes) { double energy = iterator.second / norm; @@ -419,35 +348,13 @@ namespace HGCalTriggerBackend{ } //end of for-loop } -#ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerSimCluster]::[run] Push clusters in cluster products"<getTriggerCellPosition( tcellId()); - - // construct it from *both* physical and integer values - //l1t::HGCalCluster cluster( reco::LeafCandidate::LorentzVector(), - // clusterEnergyHw, clusterEtaHw, 0); - // for (auto& p : cluster_container) { //std::unordered_map > - #ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerSimCluster]::[run] Cluster: pid="<< p.first<push_back(p.second.first,p.second.second); // bx,cluster } -#ifdef HGCAL_DEBUG - cout<<"[HGCalTriggerSimCluster]::["<<__FUNCTION__<<"] cluster product (!=NULL) "< HGCalTriggerSimClusterBestChoice; -//DEFINE_EDM_PLUGIN(HGCalTriggerBackendAlgorithmFactory, HGCalTriggerSimClusterBestChoice,"HGCalTriggerSimClusterBestChoice"); +// typedef HGCalTriggerBackend::HGCalTriggerSimCluster HGCalTriggerSimClusterBestChoice; DEFINE_EDM_PLUGIN(HGCalTriggerBackendAlgorithmFactory, HGCalTriggerSimClusterBestChoice,"HGCalTriggerSimClusterBestChoice"); -//DEFINE_EDM_PLUGIN(HGCalTriggerBackendAlgorithmFactory, HGCalTriggerSimCluster,"HGCalTriggerSimCluster"); // Local Variables: // mode:c++ diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerPrimitivesSimClustering_cff.py b/L1Trigger/L1THGCal/python/hgcalTriggerPrimitivesSimClustering_cff.py new file mode 100644 index 0000000000000..5cb88d2e34156 --- /dev/null +++ b/L1Trigger/L1THGCal/python/hgcalTriggerPrimitivesSimClustering_cff.py @@ -0,0 +1,17 @@ +import FWCore.ParameterSet.Config as cms + +from L1Trigger.L1HGCal.hgcalTriggerPrimitives_cff import * + +cluster_algo_all = cms.PSet( AlgorithmName = cms.string('HGCalTriggerSimClusterBestChoice'), + FECodec = process.hgcalTriggerPrimitiveDigiProducer.FECodec , + HGCalEESensitive_tag = cms.string('HGCalEESensitive'), + HGCalHESiliconSensitive_tag = cms.string('HGCalHESiliconSensitive'), + + calib_parameters = process.hgcalTriggerPrimitiveDigiProducer.BEConfiguration.algorithms[0].calib_parameters, + simcollection=cms.InputTag("mix","MergedCaloTruth",""), + simhitsee = cms.InputTag("g4SimHits","HGCHitsEE",""), + simhitsfh = cms.InputTag("g4SimHits","HGCHitsHEfront",""), + ) + + +#process.hgcl1tpg_step = cms.Path(process.hgcalTriggerPrimitives)