diff --git a/L1Trigger/L1THGCal/src/fe_codecs/HGCalBestChoiceCodecImpl.cc b/L1Trigger/L1THGCal/src/fe_codecs/HGCalBestChoiceCodecImpl.cc index d623eed81d5a6..599e28ce72910 100644 --- a/L1Trigger/L1THGCal/src/fe_codecs/HGCalBestChoiceCodecImpl.cc +++ b/L1Trigger/L1THGCal/src/fe_codecs/HGCalBestChoiceCodecImpl.cc @@ -48,7 +48,7 @@ std::vector HGCalBestChoiceCodecImpl::encode(const HGCalBestChoiceCodecImp << " : Number of energy values = "<(0x1u<(0x1u<(value & (0x1<<(i+triggerCellTruncationBits_)));// remove the lowest bits (=triggerCellTruncationBits_) diff --git a/L1Trigger/L1THGCal/test/BuildFile.xml b/L1Trigger/L1THGCal/test/BuildFile.xml index 4dfe3d832a036..47b20febb84a8 100644 --- a/L1Trigger/L1THGCal/test/BuildFile.xml +++ b/L1Trigger/L1THGCal/test/BuildFile.xml @@ -9,6 +9,7 @@ + diff --git a/L1Trigger/L1THGCal/test/HGCalTriggerBestChoiceTester.cc b/L1Trigger/L1THGCal/test/HGCalTriggerBestChoiceTester.cc index 69a9da7b1ab47..e536c18ec2650 100644 --- a/L1Trigger/L1THGCal/test/HGCalTriggerBestChoiceTester.cc +++ b/L1Trigger/L1THGCal/test/HGCalTriggerBestChoiceTester.cc @@ -26,14 +26,18 @@ #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h" #include "DataFormats/ForwardDetId/interface/HGCTriggerDetId.h" +#include "SimDataFormats/CaloHit/interface/PCaloHit.h" +#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" +#include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h" + #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerFECodecBase.h" #include "L1Trigger/L1THGCal/interface/fe_codecs/HGCalBestChoiceCodecImpl.h" #include +#include #include "TH2.h" -using namespace std; class HGCalTriggerBestChoiceTester : public edm::EDAnalyzer { @@ -48,27 +52,44 @@ class HGCalTriggerBestChoiceTester : public edm::EDAnalyzer private: void checkSelectedCells(const edm::Event&, const edm::EventSetup&); void rerunBestChoiceFragments(const edm::Event&, const edm::EventSetup&); - void fillModule(const std::vector>&, - const HGCalBestChoiceDataPayload&, - const vector >& ); + void fillModule(const std::vector>&, const std::vector >&, const HGCalBestChoiceDataPayload&, const HGCalBestChoiceDataPayload&,const HGCalBestChoiceDataPayload&, const std::map &, const std::map& ); + // inputs edm::EDGetToken inputee_, inputfh_, inputbh_, inputbeall_, inputbeselect_; + bool is_Simhit_comp_; + edm::EDGetToken SimHits_inputee_, SimHits_inputfh_, SimHits_inputbh_; // std::unique_ptr triggerGeometry_; std::unique_ptr codec_; edm::Service fs_; + HGCalTriggerGeometryBase::es_info info_; + // histos - TH1F* hgcCellsPerModule_; TH1F* hgcCellData_; - TH1F* hgcCellData_linampl_; + TH1F* hgcCellData_SimHitasso_; + TH1F* hgcCellSimHits_; + TH2F* hgcCellData_vsSimHits_; + TH1F* hgcCellsPerModule_; TH1F* hgcCellModuleSum_; - TH1F* triggerCellsPerModule_; + // + TH1F* hgcCellData_linampl_; + TH2F* hgcCellData_linampl_vsSimHits_; + TH2F* hgcCellData_linampl_vsSimHits_zoom_; + TH1F* triggerCellData_noBestChoice_; + TH2F* triggerCellData_noBestChoice_vsSimHits_; + TH1F* triggerCellSimHits_noBestChoice_; + TH1F* triggerCellData_BestChoice_; + TH2F* triggerCellData_BestChoice_vsSimHits_; TH1F* triggerCellData_; + TH2F* triggerCellData_vsSimHits_; + TH1F* triggerCellsPerModule_; TH1F* triggerCellModuleSum_; + // TH2F* selectedCellsVsAllCells_ee_; TH2F* energyLossVsNCells_ee_; TH2F* selectedCellsVsAllCells_fh_; TH2F* energyLossVsNCells_fh_; + // }; @@ -79,7 +100,11 @@ HGCalTriggerBestChoiceTester::HGCalTriggerBestChoiceTester(const edm::ParameterS inputfh_(consumes(conf.getParameter("fhDigis"))), //inputbh_(consumes(conf.getParameter("bhDigis"))), inputbeall_(consumes(conf.getParameter("beClustersAll"))), - inputbeselect_(consumes(conf.getParameter("beClustersSelect"))) + inputbeselect_(consumes(conf.getParameter("beClustersSelect"))), + is_Simhit_comp_(conf.getParameter("isSimhitComp")), + SimHits_inputee_(consumes(conf.getParameter("eeSimHits"))), + SimHits_inputfh_(consumes(conf.getParameter("fhSimHits"))) + // SimHits_inputbh_(consumes(conf.getParameter("bhSimHits"))) /*****************************************************************/ { //setup geometry @@ -93,20 +118,40 @@ HGCalTriggerBestChoiceTester::HGCalTriggerBestChoiceTester(const edm::ParameterS codec_.reset( new HGCalBestChoiceCodecImpl(feCodecConfig) ); // initialize output trees - hgcCellsPerModule_ = fs_->make("hgcCellsPerModule","Number of cells per module", 64, 0., 64.); - hgcCellData_ = fs_->make("hgcCellData","Cell values", 500, 0., 500.); - // - hgcCellData_linampl_ = fs_->make("hgcCellData_linampl_","Cell linearized amplitudes values All", 1250, 0, 25000); + // HGC Cells + hgcCellData_ = fs_->make("hgcCellData","Cell values", 1000, 0., 2000.); + if (is_Simhit_comp_) { + hgcCellData_SimHitasso_ = fs_->make("hgcCellData_SimHitasso_","Cell values with an associated SimHit", 1000, 0, 2000.); + hgcCellSimHits_ = fs_->make("hgcCellSimHits_","Cell simhit energies", 500, 0, 0.16); + hgcCellData_vsSimHits_ = fs_->make("hgcCellData_vsSimHits_","Cell values vs simhit energies", 500,0,0.16,1000, 0., 2000.); + } + hgcCellsPerModule_ = fs_->make("hgcCellsPerModule","Number of cells per module", 128, 0., 128.); + hgcCellModuleSum_ = fs_->make("hgcCellModuleSum","Cell sum in modules", 1000, 0., 3000.); // - hgcCellModuleSum_ = fs_->make("hgcCellModuleSum","Cell sum in modules", 1000, 0., 1000.); - triggerCellsPerModule_ = fs_->make("TriggerCellsPerModule","Number of trigger cells per module", 64, 0., 64.); - triggerCellData_ = fs_->make("TriggerCellData","Trigger cell values", 500, 0., 500.); - triggerCellModuleSum_ = fs_->make("TriggerCellModuleSum","Trigger cell sum in modules", 1000, 0., 1000.); + hgcCellData_linampl_ = fs_->make("hgcCellData_linampl_","Cell linearized amplitudes values All", 1000, 0, 70000); + if (is_Simhit_comp_) { + hgcCellData_linampl_vsSimHits_ = fs_->make("hgcCellData_linampl_vsSimHits_","Cell linearized amplitudes vs simhit energies",500,0,0.16,1000,0,70000); + hgcCellData_linampl_vsSimHits_zoom_ = fs_->make("hgcCellData_linampl_vsSimHits_zoom_","Cell linearized amplitudes vssimhit energies, zoomed",1000,0,0.002,1000,0,1000); + } + + // HGC Trigger cells + triggerCellData_noBestChoice_ = fs_->make("triggerCellData_noBestChoice_","Trigger cell values, no best choice", 1000, 0., 70000.); + if (is_Simhit_comp_){ + triggerCellData_noBestChoice_vsSimHits_ = fs_->make("triggerCellData_noBestChoice_vsSimHits_","Trigger cell values vs simhit energies, no best choice", 500,0,0.16,1000, 0., 70000.); + triggerCellSimHits_noBestChoice_ = fs_->make("triggerCellSimHits_noBestChoice","Trigger cell simhit energies, no best choice", 500, 0, 0.16); + } + triggerCellData_BestChoice_ = fs_->make("triggerCellData_BestChoice_","Trigger cell values, best choice", 1000, 0., 70000.); + if (is_Simhit_comp_) triggerCellData_BestChoice_vsSimHits_ = fs_->make("triggerCellData_BestChoice_vsSimHits_","Trigger cell values vs simhit energies, best choice", 500,0,0.16,1000, 0., 70000.); + triggerCellData_ = fs_->make("triggerCellData","Trigger cell values", 1100, 0., 1100.); + if (is_Simhit_comp_) triggerCellData_vsSimHits_ = fs_->make("triggerCellData_vsSimHits_","Trigger cell values vs simhit energies", 500,0,0.16,1100, 0., 1100.); + triggerCellsPerModule_ = fs_->make("triggerCellsPerModule","Number of trigger cells per module", 64, 0., 64.); + triggerCellModuleSum_ = fs_->make("TriggerCellModuleSum","Trigger cell sum in modules", 1000, 0., 10000.); // selectedCellsVsAllCells_ee_ = fs_->make("selectedCellsVsAllCells_ee","Number of selected cells vs number of cell", 128, 0, 128, 128, 0., 128.); energyLossVsNCells_ee_ = fs_->make("energyLossVsNCells_ee","Relative energy loss after selection vs number of cell", 128, 0., 128., 101, 0, 1.01); selectedCellsVsAllCells_fh_ = fs_->make("selectedCellsVsAllCells_fh","Number of selected cells vs number of cell", 128, 0, 128, 128, 0., 128.); energyLossVsNCells_fh_ = fs_->make("energyLossVsNCells_fh","Relative energy loss after selection vs number of cell", 128, 0., 128., 101, 0, 1.01); + } @@ -123,18 +168,18 @@ void HGCalTriggerBestChoiceTester::beginRun(const edm::Run& /*run*/, /*****************************************************************/ { triggerGeometry_->reset(); - HGCalTriggerGeometryBase::es_info info; const std::string& ee_sd_name = triggerGeometry_->eeSDName(); const std::string& fh_sd_name = triggerGeometry_->fhSDName(); const std::string& bh_sd_name = triggerGeometry_->bhSDName(); - es.get().get(ee_sd_name,info.geom_ee); - es.get().get(fh_sd_name,info.geom_fh); - es.get().get(bh_sd_name,info.geom_bh); - es.get().get(ee_sd_name,info.topo_ee); - es.get().get(fh_sd_name,info.topo_fh); - es.get().get(bh_sd_name,info.topo_bh); - triggerGeometry_->initialize(info); -} + es.get().get(ee_sd_name,info_.geom_ee); + es.get().get(fh_sd_name,info_.geom_fh); + es.get().get(bh_sd_name,info_.geom_bh); + es.get().get(ee_sd_name,info_.topo_ee); + es.get().get(fh_sd_name,info_.topo_fh); + es.get().get(bh_sd_name,info_.topo_bh); + triggerGeometry_->initialize(info_); + + } /*****************************************************************/ void HGCalTriggerBestChoiceTester::analyze(const edm::Event& e, @@ -223,88 +268,226 @@ void HGCalTriggerBestChoiceTester::rerunBestChoiceFragments(const edm::Event& e, HGCalBestChoiceDataPayload data; - //loop on modules - for( const auto& module : triggerGeometry_->modules() ) { - HGCalDetId moduleId(module.first); - // prepare input data - std::vector> dataframes; - vector > linearized_dataframes; - - // loop over EE or FH digis and fill digis belonging to that module - if(moduleId.subdetId()==ForwardSubdetector::HGCEE) { - for(const auto& eedata : ee_digis) { - if(module.second->containsCell(eedata.id())) { - dataframes.emplace_back(eedata.id()); - for(int i=0; i simhit_energies; + if (is_Simhit_comp_) { + edm::Handle ee_simhits_h; + e.getByToken(SimHits_inputee_,ee_simhits_h); + const edm::PCaloHitContainer& ee_simhits = *ee_simhits_h; + edm::Handle fh_simhits_h; + e.getByToken(SimHits_inputfh_,fh_simhits_h); + const edm::PCaloHitContainer& fh_simhits = *fh_simhits_h; + + // simhit/digi association EE + HGCalDetId digiid, simid; + int layer=0,cell=0, sec=0, subsec=0, zp=0,subdet=0; + ForwardSubdetector mysubdet; + HGCalDetId recoDetId ; + int n_hits_asso=0; + + for(const auto& eedata : ee_digis) { + digiid= (HGCalDetId) eedata.id(); + bool is_hitasso=false; + double hit_energy=0; + + for( const auto& eesimhit : ee_simhits ) { + simid = (HGCalDetId) eesimhit.id(); + HGCalTestNumbering::unpackHexagonIndex(simid, subdet, zp, layer, sec, subsec, cell); + mysubdet = (ForwardSubdetector)(subdet); + std::pair recoLayerCell=info_.topo_ee->dddConstants().simToReco(cell,layer,sec,info_.topo_ee->detectorType()); + cell = recoLayerCell.first; + layer = recoLayerCell.second; + if (layer<0 || cell<0) { + continue; + } + recoDetId = HGCalDetId(mysubdet,zp,layer,subsec,sec,cell); + if (recoDetId==digiid) { + hit_energy+=eesimhit.energy(); + is_hitasso=true; + } } - else if(moduleId.subdetId()==ForwardSubdetector::HGCHEF) { - for(const auto& fhdata : fh_digis) { - if(module.second->containsCell(fhdata.id())) { - dataframes.emplace_back(fhdata.id()); - for(int i=0; i recoLayerCell=info_.topo_fh->dddConstants().simToReco(cell_fh,layer_fh,sec_fh,info_.topo_fh->detectorType()); + cell_fh = recoLayerCell.first; + layer_fh = recoLayerCell.second; + if (layer_fh<0 || cell_fh<0) { + continue; + } + recoDetId_fh = HGCalDetId(mysubdet_fh,zp_fh,layer_fh,subsec_fh,sec_fh,cell_fh); + if (recoDetId_fh==digiid_fh) { + hit_energy+=fhsimhit.energy(); + is_hitasso=true; + } } + if (is_hitasso){ n_hits_asso_fh++; n_hits_asso++;} + simhit_energies[digiid_fh] = hit_energy; + } + } - // Best choice encoding - data.reset(); - codec_->linearize(*(module.second), dataframes, linearized_dataframes); - codec_->triggerCellSums(*(module.second), linearized_dataframes, data); - codec_->bestChoiceSelect(data); - std::vector dataword = codec_->encode(data); - HGCalBestChoiceDataPayload datadecoded = codec_->decode(dataword); - fillModule(dataframes, datadecoded, linearized_dataframes); + //loop on modules + for( const auto& module : triggerGeometry_->modules() ) { + HGCalDetId moduleId(module.first); + // prepare input data + std::vector> dataframes; + std::vector > linearized_dataframes; + + // loop over EE or FH digis and fill digis belonging to that module + if(moduleId.subdetId()==ForwardSubdetector::HGCEE) { + for(const auto& eedata : ee_digis) { + if(module.second->containsCell(eedata.id())) { + dataframes.emplace_back(eedata.id()); + for(int i=0; icontainsCell(fhdata.id())) { + dataframes.emplace_back(fhdata.id()); + for(int i=0; i TC_simhit_energies; + if (is_Simhit_comp_) { + for(const auto& tc_c : module.second->triggerCellComponents()){ + HGCalDetId triggercellid( tc_c.first ); + uint32_t cellid = triggercellid.cell(); + TC_simhit_energies.insert( std::make_pair(cellid, 0) ); + double simenergy = simhit_energies[tc_c.second]; + TC_simhit_energies[cellid]+=simenergy; + } + } + + // Best choice encoding + data.reset(); + codec_->linearize(*(module.second), dataframes, linearized_dataframes); + codec_->triggerCellSums(*(module.second), linearized_dataframes, data); + HGCalBestChoiceDataPayload data_TCsums_woBestChoice = data; + codec_->bestChoiceSelect(data); + HGCalBestChoiceDataPayload data_TCsums_BestChoice = data; + std::vector dataword = codec_->encode(data); + HGCalBestChoiceDataPayload datadecoded = codec_->decode(dataword); + fillModule(dataframes, linearized_dataframes, data_TCsums_woBestChoice,data_TCsums_BestChoice, datadecoded, simhit_energies, TC_simhit_energies); + } //end loop on modules - + } /*****************************************************************/ -void HGCalTriggerBestChoiceTester::fillModule( const std::vector>& dataframes, - const HGCalBestChoiceDataPayload& fe_payload, - const vector >& linearized_dataframes) +void HGCalTriggerBestChoiceTester::fillModule( const std::vector>& dataframes, const std::vector >& linearized_dataframes, const HGCalBestChoiceDataPayload& fe_payload_TCsums_woBestChoice, const HGCalBestChoiceDataPayload& fe_payload_TCsums_BestChoice, const HGCalBestChoiceDataPayload& fe_payload, const std::map & simhit_energies, const std::map& TC_simhit_energies) + /*****************************************************************/ { // HGC cells part size_t nHGCDigi = 0; unsigned hgcCellModuleSum = 0; + // digis, cell based info for(const auto& frame : dataframes) - { + { uint32_t value = frame[2].data(); - if(value>0) - { - nHGCDigi++; - hgcCellModuleSum += value; - hgcCellData_->Fill(value); + nHGCDigi++; + hgcCellModuleSum += value; + hgcCellData_->Fill(value); + if (is_Simhit_comp_){ + double sim_energy= simhit_energies.at(frame.id()); + if (sim_energy >0){ + hgcCellData_SimHitasso_->Fill(value); + hgcCellSimHits_->Fill(sim_energy); + hgcCellData_vsSimHits_->Fill(sim_energy,value); + } } - } + } hgcCellsPerModule_->Fill(nHGCDigi); hgcCellModuleSum_->Fill(hgcCellModuleSum); - + + // linearized samples, cell based info for(const auto& frame : linearized_dataframes){ - hgcCellData_linampl_-> Fill(frame.second); + hgcCellData_linampl_-> Fill(frame.second); + if (is_Simhit_comp_){ + double sim_energy= simhit_energies.at(frame.first); + if (sim_energy >0){ + hgcCellData_linampl_vsSimHits_-> Fill(sim_energy,frame.second); + hgcCellData_linampl_vsSimHits_zoom_-> Fill(sim_energy,frame.second); + } + } } - + // trigger cells part + // after sum, no best choice, no encode/decode + int icell_noBestChoice = 0; + for(const auto& tc : fe_payload_TCsums_woBestChoice.payload) + { + if(tc>0) + { + triggerCellData_noBestChoice_->Fill(tc); + if (is_Simhit_comp_){ + if (TC_simhit_energies.at(icell_noBestChoice) >0){ + triggerCellSimHits_noBestChoice_->Fill(TC_simhit_energies.at(icell_noBestChoice)); + triggerCellData_noBestChoice_vsSimHits_->Fill(TC_simhit_energies.at(icell_noBestChoice),tc); + } + } + } + icell_noBestChoice++; + } + + // after sum, best choice, no encode/decode + int icell_BestChoice = 0; + for(const auto& tc : fe_payload_TCsums_BestChoice.payload) + { + if(tc>0) + { + triggerCellData_BestChoice_->Fill(tc); + if (is_Simhit_comp_){ + if (TC_simhit_energies.at(icell_BestChoice)>0) triggerCellData_BestChoice_vsSimHits_->Fill(TC_simhit_energies.at(icell_BestChoice),tc); + } + } + icell_BestChoice++; + } + // after sum, best choice, encode/decode + int icell = 0; size_t nFEDigi = 0; unsigned triggerCellModuleSum = 0; for(const auto& tc : fe_payload.payload) - { + { uint32_t tcShifted = (tc<<2); if(tcShifted>0) - { + { nFEDigi++; triggerCellModuleSum += tcShifted; triggerCellData_->Fill(tcShifted); - } - } + if (is_Simhit_comp_){ + if (TC_simhit_energies.at(icell)>0) triggerCellData_vsSimHits_->Fill(TC_simhit_energies.at(icell),tcShifted); + } + } + icell++; + } triggerCellsPerModule_->Fill(nFEDigi); triggerCellModuleSum_->Fill(triggerCellModuleSum); } diff --git a/L1Trigger/L1THGCal/test/testHGCalL1TBestChoice_cfg.py b/L1Trigger/L1THGCal/test/testHGCalL1TBestChoice_cfg.py index c4903b211c9a2..82dc62c2aac24 100644 --- a/L1Trigger/L1THGCal/test/testHGCalL1TBestChoice_cfg.py +++ b/L1Trigger/L1THGCal/test/testHGCalL1TBestChoice_cfg.py @@ -129,6 +129,10 @@ eeDigis = cms.InputTag('mix:HGCDigisEE'), fhDigis = cms.InputTag('mix:HGCDigisHEfront'), #bhDigis = cms.InputTag('mix:HGCDigisHEback'), + isSimhitComp = cms.bool(False), + eeSimHits = cms.InputTag('g4SimHits:HGCHitsEE'), + fhSimHits = cms.InputTag('g4SimHits:HGCHitsHEfront'), + #bhSimHits = cms.InputTag('g4SimHits:HGCHitsHEback'), beClustersAll = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:SingleCellClusterAlgo'), beClustersSelect = cms.InputTag('hgcalTriggerPrimitiveDigiFEReproducer:SingleCellClusterAlgo'), TriggerGeometry = cms.PSet(