interface/jetmet/L1SeedConePFJetEmulator.h

0001 #ifndef L1Trigger_Phase2L1ParticleFlow_L1SeedConePFJetEmulator_h
0002 #define L1Trigger_Phase2L1ParticleFlow_L1SeedConePFJetEmulator_h
0003
0004 #define NCONSTITSFW 32  // DEFINE THE MAXIMUM NUMBER OF CONSTITUENTS USED TO CALCULATE THE JET MASS
0005
0006 #include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h"
0007 #include "DataFormats/L1TParticleFlow/interface/jets.h"
0008
0009 #include <iostream>
0010 #include <vector>
0011 #include <numeric>
0012 #include <algorithm>
0013
0014 #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
0015
0016 class L1SCJetEmu {
0017 public:
0018   // Data types and constants used in the FPGA and FPGA-optimized functions
0019   // This header file is also for use in the standalone FPGA-tools simulation
0020   // and thus contains no CMSSW/EDM specific content
0021   typedef l1ct::pt_t pt_t;
0022   typedef l1ct::glbeta_t etaphi_t;       // Type for eta & phi
0023   typedef ap_int<13> detaphi_t;          // Type for deta & dphi
0024   typedef ap_fixed<18, 23> detaphi2_t;   // Type for deta^2 & dphi^2
0025   typedef ap_fixed<22, 22> pt_etaphi_t;  // Type for product of pt with deta & dphi
0026
0027   typedef ap_ufixed<13, 1, AP_RND, AP_SAT>
0028       eventrig_t;  // stores values between 0 and 2, - 0 bit for sign, - 1 bit for integer, leaves 12 for frac
0029   typedef ap_fixed<13, 1, AP_RND, AP_SAT>
0030       oddtrig_t;  // stores values between -1 and 1, 13 - 1 bit for sign, - 0 bits for integer, leaves 12 for frac
0031
0032   // 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
0033   typedef l1ct::mass2_t mass2_t;
0034
0035   typedef ap_ufixed<20, 12, AP_TRN, AP_SAT> ppt_t;  // stores values between -1 and 1
0036   typedef ap_fixed<22, 14, AP_TRN, AP_SAT> npt_t;   // stores values between -1 and 1        JUST REDUCED BY 2
0037
0038   typedef l1ct::PuppiObjEmu Particle;
0039
0040   class Jet : public l1ct::Jet {
0041   public:
0042     std::vector<l1ct::PuppiObjEmu> constituents;
0043   };
0044
0045   L1SCJetEmu(bool debug, float coneSize, unsigned nJets);
0046
0047   std::vector<Jet> emulateEvent(std::vector<Particle>& parts) const;
0048
0049 private:
0050   // Configuration settings
0051   bool debug_;
0052   float coneSize_;
0053   unsigned nJets_;
0054   detaphi2_t rCone2_;
0055
0056   // constants for the axis update
0057   typedef ap_ufixed<18, -2> inv_pt_t;
0058   static constexpr int N_table_inv_pt = 1024;
0059   inv_pt_t inv_pt_table_[N_table_inv_pt];
0060   static constexpr int hwEtaPhi_steps = 185;  // corresponds to eta/phi range of 0 to 0.8
0061
0062   static constexpr int ceillog2(int x) { return (x <= 2) ? 1 : 1 + ceillog2((x + 1) / 2); }
0063
0064   static constexpr int floorlog2(int x) { return (x < 2) ? 0 : 1 + floorlog2(x / 2); }
0065
0066   template <int B>
0067   static constexpr int pow(int x) {
0068     return x == 0 ? 1 : B * pow<B>(x - 1);
0069   }
0070
0071   static constexpr int pow2(int x) { return pow<2>(x); }
0072
0073   /* ---
0074   * Balanced tree reduce implementation.
0075   * Reduces an array of inputs to a single value using the template binary operator 'Op',
0076   * for example summing all elements with Op_add, or finding the maximum with Op_max
0077   * Use only when the input array is fully unrolled. Or, slice out a fully unrolled section
0078   * before applying and accumulate the result over the rolled dimension.
0079   * Required for emulation to guarantee equality of ordering.
0080   * --- */
0081   template <class T, class Op>
0082   static T reduce(std::vector<T> x, Op op) {
0083     int N = x.size();
0084     int leftN = pow2(floorlog2(N - 1)) > 0 ? pow2(floorlog2(N - 1)) : 0;
0085     //static constexpr int rightN = N - leftN > 0 ? N - leftN : 0;
0086     if (N == 1) {
0087       return x.at(0);
0088     } else if (N == 2) {
0089       return op(x.at(0), x.at(1));
0090     } else {
0091       std::vector<T> left(x.begin(), x.begin() + leftN);
0092       std::vector<T> right(x.begin() + leftN, x.end());
0093       return op(reduce<T, Op>(left, op), reduce<T, Op>(right, op));
0094     }
0095   }
0096
0097   class OpPuppiObjMax {
0098   public:
0099     Particle const& operator()(Particle const& a, Particle const& b) const { return a.hwPt >= b.hwPt ? a : b; }
0100   };
0101
0102   static constexpr OpPuppiObjMax op_max{};
0103
0104   template <class data_T, int N>
0105   static inline float real_val_from_idx(unsigned i) {
0106     // Treat the index as the top N bits
0107     static constexpr int NB = ceillog2(N);  // number of address bits for table
0108     data_T x(0);
0109     // The MSB of 1 is implicit in the table
0110     x[x.width - 1] = 1;
0111     // So we can use the next NB bits for real data
0112     x(x.width - 2, x.width - NB - 1) = i;
0113     return (float)x;
0114   }
0115
0116   template <class data_T, int N>
0117   static inline unsigned idx_from_real_val(data_T x) {
0118     // Slice the top N bits to get an index into the table
0119     static constexpr int NB = ceillog2(N);  // number of address bits for table
0120     // Slice the top-1 NB bits of the value
0121     // the MSB of '1' is implicit, so only slice below that
0122     ap_uint<NB> y = x(x.width - 2, x.width - NB - 1);
0123     return (unsigned)y(NB - 1, 0);
0124   }
0125
0126   template <class data_T, class table_T, int N>
0127   static void init_invert_table(table_T table_out[N]) {
0128     // The template data_T is the data type used to address the table
0129     for (unsigned i = 0; i < N; i++) {
0130       float x = real_val_from_idx<data_T, N>(i);
0131       table_T inv_x = 1 / x;
0132       table_out[i] = inv_x;
0133     }
0134   }
0135
0136   template <class in_t, class table_t, int N>
0137   static table_t invert_with_shift(const in_t in, const table_t inv_table[N], bool debug = false) {
0138     // find the first '1' in the denominator
0139     int msb = 0;
0140     for (int b = 0; b < in.width; b++) {
0141       if (in[b])
0142         msb = b;
0143     }
0144     // shift up the denominator such that the left-most bit (msb) is '1'
0145     in_t in_shifted = in << (in.width - msb - 1);
0146     // lookup the inverse of the shifted input
0147     int idx = idx_from_real_val<in_t, N>(in_shifted);
0148     table_t inv_in = inv_table[idx];
0149     // shift the output back
0150     table_t out = inv_in << (in.width - msb - 1);
0151     if (debug) {
0152       dbgCout() << "           x " << in << ", msb = " << msb << ", shift = " << (in.width - msb) << ", idx = " << idx
0153                 << std::endl;
0154       dbgCout() << "     pre 1 / " << in_shifted << " = " << inv_in << "(" << 1 / (float)in_shifted << ")" << std::endl;
0155       dbgCout() << "    post 1 / " << in << " = " << out << "(" << 1 / (float)in << ")" << std::endl;
0156     }
0157     return out;
0158   }
0159
0160   template <typename lut_T, int N>
0161   static std::array<lut_T, N> init_trig_lut(lut_T (*func)(float)) {
0162     std::array<lut_T, N> lut;
0163     for (int hwEtaPhi = 0; hwEtaPhi < N; hwEtaPhi++) {
0164       lut[hwEtaPhi] = func(hwEtaPhi * l1ct::Scales::ETAPHI_LSB);
0165     }
0166     return lut;
0167   }
0168
0169   static detaphi_t deltaPhi(Particle a, Particle b);
0170   bool inCone(Particle seed, Particle part) const;
0171   std::vector<Particle> sortConstituents(const std::vector<Particle>& parts, const Particle seed) const;
0172   mass2_t jetMass_HW(const std::vector<Particle>& parts) const;
0173   Jet makeJet_HW(const std::vector<Particle>& parts, const Particle seed) const;
0174 };  // class L1SCJetEmu
0175
0176 #endif