HGCAL trigger updates

DataFormats/L1THGCal/interface/HGCalClusterT.h

@@ -95,6 +95,7 @@ namespace l1t
       double seedMipPt() const { return seedMipPt_; }
       uint32_t detId() const { return detId_.rawId(); }
 
+
       /* distance in 'cm' */
       double distance( const l1t::HGCalTriggerCell &tc ) const 
       {
@@ -139,8 +140,38 @@ namespace l1t
       uint32_t subdetId() const {return detId_.subdetId();} 
       uint32_t layer() const {return detId_.layer();}
       int32_t zside() const {return detId_.zside();}
-
 
+
+      //shower shape
+
+      int showerLength() const { return showerLength_; }
+      int coreShowerLength() const { return coreShowerLength_; }
+      int firstLayer() const { return firstLayer_; }
+      int maxLayer() const { return maxLayer_; }
+      float eMax() const { return eMax_; }
+      float sigmaEtaEtaMax() const { return sigmaEtaEtaMax_; }
+      float sigmaPhiPhiMax() const { return sigmaPhiPhiMax_; }
+      float sigmaEtaEtaTot() const { return sigmaEtaEtaTot_; }
+      float sigmaPhiPhiTot() const { return sigmaPhiPhiTot_; }
+      float sigmaZZ() const { return sigmaZZ_; }
+      float sigmaRRTot() const { return sigmaRRTot_; }
+      float sigmaRRMax() const { return sigmaRRMax_; }
+      float sigmaRRMean() const { return sigmaRRMean_; }
+
+      void set_showerLength(int showerLength) { showerLength_ = showerLength;}
+      void set_coreShowerLength(int coreShowerLength) { coreShowerLength_ = coreShowerLength;}
+      void set_firstLayer(int firstLayer) { firstLayer_ = firstLayer;}
+      void set_maxLayer(int maxLayer) { maxLayer_ = maxLayer;}
+      void set_eMax(float eMax) { eMax_ = eMax;}
+      void set_sigmaEtaEtaMax(float sigmaEtaEtaMax) { sigmaEtaEtaMax_ = sigmaEtaEtaMax;}
+      void set_sigmaEtaEtaTot(float sigmaEtaEtaTot) { sigmaEtaEtaTot_ = sigmaEtaEtaTot;}
+      void set_sigmaPhiPhiMax(float sigmaPhiPhiMax) { sigmaPhiPhiMax_ = sigmaPhiPhiMax;}
+      void set_sigmaPhiPhiTot(float sigmaPhiPhiTot) { sigmaPhiPhiTot_ = sigmaPhiPhiTot;}
+      void set_sigmaRRMax(float sigmaRRMax) { sigmaRRMax_ = sigmaRRMax;}
+      void set_sigmaRRTot(float sigmaRRTot) { sigmaRRTot_ = sigmaRRTot;}
+      void set_sigmaRRMean(float sigmaRRMean) { sigmaRRMean_ = sigmaRRMean;}
+      void set_sigmaZZ(float sigmaZZ) { sigmaZZ_ = sigmaZZ;}
+
       /* operators */
       bool operator<(const HGCalClusterT<C>& cl) const {return mipPt() < cl.mipPt();}
       bool operator>(const HGCalClusterT<C>& cl) const  { return  cl<*this;   }
@@ -158,6 +189,22 @@ namespace l1t
       double mipPt_;
       double seedMipPt_;
 
+      //shower shape
+
+      int showerLength_;
+      int coreShowerLength_;
+      int firstLayer_;
+      int maxLayer_;
+      float eMax_;
+      float sigmaEtaEtaMax_;
+      float sigmaPhiPhiMax_;
+      float sigmaRRMax_;
+      float sigmaEtaEtaTot_;
+      float sigmaPhiPhiTot_;
+      float sigmaRRTot_;
+      float sigmaRRMean_;
+      float sigmaZZ_;
+
       ClusterShapes shapes_;
 
   };

DataFormats/L1THGCal/src/classes_def.xml

@@ -1,4 +1,3 @@
-
 <lcgdict>
   <class name="l1t::HGCFETriggerDigi" ClassVersion="10">
     <version ClassVersion="10" checksum="3197268020"/>
@@ -21,16 +20,22 @@
   <class name="edm::Wrapper<l1t::HGCalTowerBxCollection>"/>
 
   <class name="l1t::HGCalClusterT<l1t::HGCalTriggerCell>" />
