File indexing completed on 2024-04-06 12:21:27
0001 #ifndef L1HPSPFTAUEMULATOR_H
0002 #define L1HPSPFTAUEMULATOR_H
0003
0004 #include "ap_int.h"
0005 #include "ap_fixed.h"
0006 #include <iostream>
0007 #include <vector>
0008 #include <numeric>
0009 #include <algorithm>
0010 #include "DataFormats/L1TParticleFlow/interface/taus.h"
0011 #include "DataFormats/L1TParticleFlow/interface/puppi.h"
0012
0013 namespace l1HPSPFTauEmu {
0014
0015
0016
0017 constexpr float etaphi_base = 720. / M_PI;
0018 constexpr float dz_base = 0.05;
0019 typedef l1ct::pt_t pt_t;
0020 typedef l1ct::glbeta_t etaphi_t;
0021
0022 typedef ap_int<13> detaphi_t;
0023
0024 typedef ap_uint<20> detaphi2_t;
0025 typedef ap_uint<5> count_t;
0026 typedef ap_uint<3> type_t;
0027 typedef l1ct::z0_t dz_t;
0028
0029 class Particle : public l1ct::PuppiObj {
0030 public:
0031 type_t pID;
0032 dz_t tempZ0;
0033 };
0034
0035 class Tau : public l1ct::Tau {
0036 public:
0037 ap_uint<5> seed_index;
0038 pt_t seed_pt;
0039 etaphi_t seed_eta;
0040 etaphi_t seed_phi;
0041 dz_t seed_z0;
0042 };
0043
0044
0045 const detaphi_t strip_phi = 0.20 * etaphi_base;
0046 const detaphi_t strip_eta = 0.05 * etaphi_base;
0047
0048 const pt_t min_leadChargedPfCand_pt = l1ct::Scales::makePtFromFloat(1.);
0049 const detaphi_t isoConeSize = 0.4 * etaphi_base;
0050 const detaphi_t delta_Rclean = 0.4 * etaphi_base;
0051 const dz_t dzCut = 0.4 / dz_base;
0052 const etaphi_t etaCutoff = 2.4 * etaphi_base;
0053
0054 template <int W, int I, ap_q_mode _AP_Q, ap_o_mode _AP_O>
0055 ap_ufixed<W, I> fp_abs(ap_fixed<W, I, _AP_Q, _AP_O> x) {
0056 ap_ufixed<W, I> result;
0057 if (x < 0) {
0058 result = -x;
0059 } else {
0060 result = x;
0061 }
0062 return result;
0063 }
0064
0065 template <int W>
0066 inline ap_uint<W> int_abs(ap_int<W> x) {
0067 ap_uint<W> result;
0068 if (x < 0) {
0069 result = -x;
0070 } else {
0071 result = x;
0072 }
0073 return result;
0074 }
0075
0076 template <class inP>
0077 inline bool is_charged(inP part) {
0078 bool charge = false;
0079 if ((part.pID == 0) || (part.pID == 1) || (part.pID == 4)) {
0080 charge = true;
0081 } else {
0082 charge = false;
0083 }
0084 return charge;
0085 }
0086
0087 inline ap_uint<20> setSConeSize2(pt_t tpt) {
0088 ap_ufixed<16, 14> total_pt = tpt * 4;
0089 ap_uint<20> SignalConeSizeSquare;
0090 if (total_pt < 115) {
0091 SignalConeSizeSquare = 23 * 23;
0092 } else if (total_pt < 120) {
0093 SignalConeSizeSquare = 22 * 22;
0094 } else if (total_pt < 126) {
0095 SignalConeSizeSquare = 21 * 21;
0096 } else if (total_pt < 132) {
0097 SignalConeSizeSquare = 20 * 20;
0098 } else if (total_pt < 139) {
0099 SignalConeSizeSquare = 19 * 19;
0100 } else if (total_pt < 147) {
0101 SignalConeSizeSquare = 18 * 18;
0102 } else if (total_pt < 156) {
0103 SignalConeSizeSquare = 17 * 17;
0104 } else if (total_pt < 166) {
0105 SignalConeSizeSquare = 16 * 16;
0106 } else if (total_pt < 178) {
0107 SignalConeSizeSquare = 15 * 15;
0108 } else if (total_pt < 191) {
0109 SignalConeSizeSquare = 14 * 14;
0110 } else if (total_pt < 206) {
0111 SignalConeSizeSquare = 13 * 13;
0112 } else if (total_pt < 224) {
0113 SignalConeSizeSquare = 12 * 12;
0114 } else {
0115 SignalConeSizeSquare = 12 * 12;
0116 }
0117 return SignalConeSizeSquare;
0118 }
0119
0120 inline bool inIsolationCone(Particle part, Particle seed) {
0121 bool inCone = false;
0122 detaphi2_t isoConeSize2 = isoConeSize * isoConeSize;
0123
0124 if (part.hwPt != 0) {
0125 detaphi_t deltaEta = part.hwEta - seed.hwEta;
0126 detaphi_t deltaPhi = part.hwPhi - seed.hwPhi;
0127 if ((deltaEta * deltaEta + deltaPhi * deltaPhi) < isoConeSize2) {
0128 inCone = true;
0129 } else {
0130 inCone = false;
0131 }
0132 } else {
0133 inCone = false;
0134 }
0135 return inCone;
0136 }
0137
0138 inline bool inSignalCone(
0139 Particle part, Particle seed, const int track_count, ap_uint<20> cone2, pt_t& iso_pt, bool& isLead) {
0140
0141
0142 bool isPotentialLead = false;
0143
0144 isPotentialLead =
0145 is_charged(part) && part.pID != 4 && part.hwPt > min_leadChargedPfCand_pt && int_abs(part.hwEta) < etaCutoff;
0146
0147
0148 bool inCone = false;
0149
0150 if (part.hwPt != 0) {
0151 detaphi_t deltaEta = part.hwEta - seed.hwEta;
0152 detaphi_t deltaPhi = part.hwPhi - seed.hwPhi;
0153 detaphi2_t deltaEta2 = deltaEta * deltaEta;
0154 detaphi2_t deltaPhi2 = deltaPhi * deltaPhi;
0155 dz_t deltaZ = 0;
0156 if (part.tempZ0 && seed.tempZ0) {
0157 deltaZ = part.tempZ0 - seed.tempZ0;
0158 }
0159
0160 if ((int_abs(deltaEta) < strip_eta) && (int_abs(deltaPhi) < strip_phi) && (part.pID == 3 || (part.pID == 1))) {
0161 if (isPotentialLead) {
0162 isLead = true;
0163 }
0164 inCone = true;
0165
0166 } else if (((deltaEta2 + deltaPhi2) < cone2) && !((part.pID == 0) && (track_count > 3)) &&
0167 !(is_charged(part) && int_abs(deltaZ) > dzCut)) {
0168 if (isPotentialLead) {
0169 isLead = true;
0170 }
0171 inCone = true;
0172 } else {
0173 if (is_charged(part) && int_abs(deltaZ) > dzCut) {
0174 iso_pt += part.hwPt;
0175 inCone = false;
0176 }
0177 }
0178 }
0179 return inCone;
0180 }
0181
0182 inline Tau makeHPSTauHW(const std::vector<Particle>& parts,
0183 const Particle seed,
0184 const pt_t total_pt ) {
0185 using namespace l1HPSPFTauEmu;
0186
0187 ap_uint<20> scone2 = setSConeSize2(total_pt);
0188
0189 pt_t isocone_pt = 0;
0190
0191 pt_t sum_pt = 0;
0192
0193 ap_fixed<22, 20> sum_eta = 0;
0194 ap_fixed<22, 20> sum_phi = 0;
0195
0196 pt_t tau_pt = 0;
0197 etaphi_t tau_eta = 0;
0198 etaphi_t tau_phi = 0;
0199
0200 pt_t chargedIsoPileup = 0;
0201 std::vector<Particle> signalParts;
0202 std::vector<Particle> outsideParts;
0203
0204 int trct = 0;
0205 bool leadCand = false;
0206 bool leadSet = false;
0207 Particle lead;
0208 for (std::vector<int>::size_type i = 0; i != parts.size(); i++) {
0209 bool isSignal = inSignalCone(parts.at(i), seed, trct, scone2, isocone_pt, leadCand);
0210 if (isSignal) {
0211 signalParts.push_back(parts.at(i));
0212 if (parts[i].pID == 0) {
0213 trct++;
0214 }
0215 if (leadCand) {
0216 if (leadSet == false) {
0217 lead = parts[i];
0218 leadCand = false;
0219 leadSet = true;
0220 } else {
0221 if (parts[i].hwPt > lead.hwPt) {
0222 lead = parts[i];
0223 leadCand = false;
0224 } else {
0225 leadCand = false;
0226 }
0227 }
0228 }
0229 } else {
0230 outsideParts.push_back(parts.at(i));
0231 }
0232 }
0233
0234 for (std::vector<int>::size_type i = 0; i != signalParts.size(); i++) {
0235 Particle sigP = signalParts.at(i);
0236 if (is_charged(sigP) || (sigP.pID == 3)) {
0237 sum_pt += sigP.hwPt;
0238 sum_eta += sigP.hwPt * sigP.hwEta;
0239 sum_phi += sigP.hwPt * sigP.hwPhi;
0240 }
0241 }
0242
0243 pt_t div_pt = 1;
0244 if (sum_pt == 0) {
0245 div_pt = 1;
0246 } else {
0247 div_pt = sum_pt;
0248 }
0249
0250 tau_pt = sum_pt;
0251 tau_eta = sum_eta / div_pt;
0252 tau_phi = sum_phi / div_pt;
0253
0254 if (tau_pt > l1ct::Scales::makePtFromFloat(20.) && int_abs(tau_eta) < etaCutoff && leadSet == true &&
0255 isocone_pt < tau_pt) {
0256 Tau tau;
0257 tau.hwPt = tau_pt;
0258 tau.hwEta = tau_eta;
0259 tau.hwPhi = tau_phi;
0260 return tau;
0261 } else {
0262 Tau tau;
0263 tau.hwPt = 0.;
0264 tau.hwEta = 0.;
0265 tau.hwPhi = 0.;
0266 return tau;
0267 }
0268 }
0269
0270 inline std::vector<Tau> emulateEvent(std::vector<Particle>& parts, std::vector<Particle>& jets, bool jEnable) {
0271 using namespace l1HPSPFTauEmu;
0272
0273 std::vector<Particle> parts_copy;
0274 parts_copy.resize(parts.size());
0275 std::transform(parts.begin(), parts.end(), parts_copy.begin(), [](const Particle& part) { return part; });
0276
0277 std::sort(
0278 parts_copy.begin(), parts_copy.end(), [](const Particle& i, const Particle& j) { return (i.hwPt > j.hwPt); });
0279
0280
0281 std::vector<Particle> jets_copy;
0282 jets_copy.resize(jets.size());
0283 std::transform(jets.begin(), jets.end(), jets_copy.begin(), [](const Particle& jet) { return jet; });
0284 std::sort(
0285 jets_copy.begin(), jets_copy.end(), [](const Particle& i, const Particle& j) { return (i.hwPt > j.hwPt); });
0286
0287 std::vector<Tau> taus;
0288 std::vector<Tau> cleaned_taus;
0289 taus.reserve(20);
0290 std::vector<Particle> preseed;
0291 preseed.reserve(144);
0292
0293
0294
0295 int jets_index = 0;
0296 int jets_max = jets_copy.size();
0297
0298 std::vector<Particle> jseeds;
0299 jseeds.reserve(4);
0300
0301 int parts_index = 0;
0302 int parts_max = parts_copy.size();
0303
0304 while (preseed.size() < 128 && parts_index != parts_max) {
0305 Particle pSeed = parts_copy.at(parts_index);
0306
0307 if (pSeed.hwPt > l1ct::Scales::makePtFromFloat(5.) && is_charged(pSeed) && int_abs(pSeed.hwEta) < etaCutoff) {
0308 preseed.push_back(pSeed);
0309 }
0310
0311 parts_index++;
0312 }
0313
0314 std::vector<Particle> seeds;
0315 seeds.reserve(16);
0316 std::vector<int>::size_type pseed_index = 0;
0317 while (seeds.size() < 16 && pseed_index < preseed.size()) {
0318 seeds.push_back(preseed.at(pseed_index));
0319 pseed_index++;
0320 }
0321
0322
0323 if (jEnable) {
0324 while (jseeds.size() < 4 && jets_index != jets_max) {
0325 Particle jSeed = jets_copy.at(jets_index);
0326
0327 if (jSeed.hwPt > l1ct::Scales::makePtFromFloat(20.) && int_abs(jSeed.hwEta) < etaCutoff) {
0328 jseeds.push_back(jSeed);
0329 }
0330 jets_index++;
0331 }
0332 }
0333 for (std::vector<int>::size_type i = 0; i != seeds.size(); i++) {
0334 Particle seed = seeds[i];
0335
0336 std::vector<Particle> iso_parts;
0337
0338 iso_parts.reserve(30);
0339 pt_t total_pt = 0;
0340 std::vector<int>::size_type iso_index = 0;
0341 while (iso_index < parts_copy.size() && iso_parts.size() < 30) {
0342 Particle isoCand = parts_copy.at(iso_index);
0343 if (inIsolationCone(isoCand, seed)) {
0344 iso_parts.push_back(isoCand);
0345 total_pt += isoCand.hwPt;
0346 }
0347 iso_index++;
0348 }
0349
0350 taus.push_back(makeHPSTauHW(iso_parts, seed, total_pt));
0351 }
0352
0353
0354 if (jEnable) {
0355 for (std::vector<int>::size_type i = 0; i != jseeds.size(); i++) {
0356 Particle jseed = jseeds[i];
0357 std::vector<Particle> iso_parts;
0358 iso_parts.reserve(30);
0359 pt_t total_pt = 0;
0360 std::vector<int>::size_type iso_index = 0;
0361
0362 pt_t max_pt_j = 0;
0363 while (iso_index < parts_copy.size()) {
0364 Particle isoCand = parts_copy.at(iso_index);
0365
0366 if (inIsolationCone(isoCand, jseed)) {
0367 if (is_charged(isoCand) && isoCand.hwPt > max_pt_j) {
0368 if (isoCand.tempZ0) {
0369 jseed.tempZ0 = isoCand.tempZ0;
0370 }
0371 max_pt_j = isoCand.hwPt;
0372 }
0373
0374 if (iso_parts.size() < 30) {
0375 iso_parts.push_back(isoCand);
0376 total_pt += isoCand.hwPt;
0377 }
0378 }
0379 iso_index++;
0380 }
0381 taus.push_back(makeHPSTauHW(iso_parts, jseed, total_pt));
0382 }
0383 }
0384
0385 std::sort(taus.begin(), taus.end(), [](const Tau& i, const Tau& j) { return (i.hwPt > j.hwPt); });
0386
0387 int taus_max = taus.size();
0388
0389 bool matrix[380];
0390
0391 for (int i = 0; i < (taus_max - 1); i++) {
0392 for (int j = i + 1; j < taus_max; j++) {
0393 etaphi_t deltaE = taus[i].hwEta - taus[j].hwEta;
0394 etaphi_t deltaP = taus[i].hwPhi - taus[j].hwPhi;
0395 if ((deltaE * deltaE + deltaP * deltaP) < (delta_Rclean * delta_Rclean)) {
0396 matrix[i * 19 + j] = true;
0397 } else {
0398 matrix[i * 19 + j] = false;
0399 }
0400 }
0401 }
0402
0403 if (!taus.empty()) {
0404 if (taus[0].hwPt > 0) {
0405 cleaned_taus.push_back(taus.at(0));
0406 }
0407
0408 bool clean[20];
0409
0410 for (int i = 0; i < 20; i++) {
0411 clean[i] = false;
0412 }
0413 for (int i = 1; i < taus_max; i++) {
0414 for (int j = i - 1; j >= 0; j--) {
0415 clean[i] |= (matrix[j * 19 + i] && !clean[j]);
0416 }
0417 if (!clean[i] && taus[i].hwPt > 0) {
0418 cleaned_taus.push_back(taus.at(i));
0419 }
0420 }
0421 }
0422
0423 return cleaned_taus;
0424 }
0425
0426 };
0427
0428 #endif