diff --git a/DataFormats/L1THGCal/interface/HGCalClusterT.h b/DataFormats/L1THGCal/interface/HGCalClusterT.h index 21497483d7590..fffdf5d6bfa4e 100644 --- a/DataFormats/L1THGCal/interface/HGCalClusterT.h +++ b/DataFormats/L1THGCal/interface/HGCalClusterT.h @@ -37,17 +37,20 @@ namespace l1t centre_(0, 0, 0), centreProj_(0., 0., 0.), mipPt_(0), - seedMipPt_(0){} + seedMipPt_(0), + sumPt_(0){} - HGCalClusterT( const edm::Ptr& c ): + + HGCalClusterT( const edm::Ptr& c, float fraction=1. ): valid_(true), detId_( c->detId() ), centre_(0., 0., 0.), centreProj_(0., 0., 0.), mipPt_(0.), - seedMipPt_(0.) + seedMipPt_(0.), + sumPt_(0.) { - addConstituent(c); + addConstituent(c, true, fraction); } ~HGCalClusterT() override {}; @@ -105,11 +108,10 @@ namespace l1t double mipPt() const { return mipPt_; } double seedMipPt() const { return seedMipPt_; } uint32_t detId() const { return detId_.rawId(); } - void setPt(double pt) { setP4( math::PtEtaPhiMLorentzVector(pt, eta(), phi(), mass() ) ); } - + double sumPt() const { return sumPt_; } /* distance in 'cm' */ double distance( const l1t::HGCalTriggerCell &tc ) const { return ( tc.position() - centre_ ).mag(); } @@ -203,6 +205,7 @@ namespace l1t double mipPt_; double seedMipPt_; + double sumPt_; //shower shape @@ -226,6 +229,7 @@ namespace l1t void updateP4AndPosition(const edm::Ptr& c, bool updateCentre=true, float fraction=1.) { double cMipt = c->mipPt()*fraction; + double cPt = c->pt()*fraction; /* update cluster positions (IF requested) */ if( updateCentre ){ Basic3DVector constituentCentre( c->position() ); @@ -246,7 +250,7 @@ namespace l1t /* update cluster energies */ mipPt_ += cMipt; - + sumPt_ += cPt; int updatedPt = hwPt() + (int)(c->hwPt()*fraction); setHwPt( updatedPt ); diff --git a/DataFormats/L1THGCal/interface/HGCalMulticluster.h b/DataFormats/L1THGCal/interface/HGCalMulticluster.h index f8e4e4f062555..5f4c59bfa1ffa 100644 --- a/DataFormats/L1THGCal/interface/HGCalMulticluster.h +++ b/DataFormats/L1THGCal/interface/HGCalMulticluster.h @@ -19,7 +19,7 @@ namespace l1t { int phi=0 ); - HGCalMulticluster( const edm::Ptr &tc ); + HGCalMulticluster( const edm::Ptr &tc, float fraction=1); ~HGCalMulticluster() override; diff --git a/DataFormats/L1THGCal/src/HGCalMulticluster.cc b/DataFormats/L1THGCal/src/HGCalMulticluster.cc index d5b7bbc98f4bb..57e4c38a52a51 100644 --- a/DataFormats/L1THGCal/src/HGCalMulticluster.cc +++ b/DataFormats/L1THGCal/src/HGCalMulticluster.cc @@ -13,8 +13,8 @@ HGCalMulticluster::HGCalMulticluster( const LorentzVector p4, } -HGCalMulticluster::HGCalMulticluster( const edm::Ptr &clusterSeed ) - : HGCalClusterT(clusterSeed), +HGCalMulticluster::HGCalMulticluster( const edm::Ptr &clusterSeed, float fraction ) + : HGCalClusterT(clusterSeed, fraction), hOverE_(-99), hOverEValid_(false) { diff --git a/DataFormats/L1THGCal/src/HGCalTriggerCell.cc b/DataFormats/L1THGCal/src/HGCalTriggerCell.cc index 599fb76c07bf2..74a3d0036703d 100644 --- a/DataFormats/L1THGCal/src/HGCalTriggerCell.cc +++ b/DataFormats/L1THGCal/src/HGCalTriggerCell.cc @@ -10,7 +10,7 @@ HGCalTriggerCell( const LorentzVector& p4, int qual, uint32_t detid): L1Candidate(p4, pt, eta, phi, qual), - detid_(detid) + detid_(detid), position_(), mipPt_(0) { } diff --git a/DataFormats/L1THGCal/src/classes_def.xml b/DataFormats/L1THGCal/src/classes_def.xml index 5de052d7f893a..659c0682c681d 100644 --- a/DataFormats/L1THGCal/src/classes_def.xml +++ b/DataFormats/L1THGCal/src/classes_def.xml @@ -30,8 +30,9 @@ - - + + + @@ -43,9 +44,11 @@ - - - + + + + + diff --git a/L1Trigger/L1THGCal/interface/backend/HGCalHistoClusteringImpl.h b/L1Trigger/L1THGCal/interface/backend/HGCalHistoClusteringImpl.h new file mode 100644 index 0000000000000..4481d05661a61 --- /dev/null +++ b/L1Trigger/L1THGCal/interface/backend/HGCalHistoClusteringImpl.h @@ -0,0 +1,73 @@ +#ifndef __L1Trigger_L1THGCal_HGCalHistoClusteringImpl_h__ +#define __L1Trigger_L1THGCal_HGCalHistoClusteringImpl_h__ + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "DataFormats/L1THGCal/interface/HGCalCluster.h" +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" +#include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" +#include "L1Trigger/L1THGCal/interface/backend/HGCalTriggerClusterIdentificationBase.h" + + +class HGCalHistoClusteringImpl{ + +public: + + HGCalHistoClusteringImpl( const edm::ParameterSet &conf); + + void eventSetup(const edm::EventSetup& es) + { + triggerTools_.eventSetup(es); + shape_.eventSetup(es); + if ( (!dr_byLayer_coefficientA_.empty() && (dr_byLayer_coefficientA_.size()-1) < triggerTools_.lastLayerBH()) + || (!dr_byLayer_coefficientB_.empty() && (dr_byLayer_coefficientB_.size()-1) < triggerTools_.lastLayerBH()) + ) { + throw cms::Exception("Configuration") << + "The per-layer dR values go up to " << (dr_byLayer_coefficientA_.size()-1) << + "(A) and " << (dr_byLayer_coefficientB_.size()-1) << "(B), while layers go up to " << triggerTools_.lastLayerBH() << "\n"; + } + } + + float dR( const l1t::HGCalCluster & clu, + const GlobalPoint & seed ) const; + + void clusterizeHisto( const std::vector> & clustersPtr, + const std::vector > & seedPositionsEnergy, + const HGCalTriggerGeometryBase & triggerGeometry, + l1t::HGCalMulticlusterBxCollection & multiclusters + ) const; + + +private: + enum ClusterAssociationStrategy{ + NearestNeighbour, + EnergySplit + }; + + std::vector clusterSeedMulticluster(const std::vector> & clustersPtrs, + const std::vector > & seeds) const; + + void finalizeClusters(std::vector&, + l1t::HGCalMulticlusterBxCollection&, + const HGCalTriggerGeometryBase&) const; + + double dr_; + std::vector dr_byLayer_coefficientA_; + std::vector dr_byLayer_coefficientB_; + double ptC3dThreshold_; + + std::string cluster_association_input_; + ClusterAssociationStrategy cluster_association_strategy_; + + HGCalShowerShape shape_; + HGCalTriggerTools triggerTools_; + std::unique_ptr id_; + + static constexpr double kMidRadius_ = 2.3; + +}; + +#endif diff --git a/L1Trigger/L1THGCal/interface/backend/HGCalMulticlusteringHistoImpl.h b/L1Trigger/L1THGCal/interface/backend/HGCalHistoSeedingImpl.h similarity index 51% rename from L1Trigger/L1THGCal/interface/backend/HGCalMulticlusteringHistoImpl.h rename to L1Trigger/L1THGCal/interface/backend/HGCalHistoSeedingImpl.h index 6d9b8715c3a54..9fbc0583903d9 100644 --- a/L1Trigger/L1THGCal/interface/backend/HGCalMulticlusteringHistoImpl.h +++ b/L1Trigger/L1THGCal/interface/backend/HGCalHistoSeedingImpl.h @@ -1,5 +1,5 @@ -#ifndef __L1Trigger_L1THGCal_HGCalMulticlusteringHistoImpl_h__ -#define __L1Trigger_L1THGCal_HGCalMulticlusteringHistoImpl_h__ +#ifndef __L1Trigger_L1THGCal_HGCalHistoSeedingImpl_h__ +#define __L1Trigger_L1THGCal_HGCalHistoSeedingImpl_h__ #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "DataFormats/L1THGCal/interface/HGCalCluster.h" @@ -12,30 +12,28 @@ #include "L1Trigger/L1THGCal/interface/backend/HGCalTriggerClusterIdentificationBase.h" -class HGCalMulticlusteringHistoImpl{ +class HGCalHistoSeedingImpl{ public: - HGCalMulticlusteringHistoImpl( const edm::ParameterSet &conf); + HGCalHistoSeedingImpl( const edm::ParameterSet &conf); void eventSetup(const edm::EventSetup& es) { triggerTools_.eventSetup(es); - shape_.eventSetup(es); } float dR( const l1t::HGCalCluster & clu, - const GlobalPoint & seed ) const; + const GlobalPoint & seed ) const; - void clusterizeHisto( const std::vector> & clustersPtr, - l1t::HGCalMulticlusterBxCollection & multiclusters, - const HGCalTriggerGeometryBase & triggerGeometry - ); + void findHistoSeeds( const std::vector> & clustersPtr, + std::vector > & seedPositionsEnergy); private: - enum MulticlusterType{ + enum SeedingType{ HistoMaxC3d, + HistoSecondaryMaxC3d, HistoThresholdC3d, HistoInterpolatedMaxC3d }; @@ -45,36 +43,29 @@ class HGCalMulticlusteringHistoImpl{ Histogram fillHistoClusters( const std::vector> & clustersPtrs ); Histogram fillSmoothPhiHistoClusters( const Histogram & histoClusters, - const vector & binSums ); + const vector & binSums ); Histogram fillSmoothRPhiHistoClusters( const Histogram & histoClusters ); - std::vector computeMaxSeeds( const Histogram & histoClusters ); + std::vector > computeMaxSeeds( const Histogram & histoClusters ); - std::vector computeInterpolatedMaxSeeds( const Histogram & histoClusters ); + std::vector > computeSecondaryMaxSeeds( const Histogram & histoClusters ); + + std::vector > computeInterpolatedMaxSeeds( const Histogram & histoClusters ); + + std::vector > computeThresholdSeeds( const Histogram & histoClusters ); - std::vector computeThresholdSeeds( const Histogram & histoClusters ); - std::vector clusterSeedMulticluster(const std::vector> & clustersPtrs, - const std::vector & seeds); + std::string seedingAlgoType_; + SeedingType seedingType_; - void finalizeClusters(std::vector&, - l1t::HGCalMulticlusterBxCollection&, - const HGCalTriggerGeometryBase&); - - double dr_; - double ptC3dThreshold_; - MulticlusterType multiclusteringAlgoType_; - std::string multiclusterAlgoType_; unsigned nBinsRHisto_ = 36; unsigned nBinsPhiHisto_ = 216; std::vector binsSumsHisto_; double histoThreshold_ = 20.; std::vector neighbour_weights_; - HGCalShowerShape shape_; HGCalTriggerTools triggerTools_; - std::unique_ptr id_; static constexpr unsigned neighbour_weights_size_ = 9; static constexpr double kROverZMin_ = 0.09; diff --git a/L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h b/L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h index 61b0f7ef6db4a..ba234be8a11c2 100644 --- a/L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h +++ b/L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h @@ -44,6 +44,8 @@ class HGCalShowerShape{ float sigmaRRMax(const l1t::HGCalMulticluster& c3d) const; float sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius=5.) const; + void fillShapes(l1t::HGCalMulticluster&, const HGCalTriggerGeometryBase&) const; + private: float meanX(const std::vector >& energy_X_tc) const; diff --git a/L1Trigger/L1THGCal/interface/concentrator/HGCalConcentratorSuperTriggerCellImpl.h b/L1Trigger/L1THGCal/interface/concentrator/HGCalConcentratorSuperTriggerCellImpl.h index efb586e7a5224..f12504dc447de 100644 --- a/L1Trigger/L1THGCal/interface/concentrator/HGCalConcentratorSuperTriggerCellImpl.h +++ b/L1Trigger/L1THGCal/interface/concentrator/HGCalConcentratorSuperTriggerCellImpl.h @@ -6,9 +6,12 @@ #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h" #include "DataFormats/L1THGCal/interface/HGCalTriggerSums.h" +#include "DataFormats/ForwardDetId/interface/HGCSiliconDetIdToROC.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" + + #include #include @@ -18,12 +21,28 @@ class HGCalConcentratorSuperTriggerCellImpl HGCalConcentratorSuperTriggerCellImpl(const edm::ParameterSet& conf); void superTriggerCellSelectImpl(const std::vector& trigCellVecInput, std::vector& trigCellVecOutput); + void eventSetup(const edm::EventSetup& es) {triggerTools_.eventSetup(es);} private: int getSuperTriggerCellId(int detid) const ; - static const int kSplit_ = 0x3a; + static std::map kSplit_; static const int kWafer_offset_ = 6; + static const int kSTCsizeCoarse_ = 16; + static const int kSTCsizeFine_ = 4; + static const int kSplit_v8_Coarse_ = 0x30; + static const int kSplit_v8_Fine_ = 0x3a; + static const int kNLayers_ = 3; + static const int kSplit_v9_ = 0x36; + + static const int kRocShift_ = 6; + static const int kRotate4_ = 4; + static const int kUShift_ = 3; + + + HGCalTriggerTools triggerTools_; + HGCSiliconDetIdToROC detIdToROC_; + std::vector stcSize_; class SuperTriggerCell { @@ -31,7 +50,7 @@ class HGCalConcentratorSuperTriggerCellImpl float sumPt_, sumMipPt_; int sumHwPt_, maxHwPt_; unsigned maxId_; - + public: SuperTriggerCell(){ sumPt_=0, sumMipPt_=0, sumHwPt_=0, maxHwPt_=0, maxId_=0 ;} void add(const l1t::HGCalTriggerCell &c) { diff --git a/L1Trigger/L1THGCal/plugins/BuildFile.xml b/L1Trigger/L1THGCal/plugins/BuildFile.xml index 9961a5d0a3bbd..028150e8acb3c 100644 --- a/L1Trigger/L1THGCal/plugins/BuildFile.xml +++ b/L1Trigger/L1THGCal/plugins/BuildFile.xml @@ -29,18 +29,3 @@ - - - - - - - - - - - - - - - diff --git a/L1Trigger/L1THGCal/plugins/backend/HGCalBackendLayer2Processor3DClustering.cc b/L1Trigger/L1THGCal/plugins/backend/HGCalBackendLayer2Processor3DClustering.cc index d1aecb1b41ad8..c0676d3995f69 100644 --- a/L1Trigger/L1THGCal/plugins/backend/HGCalBackendLayer2Processor3DClustering.cc +++ b/L1Trigger/L1THGCal/plugins/backend/HGCalBackendLayer2Processor3DClustering.cc @@ -6,7 +6,8 @@ #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" #include "L1Trigger/L1THGCal/interface/backend/HGCalMulticlusteringImpl.h" -#include "L1Trigger/L1THGCal/interface/backend/HGCalMulticlusteringHistoImpl.h" +#include "L1Trigger/L1THGCal/interface/backend/HGCalHistoSeedingImpl.h" +#include "L1Trigger/L1THGCal/interface/backend/HGCalHistoClusteringImpl.h" class HGCalBackendLayer2Processor3DClustering : public HGCalBackendLayer2ProcessorBase { @@ -23,7 +24,8 @@ class HGCalBackendLayer2Processor3DClustering : public HGCalBackendLayer2Process multiclustering_ = std::make_unique( conf.getParameterSet("C3d_parameters") ); }else if(typeMulticluster.find("Histo")!=std::string::npos){ multiclusteringAlgoType_ = HistoC3d; - multiclusteringHisto_ = std::make_unique( conf.getParameterSet("C3d_parameters") ); + multiclusteringHistoSeeding_ = std::make_unique( conf.getParameterSet("C3d_parameters") ); + multiclusteringHistoClustering_ = std::make_unique( conf.getParameterSet("C3d_parameters") ); }else { throw cms::Exception("HGCTriggerParameterError") << "Unknown Multiclustering type '" << typeMulticluster << "'"; @@ -36,7 +38,8 @@ class HGCalBackendLayer2Processor3DClustering : public HGCalBackendLayer2Process { es.get().get("", triggerGeometry_); if(multiclustering_) multiclustering_->eventSetup(es); - if(multiclusteringHisto_) multiclusteringHisto_->eventSetup(es); + if(multiclusteringHistoSeeding_) multiclusteringHistoSeeding_->eventSetup(es); + if(multiclusteringHistoClustering_) multiclusteringHistoClustering_->eventSetup(es); /* create a persistent vector of pointers to the trigger-cells */ std::vector> clustersPtrs; @@ -45,6 +48,9 @@ class HGCalBackendLayer2Processor3DClustering : public HGCalBackendLayer2Process clustersPtrs.push_back(ptr); } + /* create a vector of seed positions and their energy*/ + std::vector > seedPositionsEnergy; + /* call to multiclustering and compute shower shape*/ switch(multiclusteringAlgoType_){ case dRC3d : @@ -54,7 +60,8 @@ class HGCalBackendLayer2Processor3DClustering : public HGCalBackendLayer2Process multiclustering_->clusterizeDBSCAN( clustersPtrs, collCluster3D, *triggerGeometry_); break; case HistoC3d : - multiclusteringHisto_->clusterizeHisto( clustersPtrs, collCluster3D, *triggerGeometry_); + multiclusteringHistoSeeding_->findHistoSeeds( clustersPtrs, seedPositionsEnergy); + multiclusteringHistoClustering_->clusterizeHisto( clustersPtrs, seedPositionsEnergy, *triggerGeometry_, collCluster3D ); break; default: // Should not happen, clustering type checked in constructor @@ -73,7 +80,8 @@ class HGCalBackendLayer2Processor3DClustering : public HGCalBackendLayer2Process /* algorithms instances */ std::unique_ptr multiclustering_; - std::unique_ptr multiclusteringHisto_; + std::unique_ptr multiclusteringHistoSeeding_; + std::unique_ptr multiclusteringHistoClustering_; /* algorithm type */ MulticlusterType multiclusteringAlgoType_; diff --git a/L1Trigger/L1THGCal/plugins/concentrator/HGCalConcentratorProcessorSelection.cc b/L1Trigger/L1THGCal/plugins/concentrator/HGCalConcentratorProcessorSelection.cc index 5fa5569d91b94..32ea6af2667f7 100644 --- a/L1Trigger/L1THGCal/plugins/concentrator/HGCalConcentratorProcessorSelection.cc +++ b/L1Trigger/L1THGCal/plugins/concentrator/HGCalConcentratorProcessorSelection.cc @@ -41,6 +41,7 @@ void HGCalConcentratorProcessorSelection::run(const edm::HandlegetModuleFromTriggerCell(trigCell.detId()); tc_modules[module].push_back(trigCell); } + if ( concentratorSTCImpl_) concentratorSTCImpl_->eventSetup(es); for( const auto& module_trigcell : tc_modules ) { std::vector trigCellVecOutput; diff --git a/L1Trigger/L1THGCal/python/customClustering.py b/L1Trigger/L1THGCal/python/customClustering.py index 75b711b10589d..c7696c0d528b9 100644 --- a/L1Trigger/L1THGCal/python/customClustering.py +++ b/L1Trigger/L1THGCal/python/customClustering.py @@ -8,6 +8,22 @@ 3, 3, 3, 3, 3, 3, 3, 3 # 28 - 35 ) +EE_DR_GROUP = 7 +FH_DR_GROUP = 6 +BH_DR_GROUP = 12 +MAX_LAYERS = 52 + +dr_layerbylayer = ([0] + # no layer 0 + [0.015]*EE_DR_GROUP + [0.020]*EE_DR_GROUP + [0.030]*EE_DR_GROUP + [0.040]*EE_DR_GROUP + # EM + [0.040]*FH_DR_GROUP + [0.050]*FH_DR_GROUP + # FH + [0.050]*BH_DR_GROUP) # BH + + +dr_layerbylayer_Bcoefficient = ([0] + # no layer 0 + [0.020]*EE_DR_GROUP + [0.020]*EE_DR_GROUP + [0.02]*EE_DR_GROUP + [0.020]*EE_DR_GROUP + # EM + [0.020]*FH_DR_GROUP + [0.020]*FH_DR_GROUP + # FH + [0.020]*BH_DR_GROUP) # BH + def custom_2dclustering_distance(process, distance=6.,# cm @@ -76,58 +92,91 @@ def custom_3dclustering_dbscan(process, def custom_3dclustering_histoMax(process, - distance = 0.01, + distance = 0.03, nBins_R = 36, nBins_Phi = 216, binSumsHisto = binSums, + seed_threshold = 10, ): parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters parameters_c3d.dR_multicluster = cms.double(distance) parameters_c3d.nBins_R_histo_multicluster = cms.uint32(nBins_R) parameters_c3d.nBins_Phi_histo_multicluster = cms.uint32(nBins_Phi) parameters_c3d.binSumsHisto = binSumsHisto + parameters_c3d.threshold_histo_multicluster = seed_threshold parameters_c3d.type_multicluster = cms.string('HistoMaxC3d') return process +def custom_3dclustering_histoSecondaryMax(process, + distance = 0.03, + threshold = 5., + nBins_R = 36, + nBins_Phi = 216, + binSumsHisto = binSums, + ): + parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters + parameters_c3d.dR_multicluster = cms.double(distance) + parameters_c3d.nBins_R_histo_multicluster = cms.uint32(nBins_R) + parameters_c3d.nBins_Phi_histo_multicluster = cms.uint32(nBins_Phi) + parameters_c3d.binSumsHisto = binSumsHisto + parameters_c3d.threshold_histo_multicluster = cms.double(threshold) + parameters_c3d.type_multicluster = cms.string('HistoSecondaryMaxC3d') + return process + + +def custom_3dclustering_histoMax_variableDr(process, + distances = dr_layerbylayer, + nBins_R = 36, + nBins_Phi = 216, + binSumsHisto = binSums, + seed_threshold = 10, + ): + process = custom_3dclustering_histoMax(process, 0, nBins_R, nBins_Phi, binSumsHisto, seed_threshold) + parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters + parameters_c3d.dR_multicluster_byLayer_coefficientA = cms.vdouble(distances) + return process + + def custom_3dclustering_histoInterpolatedMax(process, - distance = 0.01, + distance = 0.03, + seed_threshold = 5., nBins_R = 36, nBins_Phi = 216, binSumsHisto = binSums, ): - process = custom_3dclustering_histoMax( process, distance, nBins_R, nBins_Phi, binSumsHisto ) + process = custom_3dclustering_histoMax( process, distance, nBins_R, nBins_Phi, binSumsHisto, seed_threshold ) parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters parameters_c3d.type_multicluster = cms.string('HistoInterpolatedMaxC3d') return process -def custom_3dclustering_histoInterpolatedMax1stOrder(process): +def custom_3dclustering_histoInterpolatedMax1stOrder(process, distance = 0.03, seed_threshold = 5. ): parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters parameters_c3d.neighbour_weights=cms.vdouble( 0 , 0.25, 0 , 0.25 , 0 , 0.25, 0 , 0.25, 0 ) - process = custom_3dclustering_histoInterpolatedMax( process ) + process = custom_3dclustering_histoInterpolatedMax( process, distance, seed_threshold ) return process -def custom_3dclustering_histoInterpolatedMax2ndOrder(process): +def custom_3dclustering_histoInterpolatedMax2ndOrder(process, distance = 0.03, seed_threshold = 5.): parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters parameters_c3d.neighbour_weights=cms.vdouble( -0.25, 0.5, -0.25, 0.5 , 0 , 0.5 , -0.25, 0.5, -0.25 ) - process = custom_3dclustering_histoInterpolatedMax( process ) + process = custom_3dclustering_histoInterpolatedMax( process, distance, seed_threshold ) return process def custom_3dclustering_histoThreshold(process, threshold = 20., - distance = 0.01, + distance = 0.03, nBins_R = 36, nBins_Phi = 216, binSumsHisto = binSums, @@ -140,3 +189,50 @@ def custom_3dclustering_histoThreshold(process, parameters_c3d.binSumsHisto = binSumsHisto parameters_c3d.type_multicluster = cms.string('HistoThresholdC3d') return process + +def custom_3dclustering_clusteringRadiusLayerbyLayerVariableEta(process, distance_coefficientA = dr_layerbylayer, distance_coefficientB = dr_layerbylayer_Bcoefficient): + + parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters + parameters_c3d.dR_multicluster_byLayer_coefficientA = distance_coefficientA + parameters_c3d.dR_multicluster_byLayer_coefficientB = distance_coefficientB + + return process + + +def custom_3dclustering_clusteringRadiusLayerbyLayerFixedEta(process, distance_coefficientA = dr_layerbylayer): + + parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters + parameters_c3d.dR_multicluster_byLayer_coefficientA = distance_coefficientA + parameters_c3d.dR_multicluster_byLayer_coefficientB = cms.vdouble( [0]*(MAX_LAYERS+1) ) + + return process + +def custom_3dclustering_clusteringRadiusNoLayerDependenceFixedEta(process, distance_coefficientA = 0.03): + + parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters + parameters_c3d.dR_multicluster_byLayer_coefficientA = cms.vdouble( [distance_coefficientA]*(MAX_LAYERS+1) ) + parameters_c3d.dR_multicluster_byLayer_coefficientB = cms.vdouble( [0]*(MAX_LAYERS+1) ) + + return process + +def custom_3dclustering_clusteringRadiusNoLayerDependenceVariableEta(process, distance_coefficientA = 0.03, distance_coefficientB = 0.02): + + parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters + parameters_c3d.dR_multicluster_byLayer_coefficientA = cms.vdouble( [distance_coefficientA]*(MAX_LAYERS+1) ) + parameters_c3d.dR_multicluster_byLayer_coefficientB = cms.vdouble( [distance_coefficientB]*(MAX_LAYERS+1) ) + + return process + + +def custom_3dclustering_nearestNeighbourAssociation(process): + + parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters + parameters_c3d.cluster_association = cms.string('NearestNeighbour') + + return process + +def custom_3dclustering_EnergySplitAssociation(process): + + parameters_c3d = process.hgcalBackEndLayer2Producer.ProcessorParameters.C3d_parameters + parameters_c3d.cluster_association = cms.string('EnergySplit') + return process diff --git a/L1Trigger/L1THGCal/python/customNtuples.py b/L1Trigger/L1THGCal/python/customNtuples.py deleted file mode 100644 index 8ada2eb9e7c5e..0000000000000 --- a/L1Trigger/L1THGCal/python/customNtuples.py +++ /dev/null @@ -1,9 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -def custom_ntuples_V9(process): - ntuples = process.hgcalTriggerNtuplizer.Ntuples - for ntuple in ntuples: - if ntuple.NtupleName=='HGCalTriggerNtupleHGCDigis' or \ - ntuple.NtupleName=='HGCalTriggerNtupleHGCTriggerCells': - ntuple.bhSimHits = cms.InputTag('g4SimHits:HGCHitsHEback') - return process diff --git a/L1Trigger/L1THGCal/python/customTriggerCellSelect.py b/L1Trigger/L1THGCal/python/customTriggerCellSelect.py index d6a66d04663ea..7473ef1b7d4ce 100644 --- a/L1Trigger/L1THGCal/python/customTriggerCellSelect.py +++ b/L1Trigger/L1THGCal/python/customTriggerCellSelect.py @@ -1,8 +1,14 @@ import FWCore.ParameterSet.Config as cms import SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi as digiparam -def custom_triggercellselect_supertriggercell(process): - process.hgcalConcentratorProducer.ProcessorParameters.Method = cms.string('superTriggerCellSelect') +def custom_triggercellselect_supertriggercell(process, + stcSize = cms.vuint32(4,4,4) + ): + + parameters = process.hgcalConcentratorProducer.ProcessorParameters + parameters.Method = cms.string('superTriggerCellSelect') + parameters.stcSize = stcSize + return process diff --git a/L1Trigger/L1THGCal/python/hgcalBackEndLayer1Producer_cfi.py b/L1Trigger/L1THGCal/python/hgcalBackEndLayer1Producer_cfi.py index d418790c84368..38dce63e5f3ce 100644 --- a/L1Trigger/L1THGCal/python/hgcalBackEndLayer1Producer_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalBackEndLayer1Producer_cfi.py @@ -6,13 +6,8 @@ import RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi as recocalibparam from . import hgcalLayersCalibrationCoefficients_cfi as layercalibparam -C2d_parValues = cms.PSet( clusterType = cms.string('dRNNC2d'), # clustering type: dRC2d--> Geometric-dR clustering; NNC2d-->Nearest Neighbors clustering - seeding_threshold_silicon = cms.double(5), # MipT - seeding_threshold_scintillator = cms.double(5), # MipT - clustering_threshold_silicon = cms.double(2), # MipT - clustering_threshold_scintillator = cms.double(2), # MipT - dR_cluster = cms.double(6.), # in cm - calibSF_cluster=cms.double(0.), +C2d_parValues = cms.PSet( clusterType = cms.string('dummyC2d'), + calibSF_cluster=cms.double(1.), layerWeights = layercalibparam.TrgLayer_weights, applyLayerCalibration = cms.bool(True) ) diff --git a/L1Trigger/L1THGCal/python/hgcalBackEndLayer2Producer_cfi.py b/L1Trigger/L1THGCal/python/hgcalBackEndLayer2Producer_cfi.py index 2b710c15217e5..1caca9a077339 100644 --- a/L1Trigger/L1THGCal/python/hgcalBackEndLayer2Producer_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalBackEndLayer2Producer_cfi.py @@ -5,17 +5,18 @@ import RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi as recocalibparam from L1Trigger.L1THGCal.egammaIdentification import egamma_identification_drnn_cone -from L1Trigger.L1THGCal.customClustering import binSums +from L1Trigger.L1THGCal.customClustering import binSums, dr_layerbylayer -C3d_parValues = cms.PSet( type_multicluster = cms.string('dRC3d'), - dR_multicluster = cms.double(0.01), +C3d_parValues = cms.PSet( type_multicluster = cms.string('HistoMaxC3d'), + dR_multicluster = cms.double(0.), + dR_multicluster_byLayer_coefficientA = cms.vdouble(dr_layerbylayer), + dR_multicluster_byLayer_coefficientB = cms.vdouble( [0]*53 ), minPt_multicluster = cms.double(0.5), # minimum pt of the multicluster (GeV) - dist_dbscan_multicluster=cms.double(0.), - minN_dbscan_multicluster=cms.uint32(0), nBins_R_histo_multicluster = cms.uint32(36), nBins_Phi_histo_multicluster = cms.uint32(216), binSumsHisto = binSums, - threshold_histo_multicluster = cms.double(20.), + threshold_histo_multicluster = cms.double(10.), + cluster_association = cms.string("NearestNeighbour"), EGIdentification=egamma_identification_drnn_cone.clone(), neighbour_weights=cms.vdouble( 0 , 0.25, 0 , 0.25 , 0 , 0.25, diff --git a/L1Trigger/L1THGCal/python/hgcalConcentratorProducer_cfi.py b/L1Trigger/L1THGCal/python/hgcalConcentratorProducer_cfi.py index 7c52c58537730..d6e19f843078f 100644 --- a/L1Trigger/L1THGCal/python/hgcalConcentratorProducer_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalConcentratorProducer_cfi.py @@ -15,7 +15,7 @@ TCThreshold_fC = cms.double(0.), TCThresholdBH_MIP = cms.double(0.), triggercell_threshold_silicon = cms.double(2.), # MipT - triggercell_threshold_scintillator = cms.double(2.) # MipT + triggercell_threshold_scintillator = cms.double(2.), # MipT ) hgcalConcentratorProducer = cms.EDProducer( diff --git a/L1Trigger/L1THGCal/python/hgcalLayersCalibrationCoefficients_cfi.py b/L1Trigger/L1THGCal/python/hgcalLayersCalibrationCoefficients_cfi.py index ebd349289ec37..f1f31913e26e5 100644 --- a/L1Trigger/L1THGCal/python/hgcalLayersCalibrationCoefficients_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalLayersCalibrationCoefficients_cfi.py @@ -55,7 +55,7 @@ 0.326339 ) -TrgLayer_weights = cms.vdouble(0.0, +TrgLayer_weights = cms.vdouble(0.0, 0.0183664, 0., 0.0305622, @@ -110,56 +110,16 @@ 0.328053 ) -TrgLayer_dEdX_weights = cms.vdouble(0.0, # there is no layer zero - 8.603+8.0675, # Mev - 0., - 8.0675*2, - 0., - 8.0675*2, - 0., - 8.0675*2., - 0., - 8.0675+8.9515, - 0., - 10.135*2, - 0., - 10.135*2., - 0., - 10.135*2., - 0., - 10.135*2., - 0., - 10.135+11.682, - 0., - 13.654*2., - 0., - 13.654*2., - 0., - 13.654*2., - 0., - 13.654+38.2005, - 0., - 55.0265, - 49.871, - 49.871, - 49.871, - 49.871, - 49.871, - 49.871, - 49.871, - 49.871, - 49.871, - 49.871, - 62.005, - 83.1675, - 92.196, - 92.196, - 92.196, - 92.196, - 92.196, - 92.196, - 92.196, - 92.196, - 92.196, - 92.196, - 46.098) +from RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi import dEdX + +EE_LAYERS = 28 +TrgLayer_dEdX_weights = cms.vdouble() +for lid, lw in enumerate(dEdX.weights): + if lid > (EE_LAYERS+1): + TrgLayer_dEdX_weights.append(lw) + else: + # Only half the layers are read in the EE at L1T + if (lid % 2) == 1: + TrgLayer_dEdX_weights.append(lw+dEdX.weights[lid-1]) + else: + TrgLayer_dEdX_weights.append(0) diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerPrimitives_cff.py b/L1Trigger/L1THGCal/python/hgcalTriggerPrimitives_cff.py index cad730b75b5d4..bf6347c9c4f2f 100644 --- a/L1Trigger/L1THGCal/python/hgcalTriggerPrimitives_cff.py +++ b/L1Trigger/L1THGCal/python/hgcalTriggerPrimitives_cff.py @@ -1,19 +1,21 @@ import FWCore.ParameterSet.Config as cms from L1Trigger.L1THGCal.hgcalTriggerGeometryESProducer_cfi import * -from L1Trigger.L1THGCal.hgcalVFEProducer_cfi import * -from L1Trigger.L1THGCal.hgcalConcentratorProducer_cfi import * -from L1Trigger.L1THGCal.hgcalBackEndLayer1Producer_cfi import * -from L1Trigger.L1THGCal.hgcalBackEndLayer2Producer_cfi import * -from L1Trigger.L1THGCal.hgcalTowerMapProducer_cfi import * -from L1Trigger.L1THGCal.hgcalTowerProducer_cfi import * +from L1Trigger.L1THGCal.hgcalVFE_cff import * +from L1Trigger.L1THGCal.hgcalConcentrator_cff import * +from L1Trigger.L1THGCal.hgcalBackEndLayer1_cff import * +from L1Trigger.L1THGCal.hgcalBackEndLayer2_cff import * +from L1Trigger.L1THGCal.hgcalTowerMap_cff import * +from L1Trigger.L1THGCal.hgcalTower_cff import * -hgcalTriggerPrimitives = cms.Sequence(hgcalVFEProducer*hgcalConcentratorProducer*hgcalBackEndLayer1Producer*hgcalBackEndLayer2Producer*hgcalTowerMapProducer*hgcalTowerProducer) +hgcalTriggerPrimitives = cms.Sequence(hgcalVFE*hgcalConcentrator*hgcalBackEndLayer1*hgcalBackEndLayer2*hgcalTowerMap*hgcalTower) from Configuration.Eras.Modifier_phase2_hgcalV9_cff import phase2_hgcalV9 from L1Trigger.L1THGCal.customTriggerGeometry import custom_geometry_V9 +from L1Trigger.L1THGCal.customCalibration import custom_cluster_calibration_global modifyHgcalTriggerPrimitivesWithV9Geometry_ = phase2_hgcalV9.makeProcessModifier(custom_geometry_V9) +modifyHgcalTriggerPrimitivesCalibWithV9Geometry_ = phase2_hgcalV9.makeProcessModifier(lambda process : custom_cluster_calibration_global(process, factor=1)) from Configuration.ProcessModifiers.convertHGCalDigisSim_cff import convertHGCalDigisSim # can't declare a producer version of simHGCalUnsuppressedDigis in the normal flow of things, diff --git a/L1Trigger/L1THGCal/src/backend/HGCalHistoClusteringImpl.cc b/L1Trigger/L1THGCal/src/backend/HGCalHistoClusteringImpl.cc new file mode 100644 index 0000000000000..753b73878f537 --- /dev/null +++ b/L1Trigger/L1THGCal/src/backend/HGCalHistoClusteringImpl.cc @@ -0,0 +1,171 @@ +#include "L1Trigger/L1THGCal/interface/backend/HGCalHistoClusteringImpl.h" +#include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h" +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/Math/interface/deltaPhi.h" + +HGCalHistoClusteringImpl::HGCalHistoClusteringImpl( const edm::ParameterSet& conf ) : + dr_(conf.getParameter("dR_multicluster")), + dr_byLayer_coefficientA_(conf.existsAs>("dR_multicluster_byLayer_coefficientA") ? conf.getParameter>("dR_multicluster_byLayer_coefficientA") : std::vector()), + dr_byLayer_coefficientB_(conf.existsAs>("dR_multicluster_byLayer_coefficientB") ? conf.getParameter>("dR_multicluster_byLayer_coefficientB") : std::vector()), + ptC3dThreshold_(conf.getParameter("minPt_multicluster")), + cluster_association_input_(conf.getParameter("cluster_association")) +{ + + if(cluster_association_input_=="NearestNeighbour"){ + cluster_association_strategy_ = NearestNeighbour; + }else if(cluster_association_input_=="EnergySplit"){ + cluster_association_strategy_ = EnergySplit; + }else { + throw cms::Exception("HGCTriggerParameterError") + << "Unknown cluster association strategy'" << cluster_association_strategy_; + } + + edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster dR: " << dr_ + <<"\nMulticluster minimum transverse-momentum: " << ptC3dThreshold_ ; + + id_ = std::unique_ptr{ HGCalTriggerClusterIdentificationFactory::get()->create("HGCalTriggerClusterIdentificationBDT") }; + id_->initialize(conf.getParameter("EGIdentification")); + +} + + + +float HGCalHistoClusteringImpl::dR( const l1t::HGCalCluster & clu, + const GlobalPoint & seed) const +{ + + Basic3DVector seed_3dv( seed ); + GlobalPoint seed_proj( seed_3dv / seed.z() ); + return (seed_proj - clu.centreProj() ).mag(); + +} + + + +std::vector HGCalHistoClusteringImpl::clusterSeedMulticluster(const std::vector> & clustersPtrs, + const std::vector > & seeds) const +{ + + std::map mapSeedMulticluster; + std::vector multiclustersTmp; + + for(const auto& clu : clustersPtrs){ + + int z_side = triggerTools_.zside(clu->detId()); + + + double radiusCoefficientA = dr_byLayer_coefficientA_.empty() ? dr_ : dr_byLayer_coefficientA_[triggerTools_.layerWithOffset(clu->detId())]; + double radiusCoefficientB = dr_byLayer_coefficientB_.empty() ? 0 : dr_byLayer_coefficientB_[triggerTools_.layerWithOffset(clu->detId())]; + + double minDist = radiusCoefficientA + radiusCoefficientB*(kMidRadius_ - std::abs(clu->eta()) ) ; + + std::vector > targetSeedsEnergy; + + for( unsigned int iseed=0; iseeddR(*clu, seeds[iseed].first); + + if ( d < minDist ){ + if ( cluster_association_strategy_ == EnergySplit ){ + targetSeedsEnergy.emplace_back( iseed, seedEnergy ); + } + else if ( cluster_association_strategy_ == NearestNeighbour ){ + + minDist = d; + + if ( targetSeedsEnergy.empty() ) { + targetSeedsEnergy.emplace_back( iseed, seedEnergy ); + } + else { + targetSeedsEnergy.at(0).first = iseed ; + targetSeedsEnergy.at(0).second = seedEnergy; + } + } + + } + + } + + if(targetSeedsEnergy.empty()) continue; + //Loop over target seeds and divide up the clusters energy + double totalTargetSeedEnergy = 0; + for (const auto& energy: targetSeedsEnergy){ + totalTargetSeedEnergy+=energy.second; + } + + for (const auto& energy: targetSeedsEnergy){ + + double seedWeight = 1; + if ( cluster_association_strategy_ == EnergySplit) seedWeight = energy.second/totalTargetSeedEnergy; + if( mapSeedMulticluster[energy.first].size()==0) { + mapSeedMulticluster[energy.first] = l1t::HGCalMulticluster(clu, seedWeight) ; + } + else{ + mapSeedMulticluster[energy.first].addConstituent(clu, true, seedWeight); + } + + } + + } + + + for(const auto& mclu : mapSeedMulticluster) multiclustersTmp.emplace_back(mclu.second); + + return multiclustersTmp; + +} + + + + +void HGCalHistoClusteringImpl::clusterizeHisto( const std::vector> & clustersPtrs, + const std::vector > & seedPositionsEnergy, + const HGCalTriggerGeometryBase & triggerGeometry, + l1t::HGCalMulticlusterBxCollection & multiclusters) const +{ + + + /* clusterize clusters around seeds */ + std::vector multiclustersTmp = clusterSeedMulticluster(clustersPtrs,seedPositionsEnergy); + /* making the collection of multiclusters */ + finalizeClusters(multiclustersTmp, multiclusters, triggerGeometry); + +} + + + + + +void +HGCalHistoClusteringImpl:: +finalizeClusters(std::vector& multiclusters_in, + l1t::HGCalMulticlusterBxCollection& multiclusters_out, + const HGCalTriggerGeometryBase& triggerGeometry) const +{ + for(auto& multicluster : multiclusters_in) { + // compute the eta, phi from its barycenter + // + pT as scalar sum of pT of constituents + double sumPt=multicluster.sumPt(); + + math::PtEtaPhiMLorentzVector multiclusterP4( sumPt, + multicluster.centre().eta(), + multicluster.centre().phi(), + 0. ); + multicluster.setP4( multiclusterP4 ); + + if( multicluster.pt() > ptC3dThreshold_ ){ + //compute shower shapes + shape_.fillShapes(multicluster, triggerGeometry); + // fill quality flag + multicluster.setHwQual(id_->decision(multicluster)); + // fill H/E + multicluster.saveHOverE(); + + multiclusters_out.push_back( 0, multicluster); + } + } +} diff --git a/L1Trigger/L1THGCal/src/backend/HGCalHistoSeedingImpl.cc b/L1Trigger/L1THGCal/src/backend/HGCalHistoSeedingImpl.cc new file mode 100644 index 0000000000000..ef36f350fb132 --- /dev/null +++ b/L1Trigger/L1THGCal/src/backend/HGCalHistoSeedingImpl.cc @@ -0,0 +1,462 @@ +#include "L1Trigger/L1THGCal/interface/backend/HGCalHistoSeedingImpl.h" +#include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h" +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/Math/interface/deltaPhi.h" + +HGCalHistoSeedingImpl::HGCalHistoSeedingImpl( const edm::ParameterSet& conf ) : + seedingAlgoType_(conf.getParameter("type_multicluster")), + nBinsRHisto_(conf.getParameter("nBins_R_histo_multicluster")), + nBinsPhiHisto_(conf.getParameter("nBins_Phi_histo_multicluster")), + binsSumsHisto_(conf.getParameter< std::vector >("binSumsHisto")), + histoThreshold_(conf.getParameter("threshold_histo_multicluster")), + neighbour_weights_(conf.getParameter< std::vector >("neighbour_weights")) +{ + + if(seedingAlgoType_=="HistoMaxC3d"){ + seedingType_ = HistoMaxC3d; + }else if(seedingAlgoType_=="HistoSecondaryMaxC3d"){ + seedingType_ = HistoSecondaryMaxC3d; + }else if(seedingAlgoType_=="HistoThresholdC3d"){ + seedingType_ = HistoThresholdC3d; + }else if(seedingAlgoType_=="HistoInterpolatedMaxC3d"){ + seedingType_ = HistoInterpolatedMaxC3d; + }else { + throw cms::Exception("HGCTriggerParameterError") + << "Unknown Multiclustering type '" << seedingAlgoType_; + } + + edm::LogInfo("HGCalMulticlusterParameters") <<"\nMulticluster number of R-bins for the histo algorithm: " << nBinsRHisto_ + <<"\nMulticluster number of Phi-bins for the histo algorithm: " << nBinsPhiHisto_ + <<"\nMulticluster MIPT threshold for histo threshold algorithm: " << histoThreshold_ + <<"\nMulticluster type of multiclustering algortihm: " << seedingAlgoType_; + + if(seedingAlgoType_.find("Histo")!=std::string::npos && nBinsRHisto_!=binsSumsHisto_.size()){ + throw cms::Exception("Inconsistent bin size") << "Inconsistent nBins_R_histo_multicluster ( " << nBinsRHisto_ << " ) and binSumsHisto ( " << binsSumsHisto_.size() << " ) size in HGCalMulticlustering\n"; + } + + if(neighbour_weights_.size()!=neighbour_weights_size_) { + throw cms::Exception("Inconsistent vector size" ) << "Inconsistent size of neighbour weights vector in HGCalMulticlustering ( " << neighbour_weights_.size() << " ). Should be " << neighbour_weights_size_ << "\n" ; + } + +} + + +HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillHistoClusters( const std::vector> & clustersPtrs ){ + + + Histogram histoClusters; //key[0] = z.side(), key[1] = bin_R, key[2] = bin_phi + + for(int z_side : {-1,1}){ + + for(int bin_R = 0; bin_RcentreProj().x(),2) + pow(clu->centreProj().y(),2) ); + int bin_R = int( (ROverZ-kROverZMin_) * nBinsRHisto_ / (kROverZMax_-kROverZMin_) ); + int bin_phi = int( (reco::reduceRange(clu->phi())+M_PI) * nBinsPhiHisto_ / (2*M_PI) ); + + histoClusters[{{triggerTools_.zside(clu->detId()), bin_R, bin_phi}}]+=clu->mipPt(); + + } + + return histoClusters; + +} + + + + +HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillSmoothPhiHistoClusters( const Histogram & histoClusters, + const vector & binSums ){ + + Histogram histoSumPhiClusters; //key[0] = z.side(), key[1] = bin_R, key[2] = bin_phi + + for(int z_side : {-1,1}){ + + for(int bin_R = 0; bin_R=int(nBinsPhiHisto_) ) binToSumRight -= nBinsPhiHisto_; + + content += histoClusters.at({{z_side,bin_R,binToSumLeft}}) / pow(2,bin_phi2); // quadratic kernel + content += histoClusters.at({{z_side,bin_R,binToSumRight}}) / pow(2,bin_phi2); // quadratic kernel + + } + + histoSumPhiClusters[{{z_side,bin_R,bin_phi}}] = content/area; + + } + + } + + } + + return histoSumPhiClusters; + +} + + + + + + +HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillSmoothRPhiHistoClusters( const Histogram & histoClusters ){ + + Histogram histoSumRPhiClusters; //key[0] = z.side(), key[1] = bin_R, key[2] = bin_phi + + for(int z_side : {-1,1}){ + + for(int bin_R = 0; bin_R0 ? histoClusters.at({{z_side,bin_R-1,bin_phi}}) : 0 ; + float contentUp = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,bin_phi}}) : 0; + + histoSumRPhiClusters[{{z_side,bin_R,bin_phi}}] = (content + 0.5*contentDown + 0.5*contentUp)/weight; + + } + + } + + } + + return histoSumRPhiClusters; + +} + + + + +std::vector > HGCalHistoSeedingImpl::computeMaxSeeds( const Histogram & histoClusters ){ + + std::vector > seedPositionsEnergy; + + for(int z_side : {-1,1}){ + + for(int bin_R = 0; bin_R histoThreshold_; + if (!isMax) continue; + + float MIPT_S = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,bin_phi}}) : 0; + float MIPT_N = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,bin_phi}}) : 0; + + int binLeft = bin_phi - 1; + if( binLeft<0 ) binLeft += nBinsPhiHisto_; + int binRight = bin_phi + 1; + if( binRight>=int(nBinsPhiHisto_) ) binRight -= nBinsPhiHisto_; + + float MIPT_W = histoClusters.at({{z_side,bin_R,binLeft}}); + float MIPT_E = histoClusters.at({{z_side,bin_R,binRight}}); + float MIPT_NW = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,binLeft}}) : 0; + float MIPT_NE = bin_R>0 ?histoClusters.at({{z_side,bin_R-1,binRight}}) : 0; + float MIPT_SW = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binLeft}}) : 0; + float MIPT_SE = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binRight}}) : 0; + + isMax &= MIPT_seed>=MIPT_S + && MIPT_seed>MIPT_N + && MIPT_seed>=MIPT_E + && MIPT_seed>=MIPT_SE + && MIPT_seed>=MIPT_NE + && MIPT_seed>MIPT_W + && MIPT_seed>MIPT_SW + && MIPT_seed>MIPT_NW; + + if(isMax){ + + float ROverZ_seed = kROverZMin_ + (bin_R+0.5) * (kROverZMax_-kROverZMin_)/nBinsRHisto_; + float phi_seed = -M_PI + (bin_phi+0.5) * 2*M_PI/nBinsPhiHisto_; + float x_seed = ROverZ_seed*cos(phi_seed); + float y_seed = ROverZ_seed*sin(phi_seed); + seedPositionsEnergy.emplace_back( GlobalPoint(x_seed,y_seed,z_side), MIPT_seed); + + } + } + + } + + } + + return seedPositionsEnergy; + +} + + +std::vector > HGCalHistoSeedingImpl::computeInterpolatedMaxSeeds( const Histogram & histoClusters ){ + + + std::vector > seedPositionsEnergy; + + for(int z_side : {-1,1}){ + + for(int bin_R = 0; bin_R0 ? histoClusters.at({{z_side,bin_R-1,bin_phi}}) : 0; + + int binLeft = bin_phi - 1; + if( binLeft<0 ) binLeft += nBinsPhiHisto_; + int binRight = bin_phi + 1; + if( binRight>=int(nBinsPhiHisto_) ) binRight -= nBinsPhiHisto_; + + float MIPT_W = histoClusters.at({{z_side,bin_R,binLeft}}); + float MIPT_E = histoClusters.at({{z_side,bin_R,binRight}}); + + float MIPT_NW = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,binLeft}}) : 0; + float MIPT_NE = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,binRight}}) : 0; + float MIPT_SW = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binLeft}}) : 0; + float MIPT_SE = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binRight}}) : 0; + + float MIPT_pred = neighbour_weights_.at(0) * MIPT_NW + neighbour_weights_.at(1) * MIPT_N + neighbour_weights_.at(2) * MIPT_NE + + neighbour_weights_.at(3) * MIPT_W + neighbour_weights_.at(5) * MIPT_E + neighbour_weights_.at(6) * MIPT_SW + + neighbour_weights_.at(7) * MIPT_S + neighbour_weights_.at(8) * MIPT_SE; + + bool isMax = MIPT_seed>=(MIPT_pred+histoThreshold_); + + if(isMax){ + float ROverZ_seed = kROverZMin_ + (bin_R+0.5) * (kROverZMax_-kROverZMin_)/nBinsRHisto_; + float phi_seed = -M_PI + (bin_phi+0.5) * 2*M_PI/nBinsPhiHisto_; + float x_seed = ROverZ_seed*cos(phi_seed); + float y_seed = ROverZ_seed*sin(phi_seed); + seedPositionsEnergy.emplace_back( GlobalPoint(x_seed,y_seed,z_side), MIPT_seed); + } + + } + + } + + } + + return seedPositionsEnergy; + +} + + +std::vector > HGCalHistoSeedingImpl::computeThresholdSeeds( const Histogram & histoClusters ){ + + + std::vector > seedPositionsEnergy; + + for(int z_side : {-1,1}){ + + for(int bin_R = 0; bin_R histoThreshold_; + + if(isSeed){ + float ROverZ_seed = kROverZMin_ + (bin_R+0.5) * (kROverZMax_-kROverZMin_)/nBinsRHisto_; + float phi_seed = -M_PI + (bin_phi+0.5) * 2*M_PI/nBinsPhiHisto_; + float x_seed = ROverZ_seed*cos(phi_seed); + float y_seed = ROverZ_seed*sin(phi_seed); + seedPositionsEnergy.emplace_back( GlobalPoint(x_seed,y_seed,z_side), MIPT_seed) ; + } + + } + + } + + } + + return seedPositionsEnergy; + +} + + + +std::vector > HGCalHistoSeedingImpl::computeSecondaryMaxSeeds( const Histogram & histoClusters ){ + + std::vector > seedPositionsEnergy; + + std::map, bool> primarySeedPositions; + std::map, bool> secondarySeedPositions; + std::map, bool> vetoPositions; + + //Search for primary seeds + for(int z_side : {-1,1}){ + + for(int bin_R = 0; bin_R histoThreshold_; + + if (!isMax) continue; + + float MIPT_S = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,bin_phi}}) : 0; + float MIPT_N = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,bin_phi}}) : 0; + + int binLeft = bin_phi - 1; + if( binLeft<0 ) binLeft += nBinsPhiHisto_; + int binRight = bin_phi + 1; + if( binRight>=int(nBinsPhiHisto_) ) binRight -= nBinsPhiHisto_; + + float MIPT_W = histoClusters.at({{z_side,bin_R,binLeft}}); + float MIPT_E = histoClusters.at({{z_side,bin_R,binRight}}); + float MIPT_NW = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,binLeft}}) : 0; + float MIPT_NE = bin_R>0 ?histoClusters.at({{z_side,bin_R-1,binRight}}) : 0; + float MIPT_SW = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binLeft}}) : 0; + float MIPT_SE = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binRight}}) : 0; + + isMax &= MIPT_seed>=MIPT_S + && MIPT_seed>MIPT_N + && MIPT_seed>=MIPT_E + && MIPT_seed>=MIPT_SE + && MIPT_seed>=MIPT_NE + && MIPT_seed>MIPT_W + && MIPT_seed>MIPT_SW + && MIPT_seed>MIPT_NW; + + if(isMax){ + + + float ROverZ_seed = kROverZMin_ + (bin_R+0.5) * (kROverZMax_-kROverZMin_)/nBinsRHisto_; + float phi_seed = -M_PI + (bin_phi+0.5) * 2*M_PI/nBinsPhiHisto_; + float x_seed = ROverZ_seed*cos(phi_seed); + float y_seed = ROverZ_seed*sin(phi_seed); + + seedPositionsEnergy.emplace_back( GlobalPoint(x_seed,y_seed,z_side), MIPT_seed ); + primarySeedPositions[std::make_tuple(bin_R,bin_phi,z_side)] = true; + + vetoPositions[std::make_tuple(bin_R,binLeft,z_side)] = true; + vetoPositions[std::make_tuple(bin_R,binRight,z_side)] = true; + if ( bin_R>0 ) { + vetoPositions[std::make_tuple(bin_R-1,bin_phi,z_side)] = true; + vetoPositions[std::make_tuple(bin_R-1,binRight,z_side)] = true; + vetoPositions[std::make_tuple(bin_R-1,binLeft,z_side)] = true; + } + if ( bin_R<(int(nBinsRHisto_)-1) ) { + vetoPositions[std::make_tuple(bin_R+1,bin_phi,z_side)] = true; + vetoPositions[std::make_tuple(bin_R+1,binRight,z_side)] = true; + vetoPositions[std::make_tuple(bin_R+1,binLeft,z_side)] = true; + } + + } + + } + + } + + } + + + //Search for secondary seeds + + for(int z_side : {-1,1}){ + + for(int bin_R = 0; bin_R histoThreshold_; + + float MIPT_S = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,bin_phi}}) : 0; + float MIPT_N = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,bin_phi}}) : 0; + + int binLeft = bin_phi - 1; + if( binLeft<0 ) binLeft += nBinsPhiHisto_; + int binRight = bin_phi + 1; + if( binRight>=int(nBinsPhiHisto_) ) binRight -= nBinsPhiHisto_; + + float MIPT_W = histoClusters.at({{z_side,bin_R,binLeft}}); + float MIPT_E = histoClusters.at({{z_side,bin_R,binRight}}); + float MIPT_NW = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,binLeft}}) : 0; + float MIPT_NE = bin_R>0 ?histoClusters.at({{z_side,bin_R-1,binRight}}) : 0; + float MIPT_SW = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binLeft}}) : 0; + float MIPT_SE = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binRight}}) : 0; + + + isMax &= ( vetoPositions[std::make_tuple(bin_R+1,bin_phi,z_side)] or MIPT_seed>=MIPT_S ) + && ( vetoPositions[std::make_tuple(bin_R-1,bin_phi,z_side)] or MIPT_seed>MIPT_N ) + && ( vetoPositions[std::make_tuple(bin_R,binRight,z_side)] or MIPT_seed>=MIPT_E ) + && ( vetoPositions[std::make_tuple(bin_R+1,binRight,z_side)] or MIPT_seed>=MIPT_SE ) + && ( vetoPositions[std::make_tuple(bin_R-1,binRight,z_side)] or MIPT_seed>=MIPT_NE ) + && ( vetoPositions[std::make_tuple(bin_R,binLeft,z_side)] or MIPT_seed>MIPT_W ) + && ( vetoPositions[std::make_tuple(bin_R+1,binLeft,z_side)] or MIPT_seed>MIPT_SW ) + && ( vetoPositions[std::make_tuple(bin_R-1,binLeft,z_side)] or MIPT_seed>MIPT_NW ); + + + + if(isMax){ + float ROverZ_seed = kROverZMin_ + (bin_R+0.5) * (kROverZMax_-kROverZMin_)/nBinsRHisto_; + float phi_seed = -M_PI + (bin_phi+0.5) * 2*M_PI/nBinsPhiHisto_; + float x_seed = ROverZ_seed*cos(phi_seed); + float y_seed = ROverZ_seed*sin(phi_seed); + seedPositionsEnergy.emplace_back( GlobalPoint(x_seed,y_seed,z_side), MIPT_seed ); + secondarySeedPositions[std::make_tuple(bin_R,bin_phi,z_side)] = true; + } + + } + + } + + } + + return seedPositionsEnergy; + +} + + +void HGCalHistoSeedingImpl::findHistoSeeds( const std::vector> & clustersPtrs, + std::vector > & seedPositionsEnergy) +{ + + /* put clusters into an r/z x phi histogram */ + Histogram histoCluster = fillHistoClusters(clustersPtrs); //key[0] = z.side(), key[1] = bin_R, key[2] = bin_phi, content = MIPTs summed along depth + + /* smoothen along the phi direction + normalize each bin to same area */ + Histogram smoothPhiHistoCluster = fillSmoothPhiHistoClusters(histoCluster,binsSumsHisto_); + + /* smoothen along the r/z direction */ + Histogram smoothRPhiHistoCluster = fillSmoothRPhiHistoClusters(histoCluster); + + /* seeds determined with local maximum criteria */ + if (seedingType_ == HistoMaxC3d) seedPositionsEnergy = computeMaxSeeds(smoothRPhiHistoCluster); + else if(seedingType_ == HistoThresholdC3d) seedPositionsEnergy = computeThresholdSeeds(smoothRPhiHistoCluster); + else if(seedingType_ == HistoInterpolatedMaxC3d) seedPositionsEnergy = computeInterpolatedMaxSeeds(smoothRPhiHistoCluster); + else if(seedingType_ == HistoSecondaryMaxC3d) seedPositionsEnergy = computeSecondaryMaxSeeds(smoothRPhiHistoCluster); + +} + diff --git a/L1Trigger/L1THGCal/src/backend/HGCalMulticlusteringHistoImpl.cc b/L1Trigger/L1THGCal/src/backend/HGCalMulticlusteringHistoImpl.cc deleted file mode 100644 index 308f124fe390e..0000000000000 --- a/L1Trigger/L1THGCal/src/backend/HGCalMulticlusteringHistoImpl.cc +++ /dev/null @@ -1,446 +0,0 @@ -#include "L1Trigger/L1THGCal/interface/backend/HGCalMulticlusteringHistoImpl.h" -#include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h" -#include "DataFormats/Math/interface/deltaR.h" -#include "DataFormats/Math/interface/deltaPhi.h" - -HGCalMulticlusteringHistoImpl::HGCalMulticlusteringHistoImpl( const edm::ParameterSet& conf ) : - dr_(conf.getParameter("dR_multicluster")), - ptC3dThreshold_(conf.getParameter("minPt_multicluster")), - multiclusterAlgoType_(conf.getParameter("type_multicluster")), - nBinsRHisto_(conf.getParameter("nBins_R_histo_multicluster")), - nBinsPhiHisto_(conf.getParameter("nBins_Phi_histo_multicluster")), - binsSumsHisto_(conf.getParameter< std::vector >("binSumsHisto")), - histoThreshold_(conf.getParameter("threshold_histo_multicluster")), - neighbour_weights_(conf.getParameter< std::vector >("neighbour_weights")) -{ - - if(multiclusterAlgoType_=="HistoMaxC3d"){ - multiclusteringAlgoType_ = HistoMaxC3d; - }else if(multiclusterAlgoType_=="HistoThresholdC3d"){ - multiclusteringAlgoType_ = HistoThresholdC3d; - }else if(multiclusterAlgoType_=="HistoInterpolatedMaxC3d"){ - multiclusteringAlgoType_ = HistoInterpolatedMaxC3d; - }else { - throw cms::Exception("HGCTriggerParameterError") - << "Unknown Multiclustering type '" << multiclusterAlgoType_; - } - - edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster dR: " << dr_ - <<"\nMulticluster minimum transverse-momentum: " << ptC3dThreshold_ - <<"\nMulticluster number of R-bins for the histo algorithm: " << nBinsRHisto_ - <<"\nMulticluster number of Phi-bins for the histo algorithm: " << nBinsPhiHisto_ - <<"\nMulticluster MIPT threshold for histo threshold algorithm: " << histoThreshold_ - <<"\nMulticluster type of multiclustering algortihm: " << multiclusterAlgoType_; - - id_ = std::unique_ptr{ HGCalTriggerClusterIdentificationFactory::get()->create("HGCalTriggerClusterIdentificationBDT") }; - id_->initialize(conf.getParameter("EGIdentification")); - if(multiclusterAlgoType_.find("Histo")!=std::string::npos && nBinsRHisto_!=binsSumsHisto_.size()){ - throw cms::Exception("Inconsistent bin size") << "Inconsistent nBins_R_histo_multicluster ( " << nBinsRHisto_ << " ) and binSumsHisto ( " << binsSumsHisto_.size() << " ) size in HGCalMulticlustering\n"; - } - - if(neighbour_weights_.size()!=neighbour_weights_size_) { - throw cms::Exception("Inconsistent vector size" ) << "Inconsistent size of neighbour weights vector in HGCalMulticlustering ( " << neighbour_weights_.size() << " ). Should be " << neighbour_weights_size_ << "\n" ; - } - -} - - - - - -float HGCalMulticlusteringHistoImpl::dR( const l1t::HGCalCluster & clu, - const GlobalPoint & seed) const -{ - - Basic3DVector seed_3dv( seed ); - GlobalPoint seed_proj( seed_3dv / seed.z() ); - return (seed_proj - clu.centreProj() ).mag(); - -} - - - - - -HGCalMulticlusteringHistoImpl::Histogram HGCalMulticlusteringHistoImpl::fillHistoClusters( const std::vector> & clustersPtrs ){ - - - Histogram histoClusters; //key[0] = z.side(), key[1] = bin_R, key[2] = bin_phi - - for(int z_side : {-1,1}){ - - for(int bin_R = 0; bin_RcentreProj().x(),2) + pow(clu->centreProj().y(),2) ); - int bin_R = int( (ROverZ-kROverZMin_) * nBinsRHisto_ / (kROverZMax_-kROverZMin_) ); - int bin_phi = int( (reco::reduceRange(clu->phi())+M_PI) * nBinsPhiHisto_ / (2*M_PI) ); - - histoClusters[{{triggerTools_.zside(clu->detId()), bin_R, bin_phi}}]+=clu->mipPt(); - - } - - return histoClusters; - -} - - - - -HGCalMulticlusteringHistoImpl::Histogram HGCalMulticlusteringHistoImpl::fillSmoothPhiHistoClusters( const Histogram & histoClusters, - const vector & binSums ){ - - Histogram histoSumPhiClusters; //key[0] = z.side(), key[1] = bin_R, key[2] = bin_phi - - for(int z_side : {-1,1}){ - - for(int bin_R = 0; bin_R=int(nBinsPhiHisto_) ) binToSumRight -= nBinsPhiHisto_; - - content += histoClusters.at({{z_side,bin_R,binToSumLeft}}) / pow(2,bin_phi2); // quadratic kernel - content += histoClusters.at({{z_side,bin_R,binToSumRight}}) / pow(2,bin_phi2); // quadratic kernel - - } - - histoSumPhiClusters[{{z_side,bin_R,bin_phi}}] = content/area; - - } - - } - - } - - return histoSumPhiClusters; - -} - - - - - - -HGCalMulticlusteringHistoImpl::Histogram HGCalMulticlusteringHistoImpl::fillSmoothRPhiHistoClusters( const Histogram & histoClusters ){ - - Histogram histoSumRPhiClusters; //key[0] = z.side(), key[1] = bin_R, key[2] = bin_phi - - for(int z_side : {-1,1}){ - - for(int bin_R = 0; bin_R0 ? histoClusters.at({{z_side,bin_R-1,bin_phi}}) : 0 ; - float contentUp = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,bin_phi}}) : 0; - - histoSumRPhiClusters[{{z_side,bin_R,bin_phi}}] = (content + 0.5*contentDown + 0.5*contentUp)/weight; - - } - - } - - } - - return histoSumRPhiClusters; - -} - - - - - -std::vector HGCalMulticlusteringHistoImpl::computeMaxSeeds( const Histogram & histoClusters ){ - - std::vector seedPositions; - - for(int z_side : {-1,1}){ - - for(int bin_R = 0; bin_R0; - - float MIPT_S = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,bin_phi}}) : 0; - float MIPT_N = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,bin_phi}}) : 0; - - int binLeft = bin_phi - 1; - if( binLeft<0 ) binLeft += nBinsPhiHisto_; - int binRight = bin_phi + 1; - if( binRight>=int(nBinsPhiHisto_) ) binRight -= nBinsPhiHisto_; - - float MIPT_W = histoClusters.at({{z_side,bin_R,binLeft}}); - float MIPT_E = histoClusters.at({{z_side,bin_R,binRight}}); - float MIPT_NW = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,binLeft}}) : 0; - float MIPT_NE = bin_R>0 ?histoClusters.at({{z_side,bin_R-1,binRight}}) : 0; - float MIPT_SW = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binLeft}}) : 0; - float MIPT_SE = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binRight}}) : 0; - - isMax &= MIPT_seed>=MIPT_S; - isMax &= MIPT_seed>MIPT_N; - isMax &= MIPT_seed>=MIPT_E; - isMax &= MIPT_seed>=MIPT_SE; - isMax &= MIPT_seed>=MIPT_NE; - isMax &= MIPT_seed>MIPT_W; - isMax &= MIPT_seed>MIPT_SW; - isMax &= MIPT_seed>MIPT_NW; - - - - if(isMax){ - - - float ROverZ_seed = kROverZMin_ + (bin_R+0.5) * (kROverZMax_-kROverZMin_)/nBinsRHisto_; - float phi_seed = -M_PI + (bin_phi+0.5) * 2*M_PI/nBinsPhiHisto_; - float x_seed = ROverZ_seed*cos(phi_seed); - float y_seed = ROverZ_seed*sin(phi_seed); - seedPositions.emplace_back(x_seed,y_seed,z_side); - } - - } - - } - - } - - return seedPositions; - -} - - -std::vector HGCalMulticlusteringHistoImpl::computeInterpolatedMaxSeeds( const Histogram & histoClusters ){ - - std::vector seedPositions; - - for(int z_side : {-1,1}){ - - for(int bin_R = 0; bin_R0; - - float MIPT_S = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,bin_phi}}) : 0; - float MIPT_N = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,bin_phi}}) : 0; - - int binLeft = bin_phi - 1; - if( binLeft<0 ) binLeft += nBinsPhiHisto_; - int binRight = bin_phi + 1; - if( binRight>=int(nBinsPhiHisto_) ) binRight -= nBinsPhiHisto_; - - float MIPT_W = histoClusters.at({{z_side,bin_R,binLeft}}); - float MIPT_E = histoClusters.at({{z_side,bin_R,binRight}}); - - float MIPT_NW = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,binLeft}}) : 0; - float MIPT_NE = bin_R>0 ? histoClusters.at({{z_side,bin_R-1,binRight}}) : 0; - float MIPT_SW = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binLeft}}) : 0; - float MIPT_SE = bin_R<(int(nBinsRHisto_)-1) ? histoClusters.at({{z_side,bin_R+1,binRight}}) : 0; - - float MIPT_pred = neighbour_weights_.at(0) * MIPT_NW + neighbour_weights_.at(1) * MIPT_N + neighbour_weights_.at(2) * MIPT_NE - + neighbour_weights_.at(3) * MIPT_W + neighbour_weights_.at(5) * MIPT_E + neighbour_weights_.at(6) * MIPT_SW - + neighbour_weights_.at(7) * MIPT_S + neighbour_weights_.at(8) * MIPT_SE; - - isMax &= MIPT_seed>=MIPT_pred; - - if(isMax){ - float ROverZ_seed = kROverZMin_ + (bin_R+0.5) * (kROverZMax_-kROverZMin_)/nBinsRHisto_; - float phi_seed = -M_PI + (bin_phi+0.5) * 2*M_PI/nBinsPhiHisto_; - float x_seed = ROverZ_seed*cos(phi_seed); - float y_seed = ROverZ_seed*sin(phi_seed); - seedPositions.emplace_back(x_seed,y_seed,z_side); - } - - } - - } - - } - - return seedPositions; - -} - - -std::vector HGCalMulticlusteringHistoImpl::computeThresholdSeeds( const Histogram & histoClusters ){ - - std::vector seedPositions; - - for(int z_side : {-1,1}){ - - for(int bin_R = 0; bin_R histoThreshold_; - - if(isSeed){ - float ROverZ_seed = kROverZMin_ + (bin_R+0.5) * (kROverZMax_-kROverZMin_)/nBinsRHisto_; - float phi_seed = -M_PI + (bin_phi+0.5) * 2*M_PI/nBinsPhiHisto_; - float x_seed = ROverZ_seed*cos(phi_seed); - float y_seed = ROverZ_seed*sin(phi_seed); - seedPositions.emplace_back(x_seed,y_seed,z_side); - } - - } - - } - - } - - return seedPositions; - -} - - - -std::vector HGCalMulticlusteringHistoImpl::clusterSeedMulticluster(const std::vector> & clustersPtrs, - const std::vector & seeds){ - - - std::map mapSeedMulticluster; - std::vector multiclustersTmp; - - for(auto & clu : clustersPtrs){ - - - int z_side = triggerTools_.zside(clu->detId()); - - double minDist = dr_; - int targetSeed = -1; - - for( unsigned int iseed=0; iseeddR(*clu, seeds[iseed]); - - if(d> & clustersPtrs, - l1t::HGCalMulticlusterBxCollection & multiclusters, - const HGCalTriggerGeometryBase & triggerGeometry) -{ - - /* put clusters into an r/z x phi histogram */ - Histogram histoCluster = fillHistoClusters(clustersPtrs); //key[0] = z.side(), key[1] = bin_R, key[2] = bin_phi, content = MIPTs summed along depth - - /* smoothen along the phi direction + normalize each bin to same area */ - Histogram smoothPhiHistoCluster = fillSmoothPhiHistoClusters(histoCluster,binsSumsHisto_); - - /* smoothen along the r/z direction */ - Histogram smoothRPhiHistoCluster = fillSmoothRPhiHistoClusters(histoCluster); - - /* seeds determined with local maximum criteria */ - std::vector seedPositions; - if (multiclusteringAlgoType_ == HistoMaxC3d) seedPositions = computeMaxSeeds(smoothRPhiHistoCluster); - else if(multiclusteringAlgoType_ == HistoThresholdC3d) seedPositions = computeThresholdSeeds(smoothRPhiHistoCluster); - else if(multiclusteringAlgoType_ == HistoInterpolatedMaxC3d) seedPositions = computeInterpolatedMaxSeeds(smoothRPhiHistoCluster); - /* clusterize clusters around seeds */ - std::vector multiclustersTmp = clusterSeedMulticluster(clustersPtrs,seedPositions); - - /* making the collection of multiclusters */ - finalizeClusters(multiclustersTmp, multiclusters, triggerGeometry); - -} - - - - - - - -void -HGCalMulticlusteringHistoImpl:: -finalizeClusters(std::vector& multiclusters_in, - l1t::HGCalMulticlusterBxCollection& multiclusters_out, - const HGCalTriggerGeometryBase& triggerGeometry) { - for(auto& multicluster : multiclusters_in) { - // compute the eta, phi from its barycenter - // + pT as scalar sum of pT of constituents - double sumPt=0.; - const std::unordered_map>& clusters = multicluster.constituents(); - for(const auto& id_cluster : clusters) sumPt += id_cluster.second->pt(); - - math::PtEtaPhiMLorentzVector multiclusterP4( sumPt, - multicluster.centre().eta(), - multicluster.centre().phi(), - 0. ); - multicluster.setP4( multiclusterP4 ); - - if( multicluster.pt() > ptC3dThreshold_ ){ - //compute shower shapes - multicluster.showerLength(shape_.showerLength(multicluster)); - multicluster.coreShowerLength(shape_.coreShowerLength(multicluster, triggerGeometry)); - multicluster.firstLayer(shape_.firstLayer(multicluster)); - multicluster.maxLayer(shape_.maxLayer(multicluster)); - multicluster.sigmaEtaEtaTot(shape_.sigmaEtaEtaTot(multicluster)); - multicluster.sigmaEtaEtaMax(shape_.sigmaEtaEtaMax(multicluster)); - multicluster.sigmaPhiPhiTot(shape_.sigmaPhiPhiTot(multicluster)); - multicluster.sigmaPhiPhiMax(shape_.sigmaPhiPhiMax(multicluster)); - multicluster.sigmaZZ(shape_.sigmaZZ(multicluster)); - multicluster.sigmaRRTot(shape_.sigmaRRTot(multicluster)); - multicluster.sigmaRRMax(shape_.sigmaRRMax(multicluster)); - multicluster.sigmaRRMean(shape_.sigmaRRMean(multicluster)); - multicluster.eMax(shape_.eMax(multicluster)); - // fill quality flag - multicluster.setHwQual(id_->decision(multicluster)); - // fill H/E - multicluster.saveHOverE(); - - multiclusters_out.push_back( 0, multicluster); - } - } -} diff --git a/L1Trigger/L1THGCal/src/backend/HGCalMulticlusteringImpl.cc b/L1Trigger/L1THGCal/src/backend/HGCalMulticlusteringImpl.cc index 6d912c26d4975..c32f8c4f91889 100644 --- a/L1Trigger/L1THGCal/src/backend/HGCalMulticlusteringImpl.cc +++ b/L1Trigger/L1THGCal/src/backend/HGCalMulticlusteringImpl.cc @@ -190,19 +190,7 @@ finalizeClusters(std::vector& multiclusters_in, if( multicluster.pt() > ptC3dThreshold_ ){ //compute shower shapes - multicluster.showerLength(shape_.showerLength(multicluster)); - multicluster.coreShowerLength(shape_.coreShowerLength(multicluster, triggerGeometry)); - multicluster.firstLayer(shape_.firstLayer(multicluster)); - multicluster.maxLayer(shape_.maxLayer(multicluster)); - multicluster.sigmaEtaEtaTot(shape_.sigmaEtaEtaTot(multicluster)); - multicluster.sigmaEtaEtaMax(shape_.sigmaEtaEtaMax(multicluster)); - multicluster.sigmaPhiPhiTot(shape_.sigmaPhiPhiTot(multicluster)); - multicluster.sigmaPhiPhiMax(shape_.sigmaPhiPhiMax(multicluster)); - multicluster.sigmaZZ(shape_.sigmaZZ(multicluster)); - multicluster.sigmaRRTot(shape_.sigmaRRTot(multicluster)); - multicluster.sigmaRRMax(shape_.sigmaRRMax(multicluster)); - multicluster.sigmaRRMean(shape_.sigmaRRMean(multicluster)); - multicluster.eMax(shape_.eMax(multicluster)); + shape_.fillShapes(multicluster, triggerGeometry); // fill quality flag multicluster.setHwQual(id_->decision(multicluster)); // fill H/E diff --git a/L1Trigger/L1THGCal/src/backend/HGCalShowerShape.cc b/L1Trigger/L1THGCal/src/backend/HGCalShowerShape.cc index ef77603b47a2d..2fe43fbdecc16 100644 --- a/L1Trigger/L1THGCal/src/backend/HGCalShowerShape.cc +++ b/L1Trigger/L1THGCal/src/backend/HGCalShowerShape.cc @@ -486,3 +486,25 @@ float HGCalShowerShape::sigmaRRTot(const l1t::HGCalCluster& c2d) const { return Srr; } + + + +void +HGCalShowerShape:: +fillShapes(l1t::HGCalMulticluster& c3d, const HGCalTriggerGeometryBase& triggerGeometry) const +{ + c3d.showerLength(showerLength(c3d)); + c3d.coreShowerLength(coreShowerLength(c3d, triggerGeometry)); + c3d.firstLayer(firstLayer(c3d)); + c3d.maxLayer(maxLayer(c3d)); + c3d.sigmaEtaEtaTot(sigmaEtaEtaTot(c3d)); + c3d.sigmaEtaEtaMax(sigmaEtaEtaMax(c3d)); + c3d.sigmaPhiPhiTot(sigmaPhiPhiTot(c3d)); + c3d.sigmaPhiPhiMax(sigmaPhiPhiMax(c3d)); + c3d.sigmaZZ(sigmaZZ(c3d)); + c3d.sigmaRRTot(sigmaRRTot(c3d)); + c3d.sigmaRRMax(sigmaRRMax(c3d)); + c3d.sigmaRRMean(sigmaRRMean(c3d)); + c3d.eMax(eMax(c3d)); +} + diff --git a/L1Trigger/L1THGCal/src/concentrator/HGCalConcentratorSuperTriggerCellImpl.cc b/L1Trigger/L1THGCal/src/concentrator/HGCalConcentratorSuperTriggerCellImpl.cc index beb52a1c73742..c5f7e107f7225 100644 --- a/L1Trigger/L1THGCal/src/concentrator/HGCalConcentratorSuperTriggerCellImpl.cc +++ b/L1Trigger/L1THGCal/src/concentrator/HGCalConcentratorSuperTriggerCellImpl.cc @@ -1,21 +1,107 @@ #include "L1Trigger/L1THGCal/interface/concentrator/HGCalConcentratorSuperTriggerCellImpl.h" +#include "DataFormats/ForwardDetId/interface/HGCalTriggerDetId.h" #include HGCalConcentratorSuperTriggerCellImpl:: -HGCalConcentratorSuperTriggerCellImpl(const edm::ParameterSet& conf){} +HGCalConcentratorSuperTriggerCellImpl(const edm::ParameterSet& conf) + : stcSize_(conf.getParameter< std::vector >("stcSize")) +{ + if ( stcSize_.size() != kNLayers_ ){ + throw cms::Exception("HGCTriggerParameterError") + << "Inconsistent size of super trigger cell size vector" << stcSize_.size() ; + } + for(auto stc : stcSize_) { + if ( stc!=kSTCsizeFine_ && stc!=kSTCsizeCoarse_ ){ + throw cms::Exception("HGCTriggerParameterError") + << "Super Trigger Cell should be of size "<< + kSTCsizeFine_ << " or " << kSTCsizeCoarse_; + } + } + +} + +std::map +HGCalConcentratorSuperTriggerCellImpl::kSplit_ = { + {kSTCsizeFine_, kSplit_v8_Fine_}, + {kSTCsizeCoarse_, kSplit_v8_Coarse_} +}; int HGCalConcentratorSuperTriggerCellImpl::getSuperTriggerCellId(int detid) const { - // FIXME: won't work in the V9 geometry - HGCalDetId TC_id(detid); - if(TC_id.subdetId()==HGCHEB) { - return TC_id.cell(); //scintillator - } else { - int TC_wafer = TC_id.wafer(); - int TC_12th = ( TC_id.cell() & kSplit_ ); - return TC_wafer< + + + + + + + + diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h b/L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h similarity index 85% rename from L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h rename to L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h index 3b661154c0c38..0b9d678bbdf1b 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h +++ b/L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h @@ -1,5 +1,5 @@ -#ifndef __L1Trigger_L1THGCal_HGCalTriggerNtupleBase_h__ -#define __L1Trigger_L1THGCal_HGCalTriggerNtupleBase_h__ +#ifndef __L1Trigger_L1THGCalUtilities_HGCalTriggerNtupleBase_h__ +#define __L1Trigger_L1THGCalUtilities_HGCalTriggerNtupleBase_h__ #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/EventSetup.h" diff --git a/L1Trigger/L1THGCalUtilities/plugins/BuildFile.xml b/L1Trigger/L1THGCalUtilities/plugins/BuildFile.xml new file mode 100644 index 0000000000000..77ee5b8ff650d --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/plugins/BuildFile.xml @@ -0,0 +1,25 @@ + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleEvent.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleEvent.cc similarity index 93% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleEvent.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleEvent.cc index aa94870fd1f28..a2bcbc042323f 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleEvent.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleEvent.cc @@ -1,4 +1,4 @@ -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" class HGCalTriggerNtupleEvent : public HGCalTriggerNtupleBase { diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGen.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGen.cc similarity index 99% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGen.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGen.cc index fc9be6adecca4..2a2914e64704c 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGen.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGen.cc @@ -4,7 +4,7 @@ #include "DataFormats/HcalDetId/interface/HcalDetId.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" #include "MagneticField/Engine/interface/MagneticField.h" diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenJet.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGenJet.cc similarity index 96% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenJet.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGenJet.cc index ba0e3ade042d0..3a2915caae244 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenJet.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGenJet.cc @@ -2,7 +2,7 @@ #include "DataFormats/JetReco/interface/GenJet.h" #include "DataFormats/JetReco/interface/GenJetCollection.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenTau.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGenTau.cc similarity index 99% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenTau.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGenTau.cc index eb15b1a32d6bd..927f0019f4d23 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenTau.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGenTau.cc @@ -1,7 +1,7 @@ #include #include "DataFormats/Math/interface/LorentzVector.h" #include "DataFormats/HepMCCandidate/interface/GenParticle.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" #include "DataFormats/Math/interface/LorentzVector.h" typedef math::XYZTLorentzVector LorentzVector; diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc similarity index 98% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc index 9e688cd18c433..ddfc840807b1a 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc @@ -4,7 +4,7 @@ #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc similarity index 99% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc index 29c10d2d1ff65..7300a61f42bee 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc @@ -2,7 +2,7 @@ #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "DataFormats/HcalDetId/interface/HcalDetId.h" #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" #include "FWCore/Framework/interface/ESHandle.h" diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc similarity index 98% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc index 4b647db46cdbc..8c01b9624c0bf 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc @@ -1,7 +1,7 @@ #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" #include "L1Trigger/L1THGCal/interface/backend/HGCalTriggerClusterIdentificationBase.h" diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCPanels.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCPanels.cc similarity index 98% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCPanels.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCPanels.cc index 7ed8e316158e4..92976057f49b9 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCPanels.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCPanels.cc @@ -3,7 +3,7 @@ #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc similarity index 99% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc index e599031d74e47..623726317ac51 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc @@ -8,7 +8,7 @@ #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleManager.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleManager.cc similarity index 96% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleManager.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleManager.cc index f9bbdda7fa263..9bfc16fa5e265 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleManager.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleManager.cc @@ -6,7 +6,7 @@ #include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/ServiceRegistry/interface/Service.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" class HGCalTriggerNtupleManager : public edm::EDAnalyzer { diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleTowers.cc b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleTowers.cc similarity index 97% rename from L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleTowers.cc rename to L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleTowers.cc index 761983d3a57b5..605093aa99da2 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleTowers.cc +++ b/L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleTowers.cc @@ -1,7 +1,7 @@ #include "DataFormats/L1THGCal/interface/HGCalTower.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" diff --git a/L1Trigger/L1THGCalUtilities/plugins/selectors/HGC3DClusterSimpleSelector.cc b/L1Trigger/L1THGCalUtilities/plugins/selectors/HGC3DClusterSimpleSelector.cc new file mode 100644 index 0000000000000..7be9149f93775 --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/plugins/selectors/HGC3DClusterSimpleSelector.cc @@ -0,0 +1,52 @@ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" +#include "CommonTools/Utils/interface/StringCutObjectSelector.h" + +namespace l1t { + class HGC3DClusterSimpleSelector : public edm::global::EDProducer<> { + public: + explicit HGC3DClusterSimpleSelector(const edm::ParameterSet&) ; + ~HGC3DClusterSimpleSelector() override {} + void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; + + private: + const edm::EDGetTokenT src_; + const StringCutObjectSelector cut_; + + }; // class +} // namespace + +l1t::HGC3DClusterSimpleSelector::HGC3DClusterSimpleSelector(const edm::ParameterSet & iConfig) : + src_(consumes(iConfig.getParameter("src"))), + cut_(iConfig.getParameter("cut")) +{ + produces(); +} + + +void +l1t::HGC3DClusterSimpleSelector::produce(edm::StreamID, edm::Event & iEvent, const edm::EventSetup &) const +{ + std::unique_ptr out = std::make_unique(); + + edm::Handle multiclusters; + iEvent.getByToken(src_, multiclusters); + + for (int bx = multiclusters->getFirstBX(); bx <= multiclusters->getLastBX(); ++bx) { + for(auto it = multiclusters->begin(bx), ed = multiclusters->end(bx); it != ed; ++it) { + const auto & c = *it; + if (cut_(c)) { + out->push_back(bx, c); + } + } + } + + iEvent.put(std::move(out)); +} +using l1t::HGC3DClusterSimpleSelector; +DEFINE_FWK_MODULE(HGC3DClusterSimpleSelector); diff --git a/L1Trigger/L1THGCalUtilities/plugins/selectors/HGC3DClusterTMVASelector.cc b/L1Trigger/L1THGCalUtilities/plugins/selectors/HGC3DClusterTMVASelector.cc new file mode 100644 index 0000000000000..c049c3935e026 --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/plugins/selectors/HGC3DClusterTMVASelector.cc @@ -0,0 +1,100 @@ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" +#include "CommonTools/Utils/interface/StringCutObjectSelector.h" +#include "CommonTools/Utils/interface/StringObjectFunction.h" +#include "CommonTools/MVAUtils/interface/TMVAZipReader.h" + +#include "TMVA/Factory.h" +#include "TMVA/Reader.h" + +namespace l1t { + class HGC3DClusterTMVASelector : public edm::stream::EDProducer<> { + public: + explicit HGC3DClusterTMVASelector(const edm::ParameterSet&) ; + ~HGC3DClusterTMVASelector() override {} + + private: + class Var { + public: + Var(const std::string & name, const std::string & expr) : + name_(name), expr_(expr) {} + void declare(TMVA::Reader & r) { r.AddVariable(name_, &val_); } + void fill(const l1t::HGCalMulticluster & c) { val_ = expr_(c); } + private: + std::string name_; + StringObjectFunction expr_; + float val_; + }; + + edm::EDGetTokenT src_; + StringCutObjectSelector preselection_; + std::vector variables_; + std::string method_, weightsFile_; + std::unique_ptr reader_; + StringObjectFunction wp_; + + + void produce(edm::Event&, const edm::EventSetup&) override; + + }; // class +} // namespace + +l1t::HGC3DClusterTMVASelector::HGC3DClusterTMVASelector(const edm::ParameterSet & iConfig) : + src_(consumes(iConfig.getParameter("src"))), + preselection_(iConfig.getParameter("preselection")), + method_(iConfig.getParameter("method")), + weightsFile_(iConfig.getParameter("weightsFile")), + reader_(new TMVA::Reader()), + wp_(iConfig.getParameter("wp")) +{ + // first create all the variables + for (const auto & psvar : iConfig.getParameter>("variables")) { + variables_.emplace_back(psvar.getParameter("name"), psvar.getParameter("value")); + } + // then declare them + for (auto & var : variables_) var.declare(*reader_); + // then read the weights + if (weightsFile_[0] != '/' && weightsFile_[0] != '.') { + weightsFile_ = edm::FileInPath(weightsFile_).fullPath(); + } + reco::details::loadTMVAWeights(&*reader_, method_, weightsFile_); + // finally, declare outputs + produces(); + produces("fail"); +} + + +void +l1t::HGC3DClusterTMVASelector::produce(edm::Event & iEvent, const edm::EventSetup &) +{ + std::unique_ptr out = std::make_unique(); + std::unique_ptr fail = std::make_unique(); + + edm::Handle multiclusters; + iEvent.getByToken(src_, multiclusters); + + for (int bx = multiclusters->getFirstBX(); bx <= multiclusters->getLastBX(); ++bx) { + for(auto it = multiclusters->begin(bx), ed = multiclusters->end(bx); it != ed; ++it) { + const auto & c = *it; + if (preselection_(c)) { + for (auto & var : variables_) var.fill(c); + float mvaOut = reader_->EvaluateMVA(method_); + if (mvaOut > wp_(c)) { + out->push_back(bx, c); + } else { + fail->push_back(bx, c); + } + } + } + } + + iEvent.put(std::move(out)); + iEvent.put(std::move(fail), "fail"); +} +using l1t::HGC3DClusterTMVASelector; +DEFINE_FWK_MODULE(HGC3DClusterTMVASelector); diff --git a/L1Trigger/L1THGCalUtilities/python/clustering2d.py b/L1Trigger/L1THGCalUtilities/python/clustering2d.py new file mode 100644 index 0000000000000..c03050115a24f --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/python/clustering2d.py @@ -0,0 +1,51 @@ +import FWCore.ParameterSet.Config as cms + + +def create_distance(process, inputs, + distance=6.,# cm + seed_threshold=5.,# MipT + cluster_threshold=2.# MipT + ): + producer = process.hgcalBackEndLayer1Producer.clone() + producer.ProcessorParameters.C2d_parameters.seeding_threshold_silicon = cms.double(seed_threshold) + producer.ProcessorParameters.C2d_parameters.seeding_threshold_scintillator = cms.double(seed_threshold) + producer.ProcessorParameters.C2d_parameters.clustering_threshold_silicon = cms.double(cluster_threshold) + producer.ProcessorParameters.C2d_parameters.clustering_threshold_scintillator = cms.double(cluster_threshold) + producer.ProcessorParameters.C2d_parameters.dR_cluster = cms.double(distance) + producer.ProcessorParameters.C2d_parameters.clusterType = cms.string('dRC2d') + producer.InputTriggerCells = cms.InputTag('{}:HGCalConcentratorProcessorSelection'.format(inputs)) + return producer + +def create_topological(process, inputs, + seed_threshold=5.,# MipT + cluster_threshold=2.# MipT + ): + producer = process.hgcalBackEndLayer1Producer.clone() + producer.ProcessorParameters.C2d_parameters.seeding_threshold_silicon = cms.double(seed_threshold) # MipT + producer.ProcessorParameters.C2d_parameters.seeding_threshold_scintillator = cms.double(seed_threshold) # MipT + producer.ProcessorParameters.C2d_parameters.clustering_threshold_silicon = cms.double(cluster_threshold) # MipT + producer.ProcessorParameters.C2d_parameters.clustering_threshold_scintillator = cms.double(cluster_threshold) # MipT + producer.ProcessorParameters.C2d_parameters.clusterType = cms.string('NNC2d') + producer.InputTriggerCells = cms.InputTag('{}:HGCalConcentratorProcessorSelection'.format(inputs)) + return producer + +def create_constrainedtopological(process, inputs, + distance=6.,# cm + seed_threshold=5.,# MipT + cluster_threshold=2.# MipT + ): + producer = process.hgcalBackEndLayer1Producer.clone() + producer.ProcessorParameters.C2d_parameters.seeding_threshold_silicon = cms.double(seed_threshold) # MipT + producer.ProcessorParameters.C2d_parameters.seeding_threshold_scintillator = cms.double(seed_threshold) # MipT + producer.ProcessorParameters.C2d_parameters.clustering_threshold_silicon = cms.double(cluster_threshold) # MipT + producer.ProcessorParameters.C2d_parameters.clustering_threshold_scintillator = cms.double(cluster_threshold) # MipT + producer.ProcessorParameters.C2d_parameters.dR_cluster = cms.double(distance) # cm + producer.ProcessorParameters.C2d_parameters.clusterType = cms.string('dRNNC2d') + producer.InputTriggerCells = cms.InputTag('{}:HGCalConcentratorProcessorSelection'.format(inputs)) + return producer + +def create_dummy(process, inputs): + producer = process.hgcalBackEndLayer1Producer.clone() + producer.ProcessorParameters.C2d_parameters.clusterType = cms.string('dummyC2d') + producer.InputTriggerCells = cms.InputTag('{}:HGCalConcentratorProcessorSelection'.format(inputs)) + return producer diff --git a/L1Trigger/L1THGCalUtilities/python/clustering3d.py b/L1Trigger/L1THGCalUtilities/python/clustering3d.py new file mode 100644 index 0000000000000..214021d070a96 --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/python/clustering3d.py @@ -0,0 +1,106 @@ +import FWCore.ParameterSet.Config as cms + +from L1Trigger.L1THGCal.customClustering import binSums, dr_layerbylayer + + +def create_distance(process, inputs, + distance=0.01 + ): + producer = process.hgcalBackEndLayer2Producer.clone() + producer.ProcessorParameters.C3d_parameters.dR_multicluster = cms.double(distance) + producer.ProcessorParameters.C3d_parameters.dist_dbscan_multicluster=cms.double(0.) + producer.ProcessorParameters.C3d_parameters.minN_dbscan_multicluster=cms.uint32(0) + producer.ProcessorParameters.C3d_parameters.type_multicluster = cms.string('dRC3d') + producer.InputCluster = cms.InputTag('{}:HGCalBackendLayer1Processor2DClustering'.format(inputs)) + return producer + + +def create_dbscan(process, inputs, + distance=0.005, + min_points=3 + ): + producer = process.hgcalBackEndLayer2Producer.clone() + producer.ProcessorParameters.C3d_parameters.dR_multicluster = cms.double(0.) + producer.ProcessorParameters.C3d_parameters.dist_dbscan_multicluster = cms.double(distance) + producer.ProcessorParameters.C3d_parameters.minN_dbscan_multicluster = cms.uint32(min_points) + producer.ProcessorParameters.C3d_parameters.type_multicluster = cms.string('DBSCANC3d') + producer.InputCluster = cms.InputTag('{}:HGCalBackendLayer1Processor2DClustering'.format(inputs)) + return producer + + +def create_histoMax(process, inputs, + distance = 0.03, + nBins_R = 36, + nBins_Phi = 216, + binSumsHisto = binSums, + seed_threshold = 0, + ): + producer = process.hgcalBackEndLayer2Producer.clone() + producer.ProcessorParameters.C3d_parameters.dR_multicluster = cms.double(distance) + producer.ProcessorParameters.C3d_parameters.nBins_R_histo_multicluster = cms.uint32(nBins_R) + producer.ProcessorParameters.C3d_parameters.nBins_Phi_histo_multicluster = cms.uint32(nBins_Phi) + producer.ProcessorParameters.C3d_parameters.binSumsHisto = binSumsHisto + producer.ProcessorParameters.C3d_parameters.threshold_histo_multicluster = seed_threshold + producer.ProcessorParameters.C3d_parameters.type_multicluster = cms.string('HistoMaxC3d') + producer.InputCluster = cms.InputTag('{}:HGCalBackendLayer1Processor2DClustering'.format(inputs)) + return producer + + +def create_histoMax_variableDr(process, inputs, + distances = dr_layerbylayer, + nBins_R = 36, + nBins_Phi = 216, + binSumsHisto = binSums, + seed_threshold = 0, + ): + producer = create_histoMax(process, inputs, 0, nBins_R, nBins_Phi, binSumsHisto, seed_threshold) + producer.ProcessorParameters.C3d_parameters.dR_multicluster_byLayer = cms.vdouble(distances) + return producer + + +def create_histoInterpolatedMax(process, inputs, + distance = 0.03, + nBins_R = 36, + nBins_Phi = 216, + binSumsHisto = binSums, + ): + producer = create_histoMax( process, inputs, distance, nBins_R, nBins_Phi, binSumsHisto ) + producer.ProcessorParameters.C3d_parameters.type_multicluster = cms.string('HistoInterpolatedMaxC3d') + return producer + +def create_histoInterpolatedMax1stOrder(process, inputs): + producer = create_histoInterpolatedMax( process, inputs ) + producer.ProcessorParameters.C3d_parameters.neighbour_weights=cms.vdouble( 0 , 0.25, 0 , + 0.25 , 0 , 0.25, + 0 , 0.25, 0 + ) + return producer + + + +def create_histoInterpolatedMax2ndOrder(process, inputs): + producer = create_histoInterpolatedMax( process,inputs ) + producer.ProcessorParameters.C3d_parameters.neighbour_weights=cms.vdouble( -0.25, 0.5, -0.25, + 0.5 , 0 , 0.5 , + -0.25, 0.5, -0.25 + ) + return producer + + + +def create_histoThreshold(process, inputs, + threshold = 20., + distance = 0.03, + nBins_R = 36, + nBins_Phi = 216, + binSumsHisto = binSums, + ): + producer = process.hgcalBackEndLayer2Producer.clone() + producer.ProcessorParameters.C3d_parameters.threshold_histo_multicluster = cms.double(threshold) + producer.ProcessorParameters.C3d_parameters.dR_multicluster = cms.double(distance) + producer.ProcessorParameters.C3d_parameters.nBins_R_histo_multicluster = cms.uint32(nBins_R) + producer.ProcessorParameters.C3d_parameters.nBins_Phi_histo_multicluster = cms.uint32(nBins_Phi) + producer.ProcessorParameters.C3d_parameters.binSumsHisto = binSumsHisto + producer.ProcessorParameters.C3d_parameters.type_multicluster = cms.string('HistoThresholdC3d') + producer.InputCluster = cms.InputTag('{}:HGCalBackendLayer1Processor2DClustering'.format(inputs)) + return producer diff --git a/L1Trigger/L1THGCalUtilities/python/concentrator.py b/L1Trigger/L1THGCalUtilities/python/concentrator.py new file mode 100644 index 0000000000000..daad7e7c1dbe6 --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/python/concentrator.py @@ -0,0 +1,47 @@ +import FWCore.ParameterSet.Config as cms +import SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi as digiparam + +def create_supertriggercell(process, inputs, + stcSize = cms.vuint32(4,4,4) + ): + producer = process.hgcalConcentratorProducer.clone() + producer.ProcessorParameters.Method = cms.string('superTriggerCellSelect') + producer.ProcessorParameters.stcSize = stcSize + producer.InputTriggerCells = cms.InputTag('{}:HGCalVFEProcessorSums'.format(inputs)) + producer.InputTriggerSums = cms.InputTag('{}:HGCalVFEProcessorSums'.format(inputs)) + return producer + + +def create_threshold(process, inputs, + threshold=2. # in mipT + ): + adcSaturationBH_MIP = digiparam.hgchebackDigitizer.digiCfg.feCfg.adcSaturation_fC + adcNbitsBH = digiparam.hgchebackDigitizer.digiCfg.feCfg.adcNbits + producer = process.hgcalConcentratorProducer.clone() + producer.ProcessorParameters.Method = cms.string('thresholdSelect') + producer.ProcessorParameters.linLSB = cms.double(100./1024.) + producer.ProcessorParameters.adcsaturationBH = adcSaturationBH_MIP + producer.ProcessorParameters.adcnBitsBH = adcNbitsBH + producer.ProcessorParameters.TCThreshold_fC = cms.double(0.) + producer.ProcessorParameters.TCThresholdBH_MIP = cms.double(0.) + producer.ProcessorParameters.triggercell_threshold_silicon = cms.double(threshold) # MipT + producer.ProcessorParameters.triggercell_threshold_scintillator = cms.double(threshold) # MipT + producer.InputTriggerCells = cms.InputTag('{}:HGCalVFEProcessorSums'.format(inputs)) + producer.InputTriggerSums = cms.InputTag('{}:HGCalVFEProcessorSums'.format(inputs)) + return producer + + +def create_bestchoice(process, inputs, + triggercells=12 + ): + adcSaturationBH_MIP = digiparam.hgchebackDigitizer.digiCfg.feCfg.adcSaturation_fC + adcNbitsBH = digiparam.hgchebackDigitizer.digiCfg.feCfg.adcNbits + producer = process.hgcalConcentratorProducer.clone() + producer.ProcessorParameters.Method = cms.string('bestChoiceSelect') + producer.ProcessorParameters.NData = cms.uint32(triggercells) + producer.ProcessorParameters.linLSB = cms.double(100./1024.) + producer.ProcessorParameters.adcsaturationBH = adcSaturationBH_MIP + producer.ProcessorParameters.adcnBitsBH = adcNbitsBH + producer.InputTriggerCells = cms.InputTag('{}:HGCalVFEProcessorSums'.format(inputs)) + producer.InputTriggerSums = cms.InputTag('{}:HGCalVFEProcessorSums'.format(inputs)) + return producer diff --git a/L1Trigger/L1THGCalUtilities/python/customNtuples.py b/L1Trigger/L1THGCalUtilities/python/customNtuples.py new file mode 100644 index 0000000000000..21de15da5b5f3 --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/python/customNtuples.py @@ -0,0 +1,40 @@ +import FWCore.ParameterSet.Config as cms + +def custom_ntuples_V9(process): + ntuples = process.hgcalTriggerNtuplizer.Ntuples + for ntuple in ntuples: + if ntuple.NtupleName=='HGCalTriggerNtupleHGCDigis' or \ + ntuple.NtupleName=='HGCalTriggerNtupleHGCTriggerCells': + ntuple.bhSimHits = cms.InputTag('g4SimHits:HGCHitsHEback') + return process + + + +def create_ntuple(process, inputs, + ntuple_list=[ + 'event', + 'gen', 'genjet', 'gentau', + 'digis', + 'triggercells', + 'clusters', 'multiclusters' + ] + ): + vpset = [] + for ntuple in ntuple_list: + pset = getattr(process, 'ntuple_'+ntuple).clone() + if ntuple=='triggercells': + pset.TriggerCells = cms.InputTag('{}:HGCalConcentratorProcessorSelection'.format(inputs[0])) + pset.Multiclusters = cms.InputTag('{}:HGCalBackendLayer2Processor3DClustering'.format(inputs[2])) + elif ntuple=='clusters': + pset.Clusters = cms.InputTag('{}:HGCalBackendLayer1Processor2DClustering'.format(inputs[1])) + pset.Multiclusters = cms.InputTag('{}:HGCalBackendLayer2Processor3DClustering'.format(inputs[2])) + elif ntuple=='multiclusters': + pset.Multiclusters = cms.InputTag('{}:HGCalBackendLayer2Processor3DClustering'.format(inputs[2])) + vpset.append(pset) + ntuplizer = process.hgcalTriggerNtuplizer.clone() + ntuplizer.Ntuples = cms.VPSet(vpset) + return ntuplizer + + + + diff --git a/L1Trigger/L1THGCalUtilities/python/hgcalTriggerChains.py b/L1Trigger/L1THGCalUtilities/python/hgcalTriggerChains.py new file mode 100644 index 0000000000000..432ae33aef58d --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/python/hgcalTriggerChains.py @@ -0,0 +1,79 @@ +import FWCore.ParameterSet.Config as cms + +class HGCalTriggerChains: + def __init__(self): + self.vfe = {} + self.concentrator = {} + self.backend1 = {} + self.backend2 = {} + self.ntuple = {} + self.chain = [] + + def register_vfe(self, name, generator): + self.vfe[name] = generator + + def register_concentrator(self, name, generator): + self.concentrator[name] = generator + + def register_backend1(self, name, generator): + self.backend1[name] = generator + + def register_backend2(self, name, generator): + self.backend2[name] = generator + + def register_ntuple(self, name, generator): + self.ntuple[name] = generator + + def register_chain(self, vfe, concentrator, backend1, backend2, ntuple=''): + if not vfe in self.vfe: + raise KeyError('{} not registered as VFE producer'.format(vfe)) + if not concentrator in self.concentrator: + raise KeyError('{} not registered as concentrator producer'.format(concentrator)) + if not backend1 in self.backend1: + raise KeyError('{} not registered as backend1 producer'.format(backend1)) + if not backend2 in self.backend2: + raise KeyError('{} not registered as backend2 producer'.format(backend2)) + if ntuple!='' and not ntuple in self.ntuple: + raise KeyError('{} not registered as ntuplizer'.format(ntuple)) + self.chain.append( (vfe, concentrator, backend1, backend2, ntuple) ) + + + def create_sequences(self, process): + tmp = cms.SequencePlaceholder("tmp") + vfe_sequence = cms.Sequence(tmp) + concentrator_sequence = cms.Sequence(tmp) + backend1_sequence = cms.Sequence(tmp) + backend2_sequence = cms.Sequence(tmp) + ntuple_sequence = cms.Sequence(tmp) + for vfe,concentrator,backend1,backend2,ntuple in self.chain: + concentrator_name = '{0}{1}'.format(vfe, concentrator) + backend1_name = '{0}{1}{2}'.format(vfe, concentrator, backend1) + backend2_name = '{0}{1}{2}{3}'.format(vfe, concentrator, backend1, backend2) + ntuple_name = '{0}{1}{2}{3}{4}'.format(vfe, concentrator, backend1, backend2, ntuple) + ntuple_inputs = [concentrator_name, backend1_name, backend2_name] + if not hasattr(process, vfe): + setattr(process, vfe, self.vfe[vfe](process)) + vfe_sequence *= getattr(process, vfe) + if not hasattr(process, concentrator_name): + setattr(process, concentrator_name, self.concentrator[concentrator](process, vfe)) + concentrator_sequence *= getattr(process, concentrator_name) + if not hasattr(process, backend1_name): + setattr(process, backend1_name, self.backend1[backend1](process, concentrator_name)) + backend1_sequence *= getattr(process, backend1_name) + if not hasattr(process, backend2_name): + setattr(process, backend2_name, self.backend2[backend2](process, backend1_name)) + backend2_sequence *= getattr(process, backend2_name) + if ntuple!='' and not hasattr(process, ntuple_name): + setattr(process, ntuple_name, self.ntuple[ntuple](process, ntuple_inputs)) + ntuple_sequence *= getattr(process, ntuple_name) + vfe_sequence.remove(tmp) + concentrator_sequence.remove(tmp) + backend1_sequence.remove(tmp) + backend2_sequence.remove(tmp) + ntuple_sequence.remove(tmp) + process.globalReplace('hgcalVFE', vfe_sequence) + process.globalReplace('hgcalConcentrator', concentrator_sequence) + process.globalReplace('hgcalBackEndLayer1', backend1_sequence) + process.globalReplace('hgcalBackEndLayer2', backend2_sequence) + process.globalReplace('hgcalTriggerNtuples', ntuple_sequence) + return process diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cff.py b/L1Trigger/L1THGCalUtilities/python/hgcalTriggerNtuples_cff.py similarity index 66% rename from L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cff.py rename to L1Trigger/L1THGCalUtilities/python/hgcalTriggerNtuples_cff.py index 486cc1656b297..4c21948a576f2 100644 --- a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cff.py +++ b/L1Trigger/L1THGCalUtilities/python/hgcalTriggerNtuples_cff.py @@ -1,11 +1,11 @@ import FWCore.ParameterSet.Config as cms -from L1Trigger.L1THGCal.hgcalTriggerNtuples_cfi import * +from L1Trigger.L1THGCalUtilities.hgcalTriggerNtuples_cfi import * hgcalTriggerNtuples = cms.Sequence(hgcalTriggerNtuplizer) from Configuration.Eras.Modifier_phase2_hgcalV9_cff import phase2_hgcalV9 -from L1Trigger.L1THGCal.customNtuples import custom_ntuples_V9 +from L1Trigger.L1THGCalUtilities.customNtuples import custom_ntuples_V9 modifyHgcalTriggerNtuplesWithV9Geometry_ = phase2_hgcalV9.makeProcessModifier(custom_ntuples_V9) diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py b/L1Trigger/L1THGCalUtilities/python/hgcalTriggerNtuples_cfi.py similarity index 95% rename from L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py rename to L1Trigger/L1THGCalUtilities/python/hgcalTriggerNtuples_cfi.py index 5866aafdda998..9b03570e3ed85 100644 --- a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py +++ b/L1Trigger/L1THGCalUtilities/python/hgcalTriggerNtuples_cfi.py @@ -4,7 +4,7 @@ import SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi as digiparam import RecoLocalCalo.HGCalRecProducers.HGCalUncalibRecHit_cfi as recoparam import RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi as recocalibparam -from . import hgcalLayersCalibrationCoefficients_cfi as layercalibparam +import L1Trigger.L1THGCal.hgcalLayersCalibrationCoefficients_cfi as layercalibparam fcPerMip = recoparam.HGCalUncalibRecHit.HGCEEConfig.fCPerMIP @@ -76,7 +76,7 @@ ) from L1Trigger.L1THGCal.egammaIdentification import egamma_identification_drnn_cone -ntuple_multicluster = cms.PSet( +ntuple_multiclusters = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleHGCMulticlusters'), Multiclusters = cms.InputTag('hgcalBackEndLayer2Producer:HGCalBackendLayer2Processor3DClustering'), EGIdentification = egamma_identification_drnn_cone.clone() @@ -87,7 +87,7 @@ TriggerCells = cms.InputTag('hgcalConcentratorProducer:HGCalConcentratorProcessorSelection') ) -ntuple_tower = cms.PSet( +ntuple_towers = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleHGCTowers'), Towers = cms.InputTag('hgcalTowerProducer:HGCalTowerProcessor') ) @@ -101,8 +101,7 @@ ntuple_gentau, ntuple_digis, ntuple_triggercells, - ntuple_clusters, - ntuple_multicluster, - ntuple_tower + ntuple_multiclusters, + ntuple_towers ) ) diff --git a/L1Trigger/L1THGCalUtilities/python/vfe.py b/L1Trigger/L1THGCalUtilities/python/vfe.py new file mode 100644 index 0000000000000..6c3913b3c4eb1 --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/python/vfe.py @@ -0,0 +1,13 @@ +import FWCore.ParameterSet.Config as cms + + +def create_compression(process, + exponent=4, + mantissa=4, + rounding=True + ): + producer = process.hgcalVFEProducer.clone() + producer.ProcessorParameters.exponentBits = cms.uint32(exponent) + producer.ProcessorParameters.mantissaBits = cms.uint32(mantissa) + producer.ProcessorParameters.rounding = cms.bool(rounding) + return producer diff --git a/L1Trigger/L1THGCal/src/HGCalTriggerNtupleBase.cc b/L1Trigger/L1THGCalUtilities/src/HGCalTriggerNtupleBase.cc similarity index 60% rename from L1Trigger/L1THGCal/src/HGCalTriggerNtupleBase.cc rename to L1Trigger/L1THGCalUtilities/src/HGCalTriggerNtupleBase.cc index ec6d7970289d2..c866ee2bdfe25 100644 --- a/L1Trigger/L1THGCal/src/HGCalTriggerNtupleBase.cc +++ b/L1Trigger/L1THGCalUtilities/src/HGCalTriggerNtupleBase.cc @@ -1,4 +1,4 @@ -#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h" EDM_REGISTER_PLUGINFACTORY(HGCalTriggerNtupleFactory, "HGCalTriggerNtupleFactory"); diff --git a/L1Trigger/L1THGCalUtilities/test/testHGCalL1T_RelValV9_cfg.py b/L1Trigger/L1THGCalUtilities/test/testHGCalL1T_RelValV9_cfg.py new file mode 100644 index 0000000000000..086c13a665caa --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/test/testHGCalL1T_RelValV9_cfg.py @@ -0,0 +1,88 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.StandardSequences.Eras import eras + +process = cms.Process('DIGI',eras.Phase2C4) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.Geometry.GeometryExtended2023D35Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2023D35_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.Generator_cff') +process.load('IOMC.EventVertexGenerators.VtxSmearedHLLHC14TeV_cfi') +process.load('GeneratorInterface.Core.genFilterSummary_cff') +process.load('Configuration.StandardSequences.SimIdeal_cff') +process.load('Configuration.StandardSequences.Digi_cff') +process.load('Configuration.StandardSequences.SimL1Emulator_cff') +process.load('Configuration.StandardSequences.DigiToRaw_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(50) +) + +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('/store/relval/CMSSW_10_4_0_pre2/RelValSinglePiFlatPt_0p7to10_pythia8_cfi/GEN-SIM-DIGI-RAW/103X_upgrade2023_realistic_v2_2023D35noPU-v1/10000/CE7F0A8C-F917-F14D-9BC4-CD46CFA8B7FD.root'), + inputCommands=cms.untracked.vstring( + 'keep *', + 'drop l1tEMTFHit2016Extras_simEmtfDigis_CSC_HLT', + 'drop l1tEMTFHit2016Extras_simEmtfDigis_RPC_HLT', + 'drop l1tEMTFHit2016s_simEmtfDigis__HLT', + 'drop l1tEMTFTrack2016Extras_simEmtfDigis__HLT', + 'drop l1tEMTFTrack2016s_simEmtfDigis__HLT', + 'drop FTLClusteredmNewDetSetVector_mtdClusters_FTLBarrel_RECO', + 'drop FTLClusteredmNewDetSetVector_mtdClusters_FTLEndcap_RECO', + 'drop MTDTrackingRecHitedmNewDetSetVector_mtdTrackingRecHits__RECO', + 'drop BTLDetIdBTLSampleFTLDataFrameTsSorted_mix_FTLBarrel_HLT', + 'drop ETLDetIdETLSampleFTLDataFrameTsSorted_mix_FTLEndcap_HLT', + ) + ) + +process.options = cms.untracked.PSet( + +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + version = cms.untracked.string('$Revision: 1.20 $'), + annotation = cms.untracked.string('SingleElectronPt10_cfi nevts:10'), + name = cms.untracked.string('Applications') +) + +# Output definition +process.TFileService = cms.Service( + "TFileService", + fileName = cms.string("ntuple.root") + ) + +# Other statements +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '') + +# load HGCAL TPG simulation +process.load('L1Trigger.L1THGCal.hgcalTriggerPrimitives_cff') +process.hgcl1tpg_step = cms.Path(process.hgcalTriggerPrimitives) +# To test V9Imp2 +#from L1Trigger.L1THGCal.customTriggerGeometry import custom_geometry_V9 +#process = custom_geometry_V9(process, implementation=2) + + +# load ntuplizer +process.load('L1Trigger.L1THGCalUtilities.hgcalTriggerNtuples_cff') +process.ntuple_step = cms.Path(process.hgcalTriggerNtuples) + +# Schedule definition +process.schedule = cms.Schedule(process.hgcl1tpg_step, process.ntuple_step) + +# Add early deletion of temporary data products to reduce peak memory need +from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete +process = customiseEarlyDelete(process) +# End adding early deletion + diff --git a/L1Trigger/L1THGCalUtilities/test/testHGCalL1T_RelVal_cfg.py b/L1Trigger/L1THGCalUtilities/test/testHGCalL1T_RelVal_cfg.py new file mode 100644 index 0000000000000..00a808e58fca8 --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/test/testHGCalL1T_RelVal_cfg.py @@ -0,0 +1,82 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.StandardSequences.Eras import eras +from Configuration.ProcessModifiers.convertHGCalDigisSim_cff import convertHGCalDigisSim + +# For old samples use the digi converter +#process = cms.Process('DIGI',eras.Phase2,convertHGCalDigisSim) +process = cms.Process('DIGI',eras.Phase2) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.Geometry.GeometryExtended2023D17Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2023D17_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.Generator_cff') +process.load('IOMC.EventVertexGenerators.VtxSmearedHLLHC14TeV_cfi') +process.load('GeneratorInterface.Core.genFilterSummary_cff') +process.load('Configuration.StandardSequences.SimIdeal_cff') +process.load('Configuration.StandardSequences.Digi_cff') +process.load('Configuration.StandardSequences.SimL1Emulator_cff') +process.load('Configuration.StandardSequences.DigiToRaw_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(50) +) + +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('/store/relval/CMSSW_10_4_0_pre2/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/103X_upgrade2023_realistic_v2_2023D21noPU-v1/20000/F4344045-AEDE-4240-B7B1-27D2CF96C34E.root'), + inputCommands=cms.untracked.vstring( + 'keep *', + 'drop l1tEMTFHit2016Extras_simEmtfDigis_CSC_HLT', + 'drop l1tEMTFHit2016Extras_simEmtfDigis_RPC_HLT', + 'drop l1tEMTFHit2016s_simEmtfDigis__HLT', + 'drop l1tEMTFTrack2016Extras_simEmtfDigis__HLT', + 'drop l1tEMTFTrack2016s_simEmtfDigis__HLT', + ) + ) + +process.options = cms.untracked.PSet( + +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + version = cms.untracked.string('$Revision: 1.20 $'), + annotation = cms.untracked.string('SingleElectronPt10_cfi nevts:10'), + name = cms.untracked.string('Applications') +) + +# Output definition +process.TFileService = cms.Service( + "TFileService", + fileName = cms.string("ntuple.root") + ) + +# Other statements +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '') + +# load HGCAL TPG simulation +process.load('L1Trigger.L1THGCal.hgcalTriggerPrimitives_cff') +process.hgcl1tpg_step = cms.Path(process.hgcalTriggerPrimitives) + +# load ntuplizer +process.load('L1Trigger.L1THGCalUtilities.hgcalTriggerNtuples_cff') +process.ntuple_step = cms.Path(process.hgcalTriggerNtuples) + +# Schedule definition +process.schedule = cms.Schedule(process.hgcl1tpg_step, process.ntuple_step) + +# Add early deletion of temporary data products to reduce peak memory need +from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete +process = customiseEarlyDelete(process) +# End adding early deletion + diff --git a/L1Trigger/L1THGCalUtilities/test/testHGCalL1T_multialgo_cfg.py b/L1Trigger/L1THGCalUtilities/test/testHGCalL1T_multialgo_cfg.py new file mode 100644 index 0000000000000..42ab1b5672792 --- /dev/null +++ b/L1Trigger/L1THGCalUtilities/test/testHGCalL1T_multialgo_cfg.py @@ -0,0 +1,127 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.StandardSequences.Eras import eras +from Configuration.ProcessModifiers.convertHGCalDigisSim_cff import convertHGCalDigisSim + +# For old samples use the digi converter +#process = cms.Process('DIGI',eras.Phase2,convertHGCalDigisSim) +process = cms.Process('DIGI',eras.Phase2) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.Geometry.GeometryExtended2023D17Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2023D17_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.Generator_cff') +process.load('IOMC.EventVertexGenerators.VtxSmearedHLLHC14TeV_cfi') +process.load('GeneratorInterface.Core.genFilterSummary_cff') +process.load('Configuration.StandardSequences.SimIdeal_cff') +process.load('Configuration.StandardSequences.Digi_cff') +process.load('Configuration.StandardSequences.SimL1Emulator_cff') +process.load('Configuration.StandardSequences.DigiToRaw_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(50) +) + +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('/store/relval/CMSSW_10_4_0_pre2/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/103X_upgrade2023_realistic_v2_2023D21noPU-v1/20000/F4344045-AEDE-4240-B7B1-27D2CF96C34E.root'), + inputCommands=cms.untracked.vstring( + 'keep *', + 'drop l1tEMTFHit2016Extras_simEmtfDigis_CSC_HLT', + 'drop l1tEMTFHit2016Extras_simEmtfDigis_RPC_HLT', + 'drop l1tEMTFHit2016s_simEmtfDigis__HLT', + 'drop l1tEMTFTrack2016Extras_simEmtfDigis__HLT', + 'drop l1tEMTFTrack2016s_simEmtfDigis__HLT', + ) + ) + +process.options = cms.untracked.PSet( + +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + version = cms.untracked.string('$Revision: 1.20 $'), + annotation = cms.untracked.string('SingleElectronPt10_cfi nevts:10'), + name = cms.untracked.string('Applications') +) + +# Output definition +process.TFileService = cms.Service( + "TFileService", + fileName = cms.string("ntuple.root") + ) + +# Other statements +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '') + +# load HGCAL TPG simulation +process.load('L1Trigger.L1THGCal.hgcalTriggerPrimitives_cff') +process.load('L1Trigger.L1THGCalUtilities.hgcalTriggerNtuples_cff') +from L1Trigger.L1THGCalUtilities.hgcalTriggerChains import HGCalTriggerChains +import L1Trigger.L1THGCalUtilities.vfe as vfe +import L1Trigger.L1THGCalUtilities.concentrator as concentrator +import L1Trigger.L1THGCalUtilities.clustering2d as clustering2d +import L1Trigger.L1THGCalUtilities.clustering3d as clustering3d +import L1Trigger.L1THGCalUtilities.customNtuples as ntuple + + +chains = HGCalTriggerChains() +# Register algorithms +## VFE +chains.register_vfe("Floatingpoint8", lambda p : vfe.create_compression(p, 4, 4, True)) +## ECON +chains.register_concentrator("Supertriggercell", concentrator.create_supertriggercell) +chains.register_concentrator("Threshold", concentrator.create_threshold) +## BE1 +chains.register_backend1("Ref2d", clustering2d.create_constrainedtopological) +chains.register_backend1("Dummy", clustering2d.create_dummy) +## BE2 +chains.register_backend2("Ref3d", clustering3d.create_distance) +chains.register_backend2("Histomax", clustering3d.create_histoMax) +chains.register_backend2("Histomaxvardrth0", lambda p,i : clustering3d.create_histoMax_variableDr(p,i,seed_threshold=0.)) +chains.register_backend2("Histomaxvardrth10", lambda p,i : clustering3d.create_histoMax_variableDr(p,i,seed_threshold=10.)) +chains.register_backend2("Histomaxvardrth20", lambda p,i : clustering3d.create_histoMax_variableDr(p,i,seed_threshold=20.)) +# Register ntuples +# Store gen info only in the reference ntuple +ntuple_list_ref = ['event', 'gen', 'multiclusters'] +ntuple_list = ['event', 'multiclusters'] +chains.register_ntuple("Genclustersntuple", lambda p,i : ntuple.create_ntuple(p,i, ntuple_list_ref)) +chains.register_ntuple("Clustersntuple", lambda p,i : ntuple.create_ntuple(p,i, ntuple_list)) + +# Register trigger chains +## Reference chain +chains.register_chain('Floatingpoint8', "Threshold", 'Ref2d', 'Ref3d', 'Genclustersntuple') +concentrator_algos = ['Supertriggercell', 'Threshold'] +backend_algos = ['Histomax', 'Histomaxvardrth0', 'Histomaxvardrth10', 'Histomaxvardrth20'] +## Make cross product fo ECON and BE algos +import itertools +for cc,be in itertools.product(concentrator_algos,backend_algos): + chains.register_chain('Floatingpoint8', cc, 'Dummy', be, 'Clustersntuple') + +process = chains.create_sequences(process) + +# Remove towers from sequence +process.hgcalTriggerPrimitives.remove(process.hgcalTowerMap) +process.hgcalTriggerPrimitives.remove(process.hgcalTower) + +process.hgcl1tpg_step = cms.Path(process.hgcalTriggerPrimitives) +process.ntuple_step = cms.Path(process.hgcalTriggerNtuples) + +# Schedule definition +process.schedule = cms.Schedule(process.hgcl1tpg_step, process.ntuple_step) + +# Add early deletion of temporary data products to reduce peak memory need +from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete +process = customiseEarlyDelete(process) +# End adding early deletion +