Merge pull request #22387 from PFCal-dev/hgc-tpg-integration-180228

HGCAL trigger updates (including new 2D clustering + towers implementation)

DataFormats/L1THGCal/interface/HGCalClusterT.h

@@ -1,22 +1,25 @@
 #ifndef DataFormats_L1Trigger_HGCalClusterT_h
 #define DataFormats_L1Trigger_HGCalClusterT_h
 
+/* CMSSW */
 #include "DataFormats/Common/interface/Ptr.h"
-#include "DataFormats/Common/interface/PtrVector.h"
 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
 #include "DataFormats/L1Trigger/interface/L1Candidate.h"
 #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h"
 #include "DataFormats/L1THGCal/interface/ClusterShapes.h"
+
+/* ROOT */
 #include "Math/Vector3D.h"
 
+#include <unordered_map>
 
 namespace l1t 
 {
   template <class C> class HGCalClusterT : public L1Candidate 
   {
 
     public:
-      typedef typename edm::PtrVector<C>::const_iterator const_iterator;
+      typedef typename std::unordered_map<uint32_t, edm::Ptr<C>>::const_iterator const_iterator;
 
     public:
       HGCalClusterT(){}
@@ -46,61 +49,63 @@ namespace l1t
 
       ~HGCalClusterT() override {};
 
-      const edm::PtrVector<C>& constituents() const {return constituents_;}        
-      const_iterator constituents_begin() const {return constituents_.begin();}
-      const_iterator constituents_end() const {return constituents_.end();}
+      const std::unordered_map<uint32_t, edm::Ptr<C>>& constituents() const { return constituents_; }
+      const_iterator constituents_begin() const { return constituents_.begin(); }
+      const_iterator constituents_end() const { return constituents_.end(); }
       unsigned size() const { return constituents_.size(); }
 
-      void addConstituent( const edm::Ptr<C>& c )
+
+      void addConstituent( const edm::Ptr<C>& c, bool updateCentre=true, float fraction=1. )
       {
-        if( constituents_.empty() )
-        { 
-          detId_ = HGCalDetId(c->detId());
-          seedMipPt_ = c->mipPt();
-        }
 
-        /* update cluster positions */
-        Basic3DVector<float> constituentCentre( c->position() );
-        Basic3DVector<float> clusterCentre( centre_ );
+        double cMipt = c->mipPt()*fraction;
 
-        clusterCentre = clusterCentre*mipPt_ + constituentCentre*c->mipPt();
-        if( mipPt_ + c->mipPt()!=0 ) 
+        if( constituents_.empty() )
         {
-          clusterCentre /= ( mipPt_ + c->mipPt() ) ;
+          detId_ = HGCalDetId( c->detId() );
+          seedMipPt_ = cMipt;
+          /* if the centre will not be dynamically calculated
+             the seed centre is considere as cluster centre */
+          if( !updateCentre )
+          {
+            centre_ = c->position();
+          }
         }
-        centre_ = GlobalPoint( clusterCentre );
+        updateP4AndPosition(c, updateCentre, fraction);
 
-        if( clusterCentre.z()!=0 ) 
-        {
-          centreProj_= GlobalPoint( clusterCentre / clusterCentre.z() );
-        }
-        /* update cluster energies */
-        mipPt_ += c->mipPt();
 
-        int updatedPt = hwPt() + c->hwPt();
-        setHwPt(updatedPt);
+        constituents_.emplace( c->detId(), c );
+        constituentsFraction_.emplace( c->detId(), fraction );
 
-        math::PtEtaPhiMLorentzVector updatedP4 ( p4() );
-        updatedP4 += c->p4(); 
-        setP4( updatedP4 );
+      }
 
-        constituents_.push_back( c );
+      void removeConstituent( const edm::Ptr<C>& c, bool updateCentre=true ){
 
+        /* remove the pointer to c from the std::vector */
+        double fraction=0;
+        const auto& constituent_itr = constituents_.find(c->detId());
+        const auto& fraction_itr = constituentsFraction_.find(c->detId());
+        if(constituent_itr!=constituents_.end())
+        {
+          // remove constituent and get its fraction in the cluster
+          fraction = fraction_itr->second;
+          constituents_.erase(constituent_itr);
+          constituentsFraction_.erase(fraction_itr);
+
+          updateP4AndPosition(c, updateCentre, -fraction);
+        }
       }
-      
-      bool valid() const { return valid_;}
-      void setValid(bool valid) { valid_ = valid;}
-      
+
+      bool valid() const { return valid_; }
+      void setValid(bool valid) { valid_ = valid; }
+
       double mipPt() const { return mipPt_; }
       double seedMipPt() const { return seedMipPt_; }
       uint32_t detId() const { return detId_.rawId(); }
 
 
       /* distance in 'cm' */
-      double distance( const l1t::HGCalTriggerCell &tc ) const 
-      {
-        return ( tc.position() - centre_ ).mag();
-      }
+      double distance( const l1t::HGCalTriggerCell &tc ) const { return ( tc.position() - centre_ ).mag(); }
 
       const GlobalPoint& position() const { return centre_; } 
       const GlobalPoint& centre() const { return centre_; }
@@ -178,11 +183,15 @@ namespace l1t
       bool operator<=(const HGCalClusterT<C>& cl) const { return !(cl>*this); }
       bool operator>=(const HGCalClusterT<C>& cl) const { return !(cl<*this); }
 
+
     private:
-        
+
       bool valid_;
-      HGCalDetId detId_;     
-      edm::PtrVector<C> constituents_;
+      HGCalDetId detId_;
+
+      std::unordered_map<uint32_t, edm::Ptr<C>> constituents_;
+      std::unordered_map<uint32_t, double> constituentsFraction_;
+
       GlobalPoint centre_;
       GlobalPoint centreProj_; // centre projected onto the first HGCal layer
 
@@ -207,6 +216,39 @@ namespace l1t
 
       ClusterShapes shapes_;
 
+
+      void updateP4AndPosition(const edm::Ptr<C>& c, bool updateCentre=true, float fraction=1.)
+      {
+        double cMipt = c->mipPt()*fraction;
+        /* update cluster positions (IF requested) */
+        if( updateCentre ){
+          Basic3DVector<float> constituentCentre( c->position() );
+          Basic3DVector<float> clusterCentre( centre_ );
+
+          clusterCentre = clusterCentre*mipPt_ + constituentCentre*cMipt;
+          if( (mipPt_ + cMipt ) > 0 )
+          {
+            clusterCentre /= ( mipPt_ +  cMipt );
+          }
+          centre_ = GlobalPoint( clusterCentre );
+
+          if( clusterCentre.z()!=0 )
+          {
+            centreProj_= GlobalPoint( clusterCentre / clusterCentre.z() );
+          }
+        }
+
+        /* update cluster energies */
+        mipPt_ += cMipt;
+
+        int updatedPt = hwPt() + (int)(c->hwPt()*fraction);
+        setHwPt( updatedPt );
+
+        math::PtEtaPhiMLorentzVector updatedP4 ( p4() );
+        updatedP4 += (c->p4()*fraction);
+        setP4( updatedP4 );
+      }
+
   };
 
 }

