L1TTrackMatch/interface/L1TkHTMissEmulatorProducer.h

0001 #ifndef L1Trigger_L1TTrackMatch_L1TkHTMissEmulatorProducer_HH
0002 #define L1Trigger_L1TTrackMatch_L1TkHTMissEmulatorProducer_HH
0003
0004 // Original Author:  Hardik Routray
0005 //         Created:  Mon, 11 Oct 2021
0006
0007 #include <ap_int.h>
0008
0009 #include <cmath>
0010 #include <cstdint>
0011 #include <filesystem>
0012 #include <fstream>
0013 #include <iostream>
0014 #include <numeric>
0015
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "DataFormats/L1Trigger/interface/TkJetWord.h"
0018
0019 // Namespace that defines constants and types used by the HTMiss Emulation
0020
0021 namespace l1tmhtemu {
0022
0023   const unsigned int kInternalPtWidth{l1t::TkJetWord::TkJetBitWidths::kPtSize};
0024   const unsigned int kInternalEtaWidth{l1t::TkJetWord::TkJetBitWidths::kGlbEtaSize};
0025   const unsigned int kInternalPhiWidth{l1t::TkJetWord::TkJetBitWidths::kGlbPhiSize};
0026
0027   // extra room for sumPx, sumPy
0028   const unsigned int kEtExtra{10};
0029   const unsigned int kValidSize{1};
0030   const unsigned int kMHTSize{16};  // For output Magnitude default 16
0031   const unsigned int kMHTIntSize{11};
0032   const unsigned int kMHTPhiSize{13};  // For output Phi default 13
0033   const unsigned int kHTSize{kInternalPtWidth + kEtExtra};
0034   const unsigned int kUnassignedSize{64 - (kHTSize + kMHTSize + kMHTPhiSize + kValidSize)};
0035
0036   enum BitLocations {
0037     // The location of the least significant bit (LSB) and most significant bit (MSB) in the sum word for different fields
0038     kValidLSB = 0,
0039     kValidMSB = kValidLSB + kValidSize - 1,
0040     kMHTLSB = kValidMSB + 1,
0041     kMHTMSB = kMHTLSB + kMHTSize - 1,
0042     kMHTPhiLSB = kMHTMSB + 1,
0043     kMHTPhiMSB = kMHTPhiLSB + kMHTPhiSize - 1,
0044     kHTLSB = kMHTPhiMSB + 1,
0045     kHTMSB = kHTLSB + kHTSize - 1,
0046     kUnassignedLSB = kHTMSB + 1,
0047     kUnassignedMSB = kUnassignedLSB + kUnassignedSize - 1
0048   };
0049
0050   const float kMaxMHT{2048};  // 2 TeV
0051   const float kMaxMHTPhi{2 * M_PI};
0052
0053   typedef ap_uint<5> ntracks_t;
0054   typedef ap_uint<kInternalPtWidth> pt_t;
0055   typedef ap_int<kInternalEtaWidth> eta_t;
0056   typedef ap_int<kInternalPhiWidth> phi_t;
0057
0058   typedef ap_int<kHTSize> Et_t;
0059   typedef ap_ufixed<kMHTSize, kMHTIntSize> MHT_t;
0060   typedef ap_uint<kMHTPhiSize> MHTphi_t;
0061
0062   const unsigned int kMHTBins = 1 << kMHTSize;
0063   const unsigned int kMHTPhiBins = 1 << kMHTPhiSize;
0064
0065   const double kStepPt{0.25};
0066   const double kStepEta{M_PI / (720)};
0067   const double kStepPhi{M_PI / (720)};
0068
0069   const double kStepMHT = (l1tmhtemu::kMaxMHT / l1tmhtemu::kMHTBins);
0070   const double kStepMHTPhi = (l1tmhtemu::kMaxMHTPhi / l1tmhtemu::kMHTPhiBins);
0071
0072   const unsigned int kPhiBins = 1 << kInternalPhiWidth;
0073
0074   const float kMaxCosLUTPhi{M_PI};
0075
0076   template <typename T>
0077   T digitizeSignedValue(double value, unsigned int nBits, double lsb) {
0078     T digitized_value = std::floor(std::abs(value) / lsb);
0079     T digitized_maximum = (1 << (nBits - 1)) - 1;  // The remove 1 bit from nBits to account for the sign
0080     if (digitized_value > digitized_maximum)
0081       digitized_value = digitized_maximum;
0082     if (value < 0)
0083       digitized_value = (1 << nBits) - digitized_value;  // two's complement encoding
0084     return digitized_value;
0085   }
0086
0087   inline std::vector<phi_t> generateCosLUT(unsigned int size) {  // Fill cosine LUT with integer values
0088     float phi = 0;
0089     std::vector<phi_t> cosLUT;
0090     for (unsigned int LUT_idx = 0; LUT_idx < size; LUT_idx++) {
0091       cosLUT.push_back(digitizeSignedValue<phi_t>(cos(phi), l1tmhtemu::kInternalPhiWidth, l1tmhtemu::kStepPhi));
0092       phi += l1tmhtemu::kStepPhi;
0093     }
0094     cosLUT.push_back((phi_t)(0));  //Prevent overflow in last bin
0095     return cosLUT;
0096   }
0097
0098   inline std::vector<MHTphi_t> generateaTanLUT(int cordicSteps) {  // Fill atan LUT with integer values
0099     std::vector<MHTphi_t> atanLUT;
0100     atanLUT.reserve(cordicSteps);
0101     for (int cordicStep = 0; cordicStep < cordicSteps; cordicStep++) {
0102       atanLUT.push_back(MHTphi_t(round((kMHTPhiBins * atan(pow(2, -1 * cordicStep))) / (2 * M_PI))));
0103     }
0104     return atanLUT;
0105   }
0106
0107   inline std::vector<Et_t> generatemagNormalisationLUT(int cordicSteps) {
0108     float val = 1.0;
0109     std::vector<Et_t> magNormalisationLUT;
0110     for (int cordicStep = 0; cordicStep < cordicSteps; cordicStep++) {
0111       val = val / (pow(1 + pow(4, -1 * cordicStep), 0.5));
0112       magNormalisationLUT.push_back(Et_t(round(kMHTBins * val)));
0113     }
0114     return magNormalisationLUT;
0115   }
0116
0117   struct EtMiss {
0118     MHT_t Et;
0119     MHTphi_t Phi;
0120   };
0121
0122   inline EtMiss cordicSqrt(Et_t x,
0123                            Et_t y,
0124                            int cordicSteps,
0125                            std::vector<l1tmhtemu::MHTphi_t> atanLUT,
0126                            std::vector<Et_t> magNormalisationLUT) {
0127     Et_t new_x = 0;
0128     Et_t new_y = 0;
0129
0130     MHTphi_t phi = 0;
0131     MHTphi_t new_phi = 0;
0132     bool sign = false;
0133
0134     EtMiss ret_etmiss;
0135
0136     if (x >= 0 && y >= 0) {
0137       phi = 0;
0138       sign = true;
0139       //x = x;
0140       //y = y;
0141     } else if (x < 0 && y >= 0) {
0142       phi = kMHTPhiBins >> 1;
0143       sign = false;
0144       x = -x;
0145       //y = y;
0146     } else if (x < 0 && y < 0) {
0147       phi = kMHTPhiBins >> 1;
0148       sign = true;
0149       x = -x;
0150       y = -y;
0151     } else {
0152       phi = kMHTPhiBins;
0153       sign = false;
0154       //x = x;
0155       y = -y;
0156     }
0157
0158     for (int step = 0; step < cordicSteps; step++) {
0159       if (y < 0) {
0160         new_x = x - (y >> step);
0161         new_y = y + (x >> step);
0162       } else {
0163         new_x = x + (y >> step);
0164         new_y = y - (x >> step);
0165       }
0166
0167       if ((y < 0) == sign) {
0168         new_phi = phi - atanLUT[step];
0169       } else {
0170         new_phi = phi + atanLUT[step];
0171       }
0172
0173       x = new_x;
0174       y = new_y;
0175       phi = new_phi;
0176     }
0177
0178     float sqrtval = (float(x * magNormalisationLUT[cordicSteps - 1]) / (float)kMHTBins) * float(kStepPt * kStepPhi);
0179     ret_etmiss.Et = sqrtval;
0180     ret_etmiss.Phi = phi;
0181
0182     return ret_etmiss;
0183   }
0184
0185 }  // namespace l1tmhtemu
0186 #endif