Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-10 01:53:49

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