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

DataFormats/L1THGCal/interface/HGCalClusterT.h

@@ -1,22 +1,23 @@
 #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"
-#include "Math/Vector3D.h"
 
+/* ROOT */
+#include "Math/Vector3D.h"
 
 namespace l1t 
 {
   template <class C> class HGCalClusterT : public L1Candidate 
   {
 
     public:
-      typedef typename edm::PtrVector<C>::const_iterator const_iterator;
+      typedef typename std::vector<edm::Ptr<C>>::const_iterator const_iterator;
 
     public:
       HGCalClusterT(){}
@@ -46,61 +47,126 @@ 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::vector<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 );
 
-        if( clusterCentre.z()!=0 ) 
-        {
-          centreProj_= GlobalPoint( clusterCentre / clusterCentre.z() );
+        /* 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_ += c->mipPt();
+        mipPt_ += cMipt;
 
-        int updatedPt = hwPt() + c->hwPt();
-        setHwPt(updatedPt);
+        int updatedPt = hwPt() + (int)(c->hwPt()*fraction);
+        setHwPt( updatedPt );
 
         math::PtEtaPhiMLorentzVector updatedP4 ( p4() );
-        updatedP4 += c->p4(); 
+        updatedP4 += (c->p4()*fraction);
         setP4( updatedP4 );
 
         constituents_.push_back( c );
+        constituentsFraction_.push_back( fraction );
 
       }
-
-      bool valid() const { return valid_;}
-      void setValid(bool valid) { valid_ = valid;}
-
+
+      void removeConstituent( const edm::Ptr<C>& c, bool updateCentre=true ){
+
+        /* remove the pointer to c from the std::vector */
+        double fraction=0;
+        bool constituentRemoved=false;
+        for( unsigned i=0; i<constituents_.size(); i++ )
+        {
+          if( constituents_[i] == c )
+          {
+            // remove constituent and get its fraction in the cluster
+            constituents_.erase( constituents_.begin()+i );
+            fraction = constituentsFraction_.at(i);
+            constituentsFraction_.erase( constituentsFraction_.begin()+i );
+            constituentRemoved=true;
+            break;
+          }
+        }
+
+        /* if a constituent has been removed update cluster info */
+        if( constituentRemoved ) {
+
+          /* update cluster positions (IF requested) */
+          double cMipt = c->mipPt()*fraction;
+          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() - ( c->hwPt()*fraction );
+          setHwPt( updatedPt );
+
+          math::PtEtaPhiMLorentzVector updatedP4 ( p4() );
+          updatedP4 -= ( c->p4()*fraction );
+          setP4( updatedP4 );
+
+        }
+
+      }
+
+      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 +244,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::vector<edm::Ptr<C>> constituents_;  /* ???? possibly change this in something like     */
+      std::vector<double> constituentsFraction_;  /*    vector<pair<edm::Ptr<C>,float>>    ????  */
+
       GlobalPoint centre_;
       GlobalPoint centreProj_; // centre projected onto the first HGCal layer
 

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;
+
+}