-  <class name="l1t::HGCalCluster" ClassVersion="10">
-    <version ClassVersion="10" checksum="607544984"/>
+  <class name="l1t::HGCalCluster" ClassVersion="13">
+  <version ClassVersion="13" checksum="3397489079"/>
+  <version ClassVersion="12" checksum="623703096"/>
+  <version ClassVersion="11" checksum="354623225"/>
+  <version ClassVersion="10" checksum="607544984"/>
   </class>
   <class name="std::vector<l1t::HGCalCluster>" />
   <class name="l1t::HGCalClusterBxCollection"/>
   <class name="edm::Wrapper<l1t::HGCalClusterBxCollection>"/>
 
   <class name="l1t::HGCalClusterT<l1t::HGCalCluster>" />
-  <class name="l1t::HGCalMulticluster" ClassVersion="10">
-    <version ClassVersion="10" checksum="1878482802"/>
+  <class name="l1t::HGCalMulticluster" ClassVersion="13">
+  <version ClassVersion="13" checksum="816077951"/>
+  <version ClassVersion="12" checksum="4221677522"/>
+  <version ClassVersion="11" checksum="1508179045"/>
+  <version ClassVersion="10" checksum="1878482802"/>
   </class>
   <class name="std::vector<l1t::HGCalMulticluster>" />
   <class name="l1t::HGCalMulticlusterBxCollection"/>

L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h

@@ -62,6 +62,7 @@ class HGCalTriggerGeometryBase
 
         virtual bool validTriggerCell( const unsigned trigger_cell_id) const = 0;
         virtual bool disconnectedModule(const unsigned module_id) const = 0;
+        virtual unsigned triggerLayer(const unsigned id) const = 0;
 
     protected:
         void setCaloGeometry(const edm::ESHandle<CaloGeometry>& geom) {calo_geometry_=geom;}

L1Trigger/L1THGCal/interface/HGCalTriggerGeometryGenericMapping.h