DataFormats/L1THGCal/interface/HGCalTower.h

@@ -43,6 +43,8 @@ namespace l1t {
     int hwEtHad()const;
     int hwEtRatio()const;
 
+    HGCalTower& operator+=(const HGCalTower& tower);
+
   private:
 
     // additional hardware quantities

DataFormats/L1THGCal/interface/HGCalTowerMap.h

@@ -0,0 +1,66 @@
+#ifndef DataFormats_L1TCalorimeter_HGCalTowerMap_h
+#define DataFormats_L1TCalorimeter_HGCalTowerMap_h
+
+#include "DataFormats/L1THGCal/interface/HGCalTower.h"
+#include "DataFormats/L1Trigger/interface/BXVector.h"
+
+#include <unordered_map>
+
+namespace l1t {
+
+  class HGCalTowerMap;
+  typedef BXVector<HGCalTowerMap> HGCalTowerMapBxCollection;
+
+  class HGCalTowerMap {
+
+  public:
+
+    HGCalTowerMap(): nEtaBins_(0), nPhiBins_(0), layer_(0) {}
+
+    HGCalTowerMap( int nEtaBins, int nPhiBins );
+
+    HGCalTowerMap( const std::vector<double>& etaBins, const std::vector<double>& phiBins );
+
+    ~HGCalTowerMap();
+
+    void setLayer( const unsigned layer ) { layer_ = layer; }
+
+
+    int nEtaBins() const { return nEtaBins_; }
+    int nPhiBins() const { return nPhiBins_; }
+    const vector<double>& etaBins() const { return etaBins_; }
+    const vector<double>& phiBins() const { return phiBins_; }
+    const l1t::HGCalTower& tower(int iEta, int iPhi) const { return towerMap_.at(bin_id(iEta,iPhi)); }
+
+    int iEta(const double eta) const;
+    int iPhi(const double phi) const;
+    int layer() const { return layer_;}
+
+    HGCalTowerMap& operator+=(const HGCalTowerMap& map);
+    void addTower(int iEta, int iPhi, const l1t::HGCalTower& tower) { towerMap_[bin_id(iEta,iPhi)] += tower; }
+
+  private:
+
+    static constexpr double kEtaMin_ = 1.479;
+    static constexpr double kEtaMax_ = 3.;
+    static constexpr double kEtaMinLoose_ = 1.401; //BH has some TC below 1.479
+    static constexpr double kEtaMaxLoose_ = 3.085; //FH has some TC above 3.0
+    static constexpr double kPhiMin_ = -M_PI;
+    static constexpr double kPhiMax_ = +M_PI;  
+
+    int nEtaBins_;
+    int nPhiBins_;
+    vector<double> etaBins_;
+    vector<double> phiBins_;
+    std::unordered_map<int,l1t::HGCalTower>  towerMap_;
+    unsigned layer_;
+
+    int bin_id(int iEta,int iPhi) const;
+
+  };
+
+}
+
+#endif
+
+

DataFormats/L1THGCal/src/HGCalTower.cc

@@ -1,4 +1,5 @@
 #include "DataFormats/L1THGCal/interface/HGCalTower.h"
+#include "FWCore/Utilities/interface/EDMException.h"
 
 using namespace l1t;
 
@@ -76,3 +77,25 @@ int HGCalTower::hwEtRatio()const
 {
   return hwEtRatio_;
 }
+
+
+
+
+HGCalTower& HGCalTower::operator+=(const HGCalTower& tower){
+
+  if(this->hwEta()!= tower.hwEta() || this->hwPhi()!= tower.hwPhi()){
+    throw edm::Exception(edm::errors::StdException, "StdException")
+      << "HGCalTower: Trying to add HGCalTowers with different coordinates"<<endl;
+  }
+
+  this->setP4(this->p4() + tower.p4());
+  this->setEtEm(this->etEm() + tower.etEm());
+  this->setEtHad(this->etHad() + tower.etHad());
+
+  this->setHwPt(this->hwPt() + tower.hwPt());
+  this->setHwEtEm(this->hwEtEm() + tower.hwEtEm());
+  this->setHwEtHad(this->hwEtHad() + tower.hwEtHad());
+
+  return *this;
+
+}