File indexing completed on 2021-07-28 03:10:07
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 #include <iostream>
0014 #include <iomanip>
0015 #include <iterator>
0016 #include <cmath>
0017 #include <map>
0018
0019 #include "L1Trigger/L1TTwinMux/interface/RPCHitCleaner.h"
0020 #include "L1Trigger/L1TTwinMux/interface/RPCtoDTTranslator.h"
0021 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0022 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0023 #include "CondFormats/DataRecord/interface/RPCEMapRcd.h"
0024 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0025
0026 using namespace std;
0027
0028 RPCtoDTTranslator::RPCtoDTTranslator(RPCDigiCollection const& inrpcDigis) : m_rpcDigis{inrpcDigis} {}
0029
0030 namespace {
0031 constexpr int max_rpc_bx = 2;
0032 constexpr int min_rpc_bx = -2;
0033
0034 struct rpc_hit {
0035 int bx;
0036 int station;
0037 int sector;
0038 int wheel;
0039 RPCDetId detid;
0040 int strip;
0041 int roll;
0042 int layer;
0043 rpc_hit(int pbx, int pstation, int psector, int pwheel, RPCDetId pdet, int pstrip, int proll, int player)
0044 : bx(pbx),
0045 station(pstation),
0046 sector(psector),
0047 wheel(pwheel),
0048 detid(pdet),
0049 strip(pstrip),
0050 roll(proll),
0051 layer(player) {}
0052 };
0053
0054
0055
0056 class BxToHit {
0057 public:
0058 BxToHit() : m_hits{} {}
0059
0060 static bool outOfRange(int iBX) { return (iBX > max_rpc_bx or iBX < min_rpc_bx); }
0061
0062 int& operator[](int iBX) { return m_hits[iBX - min_rpc_bx]; }
0063
0064 size_t size() const { return m_hits.size(); }
0065
0066 private:
0067 std::array<int, max_rpc_bx - min_rpc_bx + 1> m_hits;
0068 };
0069 }
0070
0071 void RPCtoDTTranslator::run(const RPCGeometry& rpcGeometry) {
0072 std::vector<L1MuDTChambPhDigi> l1ttma_out;
0073 std::vector<L1MuDTChambPhDigi> l1ttma_hits_out;
0074
0075 std::vector<rpc_hit> vrpc_hit_layer1, vrpc_hit_layer2, vrpc_hit_st3, vrpc_hit_st4;
0076
0077
0078 for (auto chamber = m_rpcDigis.begin(); chamber != m_rpcDigis.end(); ++chamber) {
0079 RPCDetId detid = (*chamber).first;
0080 for (auto digi = (*chamber).second.first; digi != (*chamber).second.second; ++digi) {
0081 if (detid.region() != 0)
0082 continue;
0083 if (BxToHit::outOfRange(digi->bx()))
0084 continue;
0085 if (detid.layer() == 1)
0086 vrpc_hit_layer1.emplace_back(digi->bx(),
0087 detid.station(),
0088 detid.sector(),
0089 detid.ring(),
0090 detid,
0091 digi->strip(),
0092 detid.roll(),
0093 detid.layer());
0094 if (detid.station() == 3)
0095 vrpc_hit_st3.emplace_back(digi->bx(),
0096 detid.station(),
0097 detid.sector(),
0098 detid.ring(),
0099 detid,
0100 digi->strip(),
0101 detid.roll(),
0102 detid.layer());
0103 if (detid.layer() == 2)
0104 vrpc_hit_layer2.emplace_back(digi->bx(),
0105 detid.station(),
0106 detid.sector(),
0107 detid.ring(),
0108 detid,
0109 digi->strip(),
0110 detid.roll(),
0111 detid.layer());
0112 if (detid.station() == 4)
0113 vrpc_hit_st4.emplace_back(digi->bx(),
0114 detid.station(),
0115 detid.sector(),
0116 detid.ring(),
0117 detid,
0118 digi->strip(),
0119 detid.roll(),
0120 detid.layer());
0121 }
0122 }
0123
0124 vector<int> vcluster_size;
0125 int cluster_id = -1;
0126 int itr = 0;
0127
0128 std::map<RPCHitCleaner::detId_Ext, int> hits;
0129 int cluster_size = 0;
0130 for (auto chamber = m_rpcDigis.begin(); chamber != m_rpcDigis.end(); ++chamber) {
0131 RPCDetId detid = (*chamber).first;
0132 int strip_n1 = -10000;
0133 int bx_n1 = -10000;
0134 if (detid.region() != 0)
0135 continue;
0136 for (auto digi = (*chamber).second.first; digi != (*chamber).second.second; ++digi) {
0137 if (fabs(digi->bx()) > 3)
0138 continue;
0139
0140
0141 if (abs(digi->strip() - strip_n1) != 1 || digi->bx() != bx_n1) {
0142 if (itr != 0)
0143 vcluster_size.push_back(cluster_size);
0144 cluster_size = 0;
0145 cluster_id++;
0146 }
0147 itr++;
0148 cluster_size++;
0149
0150
0151 RPCHitCleaner::detId_Ext tmp{detid, digi->bx(), digi->strip()};
0152 hits[tmp] = cluster_id;
0153
0154 strip_n1 = digi->strip();
0155 bx_n1 = digi->bx();
0156 }
0157 }
0158 vcluster_size.push_back(cluster_size);
0159
0160 for (int wh = -2; wh <= 2; wh++) {
0161 for (int sec = 1; sec <= 12; sec++) {
0162 for (int st = 1; st <= 4; st++) {
0163 int rpcbx = 0;
0164 std::vector<int> delta_phib;
0165 bool found_hits = false;
0166 std::vector<int> rpc2dt_phi, rpc2dt_phib;
0167
0168 int itr1 = 0;
0169 for (unsigned int l1 = 0; l1 < vrpc_hit_layer1.size(); l1++) {
0170 RPCHitCleaner::detId_Ext tmp{vrpc_hit_layer1[l1].detid, vrpc_hit_layer1[l1].bx, vrpc_hit_layer1[l1].strip};
0171 int id = hits[tmp];
0172 int phi1 = radialAngle(vrpc_hit_layer1[l1].detid, rpcGeometry, vrpc_hit_layer1[l1].strip);
0173 if (vcluster_size[id] == 2 && itr1 == 0) {
0174 itr1++;
0175 continue;
0176 }
0177 if (vcluster_size[id] == 2 && itr1 == 1) {
0178 itr1 = 0;
0179 phi1 = phi1 + (radialAngle(vrpc_hit_layer1[l1 - 1].detid, rpcGeometry, vrpc_hit_layer1[l1 - 1].strip));
0180 phi1 /= 2;
0181 }
0182 int itr2 = 0;
0183 for (unsigned int l2 = 0; l2 < vrpc_hit_layer2.size(); l2++) {
0184 if (vrpc_hit_layer1[l1].station != st || vrpc_hit_layer2[l2].station != st)
0185 continue;
0186 if (vrpc_hit_layer1[l1].sector != sec || vrpc_hit_layer2[l2].sector != sec)
0187 continue;
0188 if (vrpc_hit_layer1[l1].wheel != wh || vrpc_hit_layer2[l2].wheel != wh)
0189 continue;
0190 if (vrpc_hit_layer1[l1].bx != vrpc_hit_layer2[l2].bx)
0191 continue;
0192 RPCHitCleaner::detId_Ext tmp{vrpc_hit_layer2[l2].detid, vrpc_hit_layer2[l2].bx, vrpc_hit_layer2[l2].strip};
0193 int id = hits[tmp];
0194
0195 if (vcluster_size[id] == 2 && itr2 == 0) {
0196 itr2++;
0197 continue;
0198 }
0199
0200
0201 int phi2 = radialAngle(vrpc_hit_layer2[l2].detid, rpcGeometry, vrpc_hit_layer2[l2].strip);
0202 if (vcluster_size[id] == 2 && itr2 == 1) {
0203 itr2 = 0;
0204 phi2 = phi2 + (radialAngle(vrpc_hit_layer2[l2 - 1].detid, rpcGeometry, vrpc_hit_layer2[l2 - 1].strip));
0205 phi2 /= 2;
0206 }
0207 int average = ((phi1 + phi2) / 2) << 2;
0208 rpc2dt_phi.push_back(average);
0209
0210
0211
0212
0213 int xin = localXX((phi1 << 2), 1, vrpc_hit_layer1[l1].station);
0214 int xout = localXX((phi2 << 2), 2, vrpc_hit_layer2[l2].station);
0215 if (vcluster_size[id] == 2 && itr2 == 1) {
0216 int phi1_n1 = radialAngle(vrpc_hit_layer1[l1 - 1].detid, rpcGeometry, vrpc_hit_layer1[l1 - 1].strip);
0217 int phi2_n1 = radialAngle(vrpc_hit_layer2[l2 - 1].detid, rpcGeometry, vrpc_hit_layer2[l2 - 1].strip);
0218 xin += localXX((phi1_n1 << 2), 1, vrpc_hit_layer1[l1].station);
0219 xout += localXX((phi2_n1 << 2), 2, vrpc_hit_layer2[l2].station);
0220 xin /= 2;
0221 xout /= 2;
0222 }
0223
0224 int phi_b = bendingAngle(xin, xout, average);
0225
0226 rpc2dt_phib.push_back(phi_b);
0227
0228 delta_phib.push_back(abs(phi_b));
0229 found_hits = true;
0230 rpcbx = vrpc_hit_layer1[l1].bx;
0231 }
0232 }
0233 if (found_hits) {
0234
0235 int min_index = std::distance(delta_phib.begin(), std::min_element(delta_phib.begin(), delta_phib.end())) + 0;
0236
0237 l1ttma_out.emplace_back(rpcbx, wh, sec - 1, st, rpc2dt_phi[min_index], rpc2dt_phib[min_index], 3, 0, 0, 2);
0238 }
0239
0240 BxToHit hit;
0241 itr1 = 0;
0242 for (unsigned int l1 = 0; l1 < vrpc_hit_layer1.size(); l1++) {
0243 if (vrpc_hit_layer1[l1].station != st || st > 2 || vrpc_hit_layer1[l1].sector != sec ||
0244 vrpc_hit_layer1[l1].wheel != wh)
0245 continue;
0246
0247 RPCHitCleaner::detId_Ext tmp{vrpc_hit_layer1[l1].detid, vrpc_hit_layer1[l1].bx, vrpc_hit_layer1[l1].strip};
0248 int id = hits[tmp];
0249 if (vcluster_size[id] == 2 && itr1 == 0) {
0250 itr1++;
0251 continue;
0252 }
0253 int phi2 = radialAngle(vrpc_hit_layer1[l1].detid, rpcGeometry, vrpc_hit_layer1[l1].strip);
0254 phi2 = phi2 << 2;
0255 if (vcluster_size[id] == 2 && itr1 == 1) {
0256 itr1 = 0;
0257 phi2 = phi2 + (radialAngle(vrpc_hit_layer1[l1 - 1].detid, rpcGeometry, vrpc_hit_layer1[l1 - 1].strip) << 2);
0258 phi2 /= 2;
0259 }
0260
0261 l1ttma_hits_out.emplace_back(
0262 vrpc_hit_layer1[l1].bx, wh, sec - 1, st, phi2, 0, 3, hit[vrpc_hit_layer1[l1].bx], 0, 2);
0263 hit[vrpc_hit_layer1[l1].bx]++;
0264 }
0265 itr1 = 0;
0266 for (unsigned int l2 = 0; l2 < vrpc_hit_layer2.size(); l2++) {
0267 if (vrpc_hit_layer2[l2].station != st || st > 2 || vrpc_hit_layer2[l2].sector != sec ||
0268 vrpc_hit_layer2[l2].wheel != wh)
0269 continue;
0270 RPCHitCleaner::detId_Ext tmp{vrpc_hit_layer2[l2].detid, vrpc_hit_layer2[l2].bx, vrpc_hit_layer2[l2].strip};
0271 int id = hits[tmp];
0272
0273 if (vcluster_size[id] == 2 && itr1 == 0) {
0274 itr1++;
0275 continue;
0276 }
0277 int phi2 = radialAngle(vrpc_hit_layer2[l2].detid, rpcGeometry, vrpc_hit_layer2[l2].strip);
0278 phi2 = phi2 << 2;
0279 if (vcluster_size[id] == 2 && itr1 == 1) {
0280 itr1 = 0;
0281 phi2 = phi2 + (radialAngle(vrpc_hit_layer2[l2 - 1].detid, rpcGeometry, vrpc_hit_layer2[l2 - 1].strip) << 2);
0282 phi2 /= 2;
0283 }
0284 l1ttma_hits_out.emplace_back(
0285 vrpc_hit_layer2[l2].bx, wh, sec - 1, st, phi2, 0, 3, hit[vrpc_hit_layer2[l2].bx], 0, 2);
0286 hit[vrpc_hit_layer2[l2].bx]++;
0287 }
0288 itr1 = 0;
0289
0290 for (unsigned int l1 = 0; l1 < vrpc_hit_st3.size(); l1++) {
0291 if (st != 3 || vrpc_hit_st3[l1].station != 3 || vrpc_hit_st3[l1].wheel != wh ||
0292 vrpc_hit_st3[l1].sector != sec)
0293 continue;
0294 RPCHitCleaner::detId_Ext tmp{vrpc_hit_st3[l1].detid, vrpc_hit_st3[l1].bx, vrpc_hit_st3[l1].strip};
0295 int id = hits[tmp];
0296
0297 if (vcluster_size[id] == 2 && itr1 == 0) {
0298 itr1++;
0299 continue;
0300 }
0301 int phi2 = radialAngle(vrpc_hit_st3[l1].detid, rpcGeometry, vrpc_hit_st3[l1].strip);
0302 phi2 = phi2 << 2;
0303 if (vcluster_size[id] == 2 && itr1 == 1) {
0304 itr1 = 0;
0305 phi2 = phi2 + (radialAngle(vrpc_hit_st3[l1 - 1].detid, rpcGeometry, vrpc_hit_st3[l1 - 1].strip) << 2);
0306 phi2 /= 2;
0307 }
0308 l1ttma_hits_out.emplace_back(
0309 vrpc_hit_st3[l1].bx, wh, sec - 1, st, phi2, 0, 3, hit[vrpc_hit_st3[l1].bx], 0, 2);
0310 hit[vrpc_hit_st3[l1].bx]++;
0311 }
0312 itr1 = 0;
0313
0314 for (unsigned int l1 = 0; l1 < vrpc_hit_st4.size(); l1++) {
0315 if (st != 4 || vrpc_hit_st4[l1].station != 4 || vrpc_hit_st4[l1].wheel != wh ||
0316 vrpc_hit_st4[l1].sector != sec)
0317 continue;
0318
0319 RPCHitCleaner::detId_Ext tmp{vrpc_hit_st4[l1].detid, vrpc_hit_st4[l1].bx, vrpc_hit_st4[l1].strip};
0320 int id = hits[tmp];
0321 if (vcluster_size[id] == 2 && itr1 == 0) {
0322 itr1++;
0323 continue;
0324 }
0325 int phi2 = radialAngle(vrpc_hit_st4[l1].detid, rpcGeometry, vrpc_hit_st4[l1].strip);
0326 phi2 = phi2 << 2;
0327 if (vcluster_size[id] == 2 && itr1 == 1) {
0328 itr1 = 0;
0329 phi2 = phi2 + (radialAngle(vrpc_hit_st4[l1 - 1].detid, rpcGeometry, vrpc_hit_st4[l1 - 1].strip) << 2);
0330 phi2 /= 2;
0331 }
0332 l1ttma_hits_out.emplace_back(
0333 vrpc_hit_st4[l1].bx, wh, sec - 1, st, phi2, 0, 3, hit[vrpc_hit_st4[l1].bx], 0, 2);
0334 hit[vrpc_hit_st4[l1].bx]++;
0335
0336
0337
0338 }
0339 }
0340 }
0341 }
0342
0343 m_rpcdt_translated.setContainer(l1ttma_out);
0344
0345 m_rpchitsdt_translated.setContainer(l1ttma_hits_out);
0346 }
0347
0348
0349 int RPCtoDTTranslator::radialAngle(RPCDetId detid, const RPCGeometry& rpcGeometry, int strip) {
0350 int radialAngle;
0351
0352
0353 const RPCRoll* roll = rpcGeometry.roll(detid);
0354 GlobalPoint stripPosition = roll->toGlobal(roll->centreOfStrip(strip));
0355
0356 double globalphi = stripPosition.phi();
0357 int sector = (roll->id()).sector();
0358 if (sector == 1)
0359 radialAngle = int(globalphi * 1024);
0360 else {
0361 if (globalphi >= 0)
0362 radialAngle = int((globalphi - (sector - 1) * Geom::pi() / 6.) * 1024);
0363 else
0364 radialAngle = int((globalphi + (13 - sector) * Geom::pi() / 6.) * 1024);
0365 }
0366 return radialAngle;
0367 }
0368
0369
0370 int RPCtoDTTranslator::localX(RPCDetId detid, const RPCGeometry& rpcGeometry, int strip) {
0371 const RPCRoll* roll = rpcGeometry.roll(detid);
0372
0373
0374 GlobalPoint p1cmPRG = roll->toGlobal(LocalPoint(1, 0, 0));
0375 GlobalPoint m1cmPRG = roll->toGlobal(LocalPoint(-1, 0, 0));
0376 float phip1cm = p1cmPRG.phi();
0377 float phim1cm = m1cmPRG.phi();
0378 int direction = (phip1cm - phim1cm) / abs(phip1cm - phim1cm);
0379
0380
0381 return direction * roll->centreOfStrip(strip).x();
0382 }
0383
0384 int RPCtoDTTranslator::bendingAngle(int xin, int xout, int phi) {
0385
0386 int atanv = (int)(atan((xout - xin) / 34.6) * 512);
0387 if (atanv > 512)
0388 return 512;
0389 int rvalue = atanv - phi / 8;
0390 return rvalue;
0391 }
0392
0393 int RPCtoDTTranslator::localXX(int phi, int layer, int station) {
0394
0395 double R[2][2] = {{410.0, 444.8}, {492.7, 527.3}};
0396 double rvalue = R[station - 1][layer - 1] * tan(phi / 4096.);
0397 return rvalue;
0398 }