@@ -139,6 +139,7 @@ class HGCalTriggerGeometryGenericMapping : public HGCalTriggerGeometryBase {
 
   bool validTriggerCell( const unsigned trigger_cell_det_id ) const final;
   bool disconnectedModule(const unsigned module_id) const final;
+  unsigned triggerLayer(const unsigned id) const final;
 
  protected:
   geom_map cells_to_trigger_cells_;

L1Trigger/L1THGCal/interface/be_algorithms/HGCalMulticlusteringImpl.h

@@ -6,6 +6,9 @@
 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 
+#include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h"
+#include "L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h"
+
 class HGCalMulticlusteringImpl{
 
 public:
@@ -16,14 +19,38 @@ class HGCalMulticlusteringImpl{
                       const l1t::HGCalMulticluster & mclu, 
                       double dR ) const;
 
-    void clusterize( const edm::PtrVector<l1t::HGCalCluster> & clustersPtr, 
-                     l1t::HGCalMulticlusterBxCollection & multiclusters);
+    void clusterizeDR( const edm::PtrVector<l1t::HGCalCluster> & clustersPtr, 
+                     l1t::HGCalMulticlusterBxCollection & multiclusters,
+                     const HGCalTriggerGeometryBase & triggerGeometry
+                     );
+
+    void clusterizeDBSCAN( const edm::PtrVector<l1t::HGCalCluster> & clustersPtr, 
+                     l1t::HGCalMulticlusterBxCollection & multiclusters,
+                     const HGCalTriggerGeometryBase & triggerGeometry
+                     );
 
 private:
+
+    void findNeighbor( const std::vector<std::pair<unsigned int,double>>&  rankedList,
+                       unsigned int searchInd,
+                       const edm::PtrVector<l1t::HGCalCluster> & clustersPtr, 
+                       std::vector<unsigned int>& neigbors);
 
     double dr_;
     double ptC3dThreshold_;
     double calibSF_;
+    std::string multiclusterAlgoType_;
+    double distDbscan_ = 0.005;
+    unsigned minNDbscan_ = 3;
+    std::vector<double> layerWeights_;
+    bool applyLayerWeights_;
+
+    HGCalShowerShape shape_;
+
+    static const int kLayersEE_=28;
+    static const int kLayersFH_=12;
+    static const int kLayersBH_=12;
+
 };
 
 #endif

L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h

@@ -0,0 +1,58 @@
+#ifndef __L1Trigger_L1THGCal_HGCalShowerShape_h__
+#define __L1Trigger_L1THGCal_HGCalShowerShape_h__
+#include <vector>
+#include <cmath>
+#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
+#include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h"
+
+class HGCalShowerShape{
+
+    public:
+    typedef math::XYZTLorentzVector LorentzVector;
+
+    HGCalShowerShape(){}
+
+    ~HGCalShowerShape(){}
+
+    int firstLayer(const l1t::HGCalMulticluster& c3d) const;
+    int lastLayer(const l1t::HGCalMulticluster& c3d) const;
+    int maxLayer(const l1t::HGCalMulticluster& c3d) const;
+    int showerLength(const l1t::HGCalMulticluster& c3d) const {return lastLayer(c3d)-firstLayer(c3d)+1; }//in number of layers
+    // Maximum number of consecutive layers in the cluster
+    int coreShowerLength(const l1t::HGCalMulticluster& c3d, const HGCalTriggerGeometryBase& triggerGeometry) const;
+
+    float eMax(const l1t::HGCalMulticluster& c3d) const;  
+
+    float sigmaZZ(const l1t::HGCalMulticluster& c3d) const;
+
+    float sigmaEtaEtaTot(const l1t::HGCalMulticluster& c3d) const;
+    float sigmaEtaEtaTot(const l1t::HGCalCluster& c2d) const;       
+    float sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const;   
+
+    float sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const;
+    float sigmaPhiPhiTot(const l1t::HGCalCluster& c2d) const;       
+    float sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const;   
+
+    float sigmaRRTot(const l1t::HGCalMulticluster& c3d) const;
+    float sigmaRRTot(const l1t::HGCalCluster& c2d) const;       
+    float sigmaRRMax(const l1t::HGCalMulticluster& c3d) const;  
+    float sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius=5.) const;
+
+    private: 
+
+    float meanX(const std::vector<pair<float,float> >& energy_X_tc) const;
+    float sigmaXX(const std::vector<pair<float,float> >& energy_X_tc, const float X_cluster) const;
+    float sigmaPhiPhi(const std::vector<pair<float,float> >& energy_phi_tc, const float phi_cluster) const;
+
+    static const int kLayersEE_=28;
+    static const int kLayersFH_=12;
+    static const int kLayersBH_=12;
+    int HGC_layer(const uint32_t subdet, const uint32_t layer) const;
+
+
+};
+
+
+#endif
+

L1Trigger/L1THGCal/plugins/BuildFile.xml

@@ -32,5 +32,6 @@
   <use   name="SimDataFormats/CaloTest"/>
   <use   name="DataFormats/JetReco"/>
   <use   name="Geometry/HcalTowerAlgo"/>
+  <use   name="SimDataFormats/PileupSummaryInfo"/>
   <flags   EDM_PLUGIN="1"/>
 </library>

L1Trigger/L1THGCal/plugins/be_algorithms/HGCClusterAlgo.cc

@@ -28,6 +28,10 @@ class HGCClusterAlgo : public Algorithm<FECODEC>
             dRC2d,
             NNC2d
         };
+        enum MulticlusterType{
+            dRC3d,
+            DBSCANC3d
+        };
 
     public:
 
@@ -38,10 +42,10 @@ class HGCClusterAlgo : public Algorithm<FECODEC>
         multicluster_product_( new l1t::HGCalMulticlusterBxCollection ),
         calibration_( conf.getParameterSet("calib_parameters") ),
         clustering_( conf.getParameterSet("C2d_parameters") ),
