interface/taus/L1HPSPFTauEmulator.h

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   //mapping the eta phi onto bits
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;  //type for multiplicity
0026   typedef ap_uint<3> type_t;   //type for particle type
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   //constants
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     //finds the signal cone candidates (including strip pt check
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     //calculate the deta and dphi
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 /*, Config config*/) {
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     //sorting by pt
0277     std::sort(
0278         parts_copy.begin(), parts_copy.end(), [](const Particle& i, const Particle& j) { return (i.hwPt > j.hwPt); });
0279
0280     //sorting jets by pt
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     //jet seeds reserve
0294     //4 for now
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     //first find the seeds
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);  //up to 16 track + 4 jet seeds right now
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     //With jets
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     //add in the jet taus
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 };  // namespace l1HPSPFTauEmu
0427
0428 #endif