Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:21:26

0001 #ifndef L1Trigger_Phase2L1ParticleFlow_HTMHT_h
0002 #define L1Trigger_Phase2L1ParticleFlow_HTMHT_h
0003 
0004 #include "DataFormats/L1TParticleFlow/interface/jets.h"
0005 #include "DataFormats/L1TParticleFlow/interface/sums.h"
0006 #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
0007 #include "L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1SeedConePFJetEmulator.h"
0008 
0009 #ifndef CMSSW_GIT_HASH
0010 #include "hls_math.h"
0011 #endif
0012 
0013 #include <vector>
0014 #include <numeric>
0015 #include <algorithm>
0016 #include "ap_int.h"
0017 #include "ap_fixed.h"
0018 
0019 namespace P2L1HTMHTEmu {
0020   typedef l1ct::pt_t pt_t;          // Type for pt/ht 1 unit = 0.25 GeV; max = 16 TeV
0021   typedef l1ct::glbeta_t etaphi_t;  // Type for eta & phi
0022 
0023   typedef ap_fixed<12, 3> radians_t;
0024   typedef ap_fixed<9, 2> cossin_t;
0025   typedef ap_fixed<16, 13> pxy_t;
0026 
0027   static constexpr int N_TABLE = 2048;
0028 
0029   // Class for intermediate variables
0030   class PtPxPy {
0031   public:
0032     pt_t pt = 0.;
0033     pxy_t px = 0.;
0034     pxy_t py = 0.;
0035 
0036     PtPxPy operator+(const PtPxPy& b) const {
0037       PtPxPy c;
0038       c.pt = this->pt + b.pt;
0039       c.px = this->px + b.px;
0040       c.py = this->py + b.py;
0041       return c;
0042     }
0043   };
0044 
0045   namespace Scales {
0046     const ap_fixed<12, -4> scale_degToRad = M_PI / 180.;
0047   };  // namespace Scales
0048 
0049   template <class data_T, class table_T, int N>
0050   void init_sinphi_table(table_T table_out[N]) {
0051     for (int i = 0; i < N; i++) {
0052       float x = i * (M_PI / 180.) / 2.;
0053       table_T sin_x = std::sin(x);
0054       table_out[i] = sin_x;
0055     }
0056   }
0057   template <class in_t, class table_t, int N>
0058   table_t sine_with_conversion(etaphi_t hwPhi) {
0059     table_t sin_table[N];
0060     init_sinphi_table<in_t, table_t, N>(sin_table);
0061     table_t out = sin_table[hwPhi];
0062     return out;
0063   }
0064 
0065   inline etaphi_t phi_cordic(pxy_t y, pxy_t x) {
0066 #ifdef CMSSW_GIT_HASH
0067     ap_fixed<12, 3> phi = atan2(y.to_double(), x.to_double());  // hls_math.h not available yet in CMSSW
0068 #else
0069     ap_fixed<12, 3> phi = hls::atan2(y, x);
0070 #endif
0071     ap_fixed<16, 9> etaphiscale = (float)l1ct::Scales::INTPHI_PI / M_PI;  // radians to hwPhi
0072     return phi * etaphiscale;
0073   }
0074 
0075   inline PtPxPy mht_compute(l1ct::Jet jet) {
0076     // Add an extra bit to px/py for the sign, and one additional bit to improve precision (pt_t is ap_ufixed<14, 12>)
0077     PtPxPy v_pxpy;
0078 
0079     //Initialize table once
0080     cossin_t sin_table[N_TABLE];
0081     init_sinphi_table<etaphi_t, cossin_t, N_TABLE>(sin_table);
0082 
0083     cossin_t sinphi;
0084     cossin_t cosphi;
0085     bool sign = jet.hwPhi.sign();
0086 
0087     etaphi_t hwphi = jet.hwPhi;
0088 
0089     // Reduce precision of hwPhi
0090     ap_int<10> phi;
0091     phi.V = hwphi(11, 1);
0092     phi = (phi > 0) ? phi : (ap_int<10>)-phi;  //Only store values for positive phi, pick up sign later
0093 
0094     sinphi = sin_table[phi];
0095 
0096     sinphi = (sign > 0) ? (cossin_t)(-sign * sinphi) : sinphi;  // Change sign bit if hwPt is negative, sin(-x)=-sin(x)
0097     cosphi = sin_table[phi + 90 * 2];  //cos(x)=sin(x+90). Do nothing with sign, cos(-θ) = cos θ,
0098 
0099     v_pxpy.pt = jet.hwPt;
0100     v_pxpy.py = jet.hwPt * sinphi;
0101     v_pxpy.px = jet.hwPt * cosphi;
0102 
0103     return v_pxpy;
0104   }
0105 }  // namespace P2L1HTMHTEmu
0106 
0107 //TODO replace with l1ct::Jet
0108 inline l1ct::Sum htmht(std::vector<l1ct::Jet> jets) {
0109   // compute jet px, py
0110   std::vector<P2L1HTMHTEmu::PtPxPy> ptpxpy;
0111   ptpxpy.resize(jets.size());
0112   std::transform(
0113       jets.begin(), jets.end(), ptpxpy.begin(), [](const l1ct::Jet& jet) { return P2L1HTMHTEmu::mht_compute(jet); });
0114 
0115   // Sum pt, px, py over jets
0116   P2L1HTMHTEmu::PtPxPy hthxhy = std::accumulate(ptpxpy.begin(), ptpxpy.end(), P2L1HTMHTEmu::PtPxPy());
0117 
0118   // Compute the MHT magnitude and direction
0119   l1ct::Sum ht;
0120   ht.hwSumPt = hthxhy.pt;
0121 #ifdef CMSSW_GIT_HASH
0122   ht.hwPt =
0123       sqrt(((hthxhy.px * hthxhy.px) + (hthxhy.py * hthxhy.py)).to_double());  // hls_math.h not available yet in CMSSW
0124 #else
0125   ht.hwPt = hls::sqrt(((hthxhy.px * hthxhy.px) + (hthxhy.py * hthxhy.py)));
0126 #endif
0127   ht.hwPhi = P2L1HTMHTEmu::phi_cordic(hthxhy.py, hthxhy.px);
0128   return ht;
0129 }
0130 
0131 #endif