src/jetmet/L1SeedConePFJetEmulator.cc

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1SeedConePFJetEmulator.h"
0002
0003 L1SCJetEmu::L1SCJetEmu(bool debug, float coneSize, unsigned nJets)
0004     : debug_(debug),
0005       coneSize_(coneSize),
0006       nJets_(nJets),
0007       rCone2_(coneSize * coneSize / l1ct::Scales::ETAPHI_LSB / l1ct::Scales::ETAPHI_LSB) {
0008   init_invert_table<pt_t, inv_pt_t, N_table_inv_pt>(inv_pt_table_);
0009 }
0010
0011 L1SCJetEmu::detaphi_t L1SCJetEmu::deltaPhi(L1SCJetEmu::Particle a, L1SCJetEmu::Particle b) {
0012   detaphi_t dphi = detaphi_t(a.hwPhi) - detaphi_t(b.hwPhi);
0013   // phi wrap
0014   detaphi_t dphi0 =
0015       dphi > detaphi_t(l1ct::Scales::INTPHI_PI) ? detaphi_t(l1ct::Scales::INTPHI_TWOPI - dphi) : detaphi_t(dphi);
0016   detaphi_t dphi1 =
0017       dphi < detaphi_t(-l1ct::Scales::INTPHI_PI) ? detaphi_t(l1ct::Scales::INTPHI_TWOPI + dphi) : detaphi_t(dphi);
0018   detaphi_t dphiw = dphi > detaphi_t(0) ? dphi0 : dphi1;
0019   return dphiw;
0020 }
0021
0022 bool L1SCJetEmu::inCone(L1SCJetEmu::Particle seed, L1SCJetEmu::Particle part) const {
0023   // scale the particle eta, phi to hardware units
0024   detaphi_t deta = detaphi_t(seed.hwEta) - detaphi_t(part.hwEta);
0025   detaphi_t dphi = deltaPhi(seed, part);
0026   bool ret = deta * deta + dphi * dphi < rCone2_;
0027   //bool ret = r2 < cone2;
0028   if (debug_) {
0029     detaphi2_t r2 = detaphi2_t(deta) * detaphi2_t(deta) + detaphi2_t(dphi) * detaphi2_t(dphi);
0030     dbgCout() << "  part eta, seed eta: " << part.hwEta << ", " << seed.hwEta << std::endl;
0031     dbgCout() << "  part phi, seed phi: " << part.hwPhi << ", " << seed.hwPhi << std::endl;
0032     dbgCout() << "  pt, deta, dphi, r2, cone2, lt: " << part.hwPt << ", " << deta << ", " << dphi << ", "
0033               << deta * deta + dphi * dphi << ", " << rCone2_ << ", " << ret << std::endl;
0034   }
0035   return ret;
0036 }
0037
0038 L1SCJetEmu::Jet L1SCJetEmu::makeJet_HW(const std::vector<Particle>& parts) const {
0039   // Seed Cone Jet algorithm with ap_fixed types and hardware emulation
0040   Particle seed = reduce(parts, op_max);
0041
0042   // Event with saturation, order of terms doesn't matter since they're all positive
0043   auto sumpt = [](pt_t(a), const Particle& b) { return a + b.hwPt; };
0044
0045   // Sum the pt
0046   pt_t pt = std::accumulate(parts.begin(), parts.end(), pt_t(0), sumpt);
0047   inv_pt_t inv_pt = invert_with_shift<pt_t, inv_pt_t, N_table_inv_pt>(pt, inv_pt_table_, false);
0048
0049   // pt weighted d eta
0050   std::vector<pt_etaphi_t> pt_deta;
0051   pt_deta.resize(parts.size());
0052   std::transform(parts.begin(), parts.end(), pt_deta.begin(), [&seed](const Particle& part) {
0053     // In the firmware we calculate the per-particle pt-weighted deta
0054     return pt_etaphi_t(part.hwPt * detaphi_t(part.hwEta - seed.hwEta));
0055   });
0056   // Accumulate the pt-weighted etas. Init to 0, include seed in accumulation
0057   pt_etaphi_t sum_pt_eta = std::accumulate(pt_deta.begin(), pt_deta.end(), pt_etaphi_t(0));
0058   etaphi_t eta = seed.hwEta + etaphi_t(sum_pt_eta * inv_pt);
0059
0060   // pt weighted d phi
0061   std::vector<pt_etaphi_t> pt_dphi;
0062   pt_dphi.resize(parts.size());
0063   std::transform(parts.begin(), parts.end(), pt_dphi.begin(), [&seed](const Particle& part) {
0064     // In the firmware we calculate the per-particle pt-weighted dphi
0065     return pt_etaphi_t(part.hwPt * deltaPhi(part, seed));
0066   });
0067   // Accumulate the pt-weighted phis. Init to 0, include seed in accumulation
0068   pt_etaphi_t sum_pt_phi = std::accumulate(pt_dphi.begin(), pt_dphi.end(), pt_etaphi_t(0));
0069   etaphi_t phi = seed.hwPhi + etaphi_t(sum_pt_phi * inv_pt);
0070
0071   Jet jet;
0072   jet.hwPt = pt;
0073   jet.hwEta = eta;
0074   jet.hwPhi = phi;
0075   jet.constituents = parts;
0076
0077   if (debug_) {
0078     std::for_each(pt_dphi.begin(), pt_dphi.end(), [](pt_etaphi_t& x) { dbgCout() << "pt_dphi: " << x << std::endl; });
0079     std::for_each(pt_deta.begin(), pt_deta.end(), [](pt_etaphi_t& x) { dbgCout() << "pt_deta: " << x << std::endl; });
0080     dbgCout() << " sum_pt_eta: " << sum_pt_eta << ", 1/pt: " << inv_pt
0081               << ", sum_pt_eta * 1/pt: " << etaphi_t(sum_pt_eta * inv_pt) << std::endl;
0082     dbgCout() << " sum_pt_phi: " << sum_pt_phi << ", 1/pt: " << inv_pt
0083               << ", sum_pt_phi * 1/pt: " << etaphi_t(sum_pt_phi * inv_pt) << std::endl;
0084     dbgCout() << " uncorr eta: " << seed.hwEta << ", phi: " << seed.hwPhi << std::endl;
0085     dbgCout() << "   corr eta: " << eta << ", phi: " << phi << std::endl;
0086     dbgCout() << "         pt: " << pt << std::endl;
0087   }
0088
0089   return jet;
0090 }
0091
0092 std::vector<L1SCJetEmu::Jet> L1SCJetEmu::emulateEvent(std::vector<Particle>& parts) const {
0093   // The fixed point algorithm emulation
0094   std::vector<Particle> work;
0095   work.resize(parts.size());
0096   std::transform(parts.begin(), parts.end(), work.begin(), [](const Particle& part) { return part; });
0097
0098   std::vector<Jet> jets;
0099   jets.reserve(nJets_);
0100   while (!work.empty() && jets.size() < nJets_) {
0101     // Take the highest pt candidate as a seed
0102     // Use the firmware reduce function to find the same seed as the firmware
0103     // in case there are multiple seeds with the same pT
0104     Particle seed = reduce(work, op_max);
0105
0106     // Get the particles within a coneSize_ of the seed
0107     std::vector<Particle> particlesInCone;
0108     std::copy_if(work.begin(), work.end(), std::back_inserter(particlesInCone), [&](const Particle& part) {
0109       return inCone(seed, part);
0110     });
0111     if (debug_) {
0112       dbgCout() << "Seed: " << seed.hwPt << ", " << seed.hwEta << ", " << seed.hwPhi << std::endl;
0113       std::for_each(particlesInCone.begin(), particlesInCone.end(), [&](Particle& part) {
0114         dbgCout() << "  Part: " << part.hwPt << ", " << part.hwEta << ", " << part.hwPhi << std::endl;
0115         inCone(seed, part);
0116       });
0117     }
0118     jets.push_back(makeJet_HW(particlesInCone));
0119     // remove the clustered particles
0120     work.erase(std::remove_if(work.begin(), work.end(), [&](const Particle& part) { return inCone(seed, part); }),
0121                work.end());
0122   }
0123   return jets;
0124 }