-        multiclustering_( conf.getParameterSet("C3d_parameters" ) )
+        multiclustering_( conf.getParameterSet("C3d_parameters" ) ),
+        triggercell_threshold_silicon_( conf.getParameter<double>("triggercell_threshold_silicon") ),
+        triggercell_threshold_scintillator_( conf.getParameter<double>("triggercell_threshold_scintillator") )
         {
-            clustering_threshold_silicon_ = conf.getParameterSet("C2d_parameters").getParameter<double>("clustering_threshold_silicon");
-            clustering_threshold_scintillator_ = conf.getParameterSet("C2d_parameters").getParameter<double>("clustering_threshold_scintillator");
             std::string type(conf.getParameterSet("C2d_parameters").getParameter<std::string>("clusterType"));
             if(type=="dRC2d"){
                 clusteringAlgorithmType_ = dRC2d;
@@ -52,6 +56,16 @@ class HGCClusterAlgo : public Algorithm<FECODEC>
                     << "'. Using nearest neighbor NNC2d instead.\n";
                 clusteringAlgorithmType_ = NNC2d;
             }
+            std::string typeMulticluster(conf.getParameterSet("C3d_parameters").getParameter<std::string>("type_multicluster"));
+            if(typeMulticluster=="dRC3d"){
+                multiclusteringAlgoType_ = dRC3d;
+            }else if(typeMulticluster=="DBSCANC3d"){
+                multiclusteringAlgoType_ = DBSCANC3d;
+            }else {
+                edm::LogWarning("ParameterError") << "Unknown Multiclustering type '" << typeMulticluster
+                    << "'. Using Cone Algorithm instead.\n";
+                multiclusteringAlgoType_ = dRC3d;
+            }
 
         }
 
@@ -95,8 +109,9 @@ class HGCClusterAlgo : public Algorithm<FECODEC>
 
         /* algorithm type */
         ClusterType clusteringAlgorithmType_;
-        double clustering_threshold_silicon_;
-        double clustering_threshold_scintillator_;
+        double triggercell_threshold_silicon_;
+        double triggercell_threshold_scintillator_;
+        MulticlusterType multiclusteringAlgoType_;
 };
 
 
@@ -122,8 +137,8 @@ void HGCClusterAlgo<FECODEC,DATA>::run(const l1t::HGCFETriggerDigiCollection & c
             {
                 l1t::HGCalTriggerCell calibratedtriggercell( triggercell );
                 calibration_.calibrateInGeV( calibratedtriggercell); 
-                double clustering_threshold = (triggercell.subdetId()==HGCHEB ? clustering_threshold_scintillator_ : clustering_threshold_silicon_);
-                if(calibratedtriggercell.mipPt()<clustering_threshold) continue;
+                double triggercell_threshold = (triggercell.subdetId()==HGCHEB ? triggercell_threshold_scintillator_ : triggercell_threshold_silicon_);
+                if(calibratedtriggercell.mipPt()<triggercell_threshold) continue;
                 trgcell_product_->push_back( 0, calibratedtriggercell );
             }           
 
@@ -169,12 +184,24 @@ void HGCClusterAlgo<FECODEC,DATA>::run(const l1t::HGCFETriggerDigiCollection & c
         clustersPtrs.push_back(ptr);
     }
 
-    /* call to multiclustering */
-    multiclustering_.clusterize( clustersPtrs, *multicluster_product_ );
+    /* call to multiclustering and compute shower shape*/
+    switch(multiclusteringAlgoType_){
+        case dRC3d : 
+            multiclustering_.clusterizeDR( clustersPtrs, *multicluster_product_, *triggerGeometry_);
+            break;
+        case DBSCANC3d:
+            multiclustering_.clusterizeDBSCAN( clustersPtrs, *multicluster_product_, *triggerGeometry_);
+            break;
+        default:
+            // Should not happen, clustering type checked in constructor
+            break;
+    }
 
     /* retrieve the orphan handle to the multiclusters collection and put the collection in the event */
     multiclustersHandle = evt.put( std::move( multicluster_product_ ), "cluster3D");
 
+
+
 }
 
 typedef HGCClusterAlgo<HGCalTriggerCellBestChoiceCodec, HGCalTriggerCellBestChoiceCodec::data_type> HGCClusterAlgoBestChoice;

L1Trigger/L1THGCal/plugins/geometries/HGCalTriggerGeometryHexImp2.cc

@@ -36,6 +36,7 @@ class HGCalTriggerGeometryHexImp2 : public HGCalTriggerGeometryBase
 
         bool validTriggerCell( const unsigned ) const final;
         bool disconnectedModule(const unsigned) const final;
+        unsigned triggerLayer(const unsigned) const final;
 
     private:
         edm::FileInPath l1tCellsMapping_;
@@ -773,6 +774,14 @@ disconnectedModule(const unsigned module_id) const
     return false;
 }
 
+
+unsigned 
+HGCalTriggerGeometryHexImp2::
+triggerLayer(const unsigned id) const
+{
+    return HGCalDetId(id).layer();
+}
+
 bool 
 HGCalTriggerGeometryHexImp2::
 validTriggerCellFromCells(const unsigned trigger_cell_id) const