[L1T] PR containing simulation and emulation related to jet mass reconstruction (Phase-2)

DataFormats/L1TParticleFlow/interface/PFJet.h

@@ -15,7 +15,7 @@ namespace l1t {
     typedef std::vector<edm::Ptr<l1t::PFCandidate>> Constituents;
 
     PFJet() {}
-    PFJet(float pt, float eta, float phi, float mass = 0, int hwpt = 0, int hweta = 0, int hwphi = 0)
+    PFJet(float pt, float eta, float phi, float mass, int hwpt = 0, int hweta = 0, int hwphi = 0)
         : L1Candidate(PolarLorentzVector(pt, eta, phi, mass), hwpt, hweta, hwphi, /*hwQuality=*/0), rawPt_(pt) {}
 
     PFJet(const LorentzVector& p4, int hwpt = 0, int hweta = 0, int hwphi = 0)
@@ -67,7 +67,7 @@ namespace l1t {
 
     // Get and set the encodedJet_ bits. The Jet is encoded in 128 bits as a 2-element array of uint64_t
     // We store encodings both for Correlator internal usage and for Global Trigger
-    enum class HWEncoding { CT, GT };
+    enum class HWEncoding { CT, GT, GTWide };
     typedef std::array<uint64_t, 2> PackedJet;
     const PackedJet& encodedJet(const HWEncoding encoding = HWEncoding::GT) const {
       return encodedJet_[static_cast<int>(encoding)];
@@ -78,6 +78,7 @@ namespace l1t {
 
     // Accessors to HW objects with ap_* types from encoded words
     const PackedJet& getHWJetGT() const { return encodedJet(HWEncoding::GT); }
+    const PackedJet& getHWJetGTWide() const { return encodedJet(HWEncoding::GTWide); }
     const PackedJet& getHWJetCT() const { return encodedJet(HWEncoding::CT); }
 
   private:
@@ -86,7 +87,7 @@ namespace l1t {
     std::vector<l1ct::JetTagClass> tagClasses_;
     std::vector<float> tagScores_;
     float ptCorrection_;
-    std::array<PackedJet, 2> encodedJet_ = {{{{0, 0}}, {{0, 0}}}};
+    std::array<PackedJet, 3> encodedJet_ = {{{{0, 0}}, {{0, 0}}, {{0, 0}}}};
   };
 
   typedef std::vector<l1t::PFJet> PFJetCollection;

DataFormats/L1TParticleFlow/interface/datatypes.h

@@ -8,6 +8,7 @@
 namespace l1ct {
 
   typedef ap_ufixed<14, 12, AP_TRN, AP_SAT> pt_t;
+  typedef ap_ufixed<17, 15, AP_TRN, AP_SAT> mass2_t;  // can store up to 256 GeV with 0.5 GeV precision
   typedef ap_ufixed<10, 8, AP_TRN, AP_SAT> pt10_t;
   typedef ap_fixed<16, 14, AP_TRN, AP_SAT> dpt_t;
   typedef ap_ufixed<28, 24, AP_TRN, AP_SAT> pt2_t;
@@ -177,6 +178,7 @@ namespace l1ct {
     inline float floatMeanZ(meanz_t meanz) { return meanz + MEANZ_OFFSET; };
     inline float floatHoe(hoe_t hoe) { return hoe.to_float(); };
     inline float floatIDScore(id_score_t score) { return score.to_float(); };
+    inline float floatMass(mass2_t mass) { return mass.to_float(); }
 
     inline pt_t makePt(int pt) { return ap_ufixed<16, 14>(pt) >> 2; }
     inline dpt_t makeDPt(int dpt) { return ap_fixed<18, 16>(dpt) >> 2; }

DataFormats/L1TParticleFlow/interface/gt_datatypes.h

@@ -27,6 +27,8 @@ namespace l1gt {
   // While bitwise identical to the l1ct::z0_t value, we store z0 in mm to profit of ap_fixed goodies
   typedef ap_fixed<10, 9, AP_RND_CONV, AP_SAT> z0_t;   // NOTE: mm instead of cm!!!
   typedef ap_ufixed<8, 1, AP_RND, AP_SAT> id_proba_t;  // for IDs bounded in range [0-1]
+  typedef ap_ufixed<10, 1, AP_RND, AP_SAT> n_prong_score_t;
+  typedef ap_ufixed<17, 15, AP_RND_CONV, AP_SAT> mass2_t;
   typedef ap_uint<1> valid_t;
 
   // E/gamma fields
@@ -48,6 +50,7 @@ namespace l1gt {
     inline float floatEta(eta_t eta) { return eta.to_float() * ETAPHI_LSB; }
     inline float floatPhi(phi_t phi) { return phi.to_float() * ETAPHI_LSB; }
     inline float floatZ0(z0_t z0) { return z0.to_float() * Z0_UNITS; }
+    inline float floatMassSq(mass2_t massSq) { return massSq.to_float(); }
   }  // namespace Scales
 
   struct ThreeVector {
@@ -151,8 +154,90 @@ namespace l1gt {
       return unpack_ap(bits);
     }
 
+    inline static Jet unpack(const std::array<long long unsigned int, 2> &src) {
+      // unpack from two 64b ints
+      ap_uint<BITWIDTH> bits;
+      bits(63, 0) = src[0];
+      bits(127, 64) = src[1];
+      return unpack_ap(bits);
+    }
+
   };  // struct Jet
 
+  struct WideJet {
+    valid_t valid;
+    ThreeVector v3;
+    z0_t z0;
+    n_prong_score_t hwNProngScore;
+    mass2_t hwMassSq;
+
+    inline bool operator==(const WideJet &other) const {
+      return valid == other.valid && z0 == other.z0 && hwNProngScore == other.hwNProngScore &&
+             hwMassSq == other.hwMassSq && v3 == other.v3;
+    }
+
+    static const int BITWIDTH = 128;
+    inline ap_uint<BITWIDTH> pack_ap() const {
+      ap_uint<BITWIDTH> ret = 0;
+      unsigned int start = 0;
+      pack_into_bits(ret, start, valid);
+      pack_into_bits(ret, start, v3.pack());
+      pack_into_bits(ret, start, z0);
+      pack_into_bits(ret, start, hwNProngScore);
+      start = 64;  // Start second word
+      pack_into_bits(ret, start, hwMassSq);
+      return ret;
+    }
+
+    inline std::array<uint64_t, 2> pack() const {
+      std::array<uint64_t, 2> packed;
+      ap_uint<BITWIDTH> bits = this->pack_ap();
+      packed[0] = bits(63, 0);
+      packed[1] = bits(127, 64);
+      return packed;
+    }
+
+    inline static WideJet unpack_ap(const ap_uint<BITWIDTH> &src) {
+      WideJet ret;
+      ret.initFromBits(src);
+      return ret;
+    }
+
+    inline void initFromBits(const ap_uint<BITWIDTH> &src) {
+      unsigned int start = 0;
+      unpack_from_bits(src, start, valid);
+      unpack_from_bits(src, start, v3.pt);
+      unpack_from_bits(src, start, v3.phi);
+      unpack_from_bits(src, start, v3.eta);
+      unpack_from_bits(src, start, z0);
+      unpack_from_bits(src, start, hwNProngScore);
+      start = 64;  // Start second word
+      unpack_from_bits(src, start, hwMassSq);
+    }
+
+    inline static WideJet unpack(const std::array<uint64_t, 2> &src) {
+      ap_uint<BITWIDTH> bits;
+      bits(63, 0) = src[0];
+      bits(127, 64) = src[1];
+      return unpack_ap(bits);
+    }
+
+    inline static WideJet unpack(long long unsigned int &src) {
+      // unpack from single 64b int
+      ap_uint<BITWIDTH> bits = src;
+      return unpack_ap(bits);
+    }
+
+    inline static WideJet unpack(const std::array<long long unsigned int, 2> &src) {
+      // unpack from two 64b ints
+      ap_uint<BITWIDTH> bits;
+      bits(63, 0) = src[0];
+      bits(127, 64) = src[1];
+      return unpack_ap(bits);
+    }
+
+  };  // struct WideJet
+
   struct Sum {
     valid_t valid;
     pt_t vector_pt;
@@ -394,6 +479,8 @@ namespace l1ct {
     return x * Scales::ETAPHI_CTtoGT_SCALE;
   }
 
+  inline l1gt::mass2_t CTtoGT_massSq(mass2_t x) { return (l1gt::mass2_t)x; }
+
 }  // namespace l1ct
 
 #endif

DataFormats/L1TParticleFlow/interface/jets.h

@@ -58,9 +58,11 @@ namespace l1ct {
 
     static const unsigned NTagFields = 8;
     jet_tag_score_t hwTagScores[NTagFields];
+    mass2_t hwMassSq;
 
     inline bool operator==(const Jet &other) const {
-      bool eq = hwPt == other.hwPt && hwEta == other.hwEta && hwPhi == other.hwPhi && hwZ0 == other.hwZ0;
+      bool eq = hwPt == other.hwPt && hwEta == other.hwEta && hwPhi == other.hwPhi && hwMassSq == other.hwMassSq &&
+                hwZ0 == other.hwZ0;
       for (unsigned i = 0; i < NTagFields; i++) {
         eq = eq && hwTagScores[i] == other.hwTagScores[i];
       }
@@ -74,6 +76,7 @@ namespace l1ct {
       hwPt = 0;
       hwEta = 0;
       hwPhi = 0;
+      hwMassSq = 0;
       hwZ0 = 0;
       for (unsigned i = 0; i < NTagFields; i++) {
         hwTagScores[i] = 0;
@@ -87,6 +90,7 @@ namespace l1ct {
     float floatEta() const { return Scales::floatEta(hwEta); }
     float floatPhi() const { return Scales::floatPhi(hwPhi); }
     float floatZ0() const { return Scales::floatZ0(hwZ0); }
+    float floatMass() const { return Scales::floatMass(hwMassSq); }
     std::vector<float> floatIDScores() const {
       std::vector<float> scores(NTagFields);
       for (unsigned i = 0; i < NTagFields; i++) {
@@ -96,13 +100,15 @@ namespace l1ct {
     }
 
     static const int BITWIDTH =
-        pt_t::width + glbeta_t::width + glbphi_t::width + z0_t::width + NTagFields * id_score_t::width;
+        pt_t::width + glbeta_t::width + glbphi_t::width + mass2_t::width + z0_t::width + NTagFields * id_score_t::width;
+
     inline ap_uint<BITWIDTH> pack_ap() const {
       ap_uint<BITWIDTH> ret;
       unsigned int start = 0;
       pack_into_bits(ret, start, hwPt);
       pack_into_bits(ret, start, hwEta);
       pack_into_bits(ret, start, hwPhi);
+      pack_into_bits(ret, start, hwMassSq);
       pack_into_bits(ret, start, hwZ0);
       for (unsigned i = 0; i < NTagFields; i++) {
         pack_into_bits(ret, start, hwTagScores[i]);
@@ -130,6 +136,7 @@ namespace l1ct {
       unpack_from_bits(src, start, hwEta);
       unpack_from_bits(src, start, hwPhi);
       unpack_from_bits(src, start, hwZ0);
+      unpack_from_bits(src, start, hwMassSq);
       for (unsigned i = 0; i < NTagFields; i++) {
         unpack_from_bits(src, start, hwTagScores[i]);
       }
@@ -149,6 +156,14 @@ namespace l1ct {
       return unpack_ap(bits);
     }
 
+    inline static Jet unpack(const std::array<long long unsigned int, 2> &src) {
+      // unpack from two 64b ints
+      ap_uint<BITWIDTH> bits;
+      bits(63, 0) = src[0];
+      // bits(127, 64) = src[1];
+      return unpack_ap(bits);
+    }
+
     l1gt::Jet toGT() const {
       l1gt::Jet j;
       j.valid = hwPt != 0;
@@ -161,6 +176,19 @@ namespace l1ct {
       }
       return j;
     }
+
+    l1gt::WideJet toGTWide() const {
+      l1gt::WideJet j;
+      j.valid = hwPt != 0;
+      j.v3.pt = CTtoGT_pt(hwPt);
+      j.v3.phi = CTtoGT_phi(hwPhi);
+      j.v3.eta = CTtoGT_eta(hwEta);
+      j.z0(l1ct::z0_t::width - 1, 0) = hwZ0(l1ct::z0_t::width - 1, 0);
+      j.hwNProngScore = 0;
+      j.hwMassSq = CTtoGT_massSq(hwMassSq);
+
+      return j;
+    }
   };
 
   inline void clear(Jet &c) { c.clear(); }

DataFormats/L1TParticleFlow/src/classes_def.xml

@@ -35,8 +35,9 @@
   <class name="l1t::RegionalOutput<l1t::PFCandidateCollection>" />
   <class name="edm::Wrapper<l1t::RegionalOutput<l1t::PFCandidateCollection>>" />
 
-  <class name="l1t::PFJet"  ClassVersion="7">
-	<version ClassVersion="7" checksum="2216717791"/>
+  <class name="l1t::PFJet"  ClassVersion="8">
+        <version ClassVersion="8" checksum="2216717794"/>
+        <version ClassVersion="7" checksum="2216717791"/>
         <version ClassVersion="6" checksum="2599349078"/>
         <version ClassVersion="5" checksum="2270932343"/>
         <version ClassVersion="4" checksum="1424452548"/>

L1Trigger/Configuration/python/SimL1Emulator_cff.py

@@ -253,6 +253,7 @@
 # PF MET
 # ########################################################################
 from L1Trigger.Phase2L1ParticleFlow.l1pfJetMet_cff import *
+_phase2_siml1emulator.add(L1TPFJetsTask)
 _phase2_siml1emulator.add(L1TPFJetsEmulationTask)
 
 from L1Trigger.Phase2L1ParticleFlow.l1tMETPFProducer_cfi import *

L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1SeedConePFJetEmulator.h

@@ -1,6 +1,8 @@
 #ifndef L1Trigger_Phase2L1ParticleFlow_L1SeedConePFJetEmulator_h
 #define L1Trigger_Phase2L1ParticleFlow_L1SeedConePFJetEmulator_h
 
+#define NCONSTITSFW 32  // DEFINE THE MAXIMUM NUMBER OF CONSTITUENTS USED TO CALCULATE THE JET MASS
+
 #include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h"
 #include "DataFormats/L1TParticleFlow/interface/jets.h"
 
@@ -21,6 +23,18 @@ class L1SCJetEmu {
   typedef ap_int<13> detaphi_t;          // Type for deta & dphi
   typedef ap_fixed<18, 23> detaphi2_t;   // Type for deta^2 & dphi^2
   typedef ap_fixed<22, 22> pt_etaphi_t;  // Type for product of pt with deta & dphi
+
+  typedef ap_ufixed<13, 1, AP_RND, AP_SAT>
+      eventrig_t;  // stores values between 0 and 2, - 0 bit for sign, - 1 bit for integer, leaves 12 for frac
+  typedef ap_fixed<13, 1, AP_RND, AP_SAT>
+      oddtrig_t;  // stores values between -1 and 1, 13 - 1 bit for sign, - 0 bits for integer, leaves 12 for frac
+
+  // typedef l1ct::mass_t mass_t;  // stores values up to ~1 TeV, 18 bits - 0 for sign, - 10 for integer, 14 total bits improves performance
+  typedef l1ct::mass2_t mass2_t;
+
+  typedef ap_ufixed<20, 12, AP_TRN, AP_SAT> ppt_t;  // stores values between -1 and 1
+  typedef ap_fixed<22, 14, AP_TRN, AP_SAT> npt_t;   // stores values between -1 and 1        JUST REDUCED BY 2
+
   typedef l1ct::PuppiObjEmu Particle;
 
   class Jet : public l1ct::Jet {
@@ -43,6 +57,7 @@ class L1SCJetEmu {
   typedef ap_ufixed<18, -2> inv_pt_t;
   static constexpr int N_table_inv_pt = 1024;
   inv_pt_t inv_pt_table_[N_table_inv_pt];
+  static constexpr int hwEtaPhi_steps = 185;  // corresponds to eta/phi range of 0 to 0.8
 
   static constexpr int ceillog2(int x) { return (x <= 2) ? 1 : 1 + ceillog2((x + 1) / 2); }
 
@@ -142,10 +157,20 @@ class L1SCJetEmu {
     return out;
   }
 
+  template <typename lut_T, int N>
+  static std::array<lut_T, N> init_trig_lut(lut_T (*func)(float)) {
+    std::array<lut_T, N> lut;
+    for (int hwEtaPhi = 0; hwEtaPhi < N; hwEtaPhi++) {
+      lut[hwEtaPhi] = func(hwEtaPhi * l1ct::Scales::ETAPHI_LSB);
+    }
+    return lut;
+  }
+
   static detaphi_t deltaPhi(Particle a, Particle b);
   bool inCone(Particle seed, Particle part) const;
-  Jet makeJet_HW(const std::vector<Particle>& parts) const;
-
+  std::vector<Particle> sortConstituents(const std::vector<Particle>& parts, const Particle seed) const;
+  mass2_t jetMass_HW(const std::vector<Particle>& parts) const;
+  Jet makeJet_HW(const std::vector<Particle>& parts, const Particle seed) const;
 };  // class L1SCJetEmu
 
 #endif