File indexing completed on 2024-04-06 12:21:04
0001 #include "L1Trigger/L1TMuonOverlap/interface/AngleConverter.h"
0002 #include "L1Trigger/L1TMuonOverlap/interface/OMTFConfiguration.h"
0003
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Utilities/interface/Transition.h"
0007
0008 #include "L1Trigger/CSCTriggerPrimitives/interface/CSCPatternBank.h"
0009
0010 #include "L1Trigger/DTUtilities/interface/DTTrigGeom.h"
0011
0012 #include "DataFormats/CSCDigi/interface/CSCConstants.h"
0013 #include "DataFormats/CSCDigi/interface/CSCCorrelatedLCTDigi.h"
0014 #include "DataFormats/L1DTTrackFinder/interface/L1MuDTChambPhDigi.h"
0015 #include "DataFormats/L1DTTrackFinder/interface/L1MuDTChambThContainer.h"
0016 #include "DataFormats/RPCDigi/interface/RPCDigi.h"
0017
0018 #include <cmath>
0019 #include <array>
0020
0021 namespace {
0022 template <typename T>
0023 int sgn(T val) {
0024 return (T(0) < val) - (val < T(0));
0025 }
0026
0027 constexpr std::array<float, 8> bounds = {{1.24, 1.14353, 1.09844, 1.05168, 1.00313, 0.952728, 0.90037, 0.8}};
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040 int etaVal2Bit(float eta) { return bounds.rend() - std::lower_bound(bounds.rbegin(), bounds.rend(), fabs(eta)); }
0041
0042 int etaBit2Code(unsigned int bit) {
0043 int code = 73;
0044 switch (bit) {
0045 case 0: {
0046 code = 115;
0047 break;
0048 }
0049 case 1: {
0050 code = 110;
0051 break;
0052 }
0053 case 2: {
0054 code = 103;
0055 break;
0056 }
0057 case 3: {
0058 code = 99;
0059 break;
0060 }
0061 case 4: {
0062 code = 94;
0063 break;
0064 }
0065 case 5: {
0066 code = 90;
0067 break;
0068 }
0069 case 6: {
0070 code = 85;
0071 break;
0072 }
0073 case 7: {
0074 code = 78;
0075 break;
0076 }
0077 case 8: {
0078 code = 73;
0079 break;
0080 }
0081 default: {
0082 code = 95;
0083 break;
0084 }
0085 }
0086 return code;
0087 }
0088
0089 int etaVal2Code(double etaVal) {
0090 int sign = sgn(etaVal);
0091 int bit = etaVal2Bit(fabs(etaVal));
0092 int code = etaBit2Code(bit);
0093 return sign * code;
0094 }
0095
0096 int etaKeyWG2Code(const CSCDetId &detId, uint16_t keyWG) {
0097 signed int etaCode = 121;
0098 if (detId.station() == 1 && detId.ring() == 2) {
0099 if (keyWG < 49)
0100 etaCode = 121;
0101 else if (keyWG <= 57)
0102 etaCode = etaBit2Code(0);
0103 else if (keyWG <= 63)
0104 etaCode = etaBit2Code(1);
0105 } else if (detId.station() == 1 && detId.ring() == 3) {
0106 if (keyWG <= 2)
0107 etaCode = etaBit2Code(2);
0108 else if (keyWG <= 8)
0109 etaCode = etaBit2Code(3);
0110 else if (keyWG <= 15)
0111 etaCode = etaBit2Code(4);
0112 else if (keyWG <= 23)
0113 etaCode = etaBit2Code(5);
0114 else if (keyWG <= 31)
0115 etaCode = etaBit2Code(6);
0116 } else if ((detId.station() == 2 || detId.station() == 3) && detId.ring() == 2) {
0117 if (keyWG < 24)
0118 etaCode = 121;
0119 else if (keyWG <= 29)
0120 etaCode = etaBit2Code(0);
0121 else if (keyWG <= 43)
0122 etaCode = etaBit2Code(1);
0123 else if (keyWG <= 49)
0124 etaCode = etaBit2Code(2);
0125 else if (keyWG <= 56)
0126 etaCode = etaBit2Code(3);
0127 else if (keyWG <= 63)
0128 etaCode = etaBit2Code(4);
0129 }
0130
0131 if (detId.endcap() == 2)
0132 etaCode *= -1;
0133 return etaCode;
0134 }
0135
0136 int fixCscOffsetGeom(int offsetLoc) {
0137
0138
0139
0140 const std::vector<int> offCSC = {-154, -133, -17, -4, 4, 17, 133, 146, 154, 167, 283, 296, 304, 317,
0141 433, 446, 454, 467, 583, 596, 604, 617, 733, 746, 754, 767, 883, 904};
0142 auto gep = std::lower_bound(offCSC.begin(), offCSC.end(), offsetLoc);
0143 int fixOff = (gep != offCSC.end()) ? *gep : *(gep - 1);
0144 if (gep != offCSC.begin() && std::abs(*(gep - 1) - offsetLoc) < std::abs(fixOff - offsetLoc))
0145 fixOff = *(gep - 1);
0146 return fixOff;
0147 }
0148
0149 }
0150
0151 AngleConverter::AngleConverter(edm::ConsumesCollector &iC, bool getDuringEvent) : _geom_cache_id(0ULL) {
0152 if (getDuringEvent) {
0153 rpcGeometryToken_ = iC.esConsumes<RPCGeometry, MuonGeometryRecord>();
0154 cscGeometryToken_ = iC.esConsumes<CSCGeometry, MuonGeometryRecord>();
0155 dtGeometryToken_ = iC.esConsumes<DTGeometry, MuonGeometryRecord>();
0156 } else {
0157 rpcGeometryToken_ = iC.esConsumes<RPCGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0158 cscGeometryToken_ = iC.esConsumes<CSCGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0159 dtGeometryToken_ = iC.esConsumes<DTGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0160 }
0161 }
0162
0163
0164 AngleConverter::~AngleConverter() {}
0165
0166
0167 void AngleConverter::checkAndUpdateGeometry(const edm::EventSetup &es, unsigned int phiBins) {
0168 const MuonGeometryRecord &geom = es.get<MuonGeometryRecord>();
0169 unsigned long long geomid = geom.cacheIdentifier();
0170 if (_geom_cache_id != geomid) {
0171 _georpc = &geom.get(rpcGeometryToken_);
0172 _geocsc = &geom.get(cscGeometryToken_);
0173 _geodt = &geom.get(dtGeometryToken_);
0174 _geom_cache_id = geomid;
0175 }
0176
0177 nPhiBins = phiBins;
0178 }
0179
0180
0181 int AngleConverter::getProcessorPhi(unsigned int iProcessor, l1t::tftype part, const L1MuDTChambPhDigi &digi) const {
0182 double hsPhiPitch = 2 * M_PI / nPhiBins;
0183 const int dummy = nPhiBins;
0184 int processor = iProcessor + 1;
0185 int posneg = (part == l1t::tftype::omtf_pos) ? 1 : -1;
0186
0187 int sector =
0188 digi.scNum() + 1;
0189 int wheel = digi.whNum();
0190 int station = digi.stNum();
0191 int phiDT = digi.phi();
0192
0193 if (posneg * 2 != wheel)
0194 return dummy;
0195 if (station > 3)
0196 return dummy;
0197
0198
0199 if (processor != 6 && !(sector >= processor * 2 - 1 && sector <= processor * 2 + 1))
0200 return dummy;
0201 if (processor == 6 && !(sector >= 11 || sector == 1))
0202 return dummy;
0203
0204
0205 int ichamber = sector - 1 - 2 * (processor - 1);
0206 if (ichamber < 0)
0207 ichamber += 12;
0208
0209 int offsetLoc = lround(((ichamber - 1) * M_PI / 6 + M_PI / 12.) / hsPhiPitch);
0210 double scale = 1. / 4096 / hsPhiPitch;
0211 int scale_coeff = lround(scale * pow(2, 11));
0212
0213 int phi = floor(phiDT * scale_coeff / pow(2, 11)) + offsetLoc;
0214
0215 return phi;
0216 }
0217
0218
0219 int AngleConverter::getProcessorPhi(unsigned int iProcessor,
0220 l1t::tftype part,
0221 const CSCDetId &csc,
0222 const CSCCorrelatedLCTDigi &digi) const {
0223 const double hsPhiPitch = 2 * M_PI / nPhiBins;
0224 const int dummy = nPhiBins;
0225 int processor = iProcessor + 1;
0226 int posneg = (part == l1t::tftype::omtf_pos) ? 1 : -1;
0227
0228
0229
0230 if (posneg != csc.zendcap())
0231 return dummy;
0232 if (csc.ring() != 3 && !(csc.ring() == 2 && (csc.station() == 2 || csc.station() == 3 || csc.station() == 1)))
0233 return dummy;
0234 if (processor != 6) {
0235 if (csc.chamber() < (processor - 1) * 6 + 2)
0236 return dummy;
0237 if (csc.chamber() > (processor - 1) * 6 + 8)
0238 return dummy;
0239 } else {
0240 if (csc.chamber() > 2 && csc.chamber() < 32)
0241 return dummy;
0242 }
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254 const CSCChamber *chamber = _geocsc->chamber(csc);
0255 const CSCChamberSpecs *cspec = chamber->specs();
0256 const CSCLayer *layer = chamber->layer(3);
0257 int order = (layer->centerOfStrip(2).phi() - layer->centerOfStrip(1).phi() > 0) ? 1 : -1;
0258 double stripPhiPitch = cspec->stripPhiPitch();
0259 double scale = fabs(stripPhiPitch / hsPhiPitch / 2.);
0260 if (fabs(scale - 1.) < 0.0002)
0261 scale = 1.;
0262 double phi15deg = M_PI / 3. * (processor - 1) + M_PI / 12.;
0263 double phiHalfStrip0 = layer->centerOfStrip(10).phi() - order * 9 * stripPhiPitch - order * stripPhiPitch / 4.;
0264 if (processor == 6 || phiHalfStrip0 < 0)
0265 phiHalfStrip0 += 2 * M_PI;
0266 int offsetLoc = lround((phiHalfStrip0 - phi15deg) / hsPhiPitch);
0267
0268 int halfStrip = digi.getStrip();
0269
0270
0271
0272
0273 int fixOff = fixCscOffsetGeom(offsetLoc);
0274
0275 int phi = fixOff + order * scale * halfStrip;
0276
0277
0278
0279
0280 return phi;
0281 }
0282
0283
0284
0285 int AngleConverter::getProcessorPhi(unsigned int iProcessor,
0286 l1t::tftype part,
0287 const RPCDetId &rollId,
0288 const unsigned int &digi1,
0289 const unsigned int &digi2) const {
0290 const double hsPhiPitch = 2 * M_PI / nPhiBins;
0291 const int dummy = nPhiBins;
0292 int processor = iProcessor + 1;
0293 const RPCRoll *roll = _georpc->roll(rollId);
0294 if (!roll)
0295 return dummy;
0296
0297 double phi15deg =
0298 M_PI / 3. * (processor - 1) + M_PI / 12.;
0299 double stripPhi1 = (roll->toGlobal(roll->centreOfStrip((int)digi1))).phi();
0300 double stripPhi2 = (roll->toGlobal(roll->centreOfStrip((int)digi2))).phi();
0301
0302
0303
0304
0305 switch (processor) {
0306 case 1:
0307 break;
0308 case 6: {
0309 phi15deg -= 2 * M_PI;
0310 break;
0311 }
0312 default: {
0313 if (stripPhi1 < 0)
0314 stripPhi1 += 2 * M_PI;
0315 if (stripPhi2 < 0)
0316 stripPhi2 += 2 * M_PI;
0317 break;
0318 }
0319 }
0320
0321
0322 int halfStrip = lround(((stripPhi1 + stripPhi2) / 2. - phi15deg) / hsPhiPitch);
0323
0324 return halfStrip;
0325 }
0326
0327 int AngleConverter::getProcessorPhi(unsigned int iProcessor,
0328 l1t::tftype part,
0329 const RPCDetId &rollId,
0330 const unsigned int &digi) const {
0331 const double hsPhiPitch = 2 * M_PI / nPhiBins;
0332 const int dummy = nPhiBins;
0333 int processor = iProcessor + 1;
0334 const RPCRoll *roll = _georpc->roll(rollId);
0335 if (!roll)
0336 return dummy;
0337
0338 double phi15deg =
0339 M_PI / 3. * (processor - 1) + M_PI / 12.;
0340 double stripPhi = (roll->toGlobal(roll->centreOfStrip((int)digi))).phi();
0341
0342
0343 switch (processor) {
0344 case 1:
0345 break;
0346 case 6: {
0347 phi15deg -= 2 * M_PI;
0348 break;
0349 }
0350 default: {
0351 if (stripPhi < 0)
0352 stripPhi += 2 * M_PI;
0353 break;
0354 }
0355 }
0356
0357
0358 int halfStrip = lround((stripPhi - phi15deg) / hsPhiPitch);
0359
0360 return halfStrip;
0361 }
0362
0363
0364 int AngleConverter::getGlobalEta(unsigned int rawid,
0365 const L1MuDTChambPhDigi &aDigi,
0366 const L1MuDTChambThContainer *dtThDigis) {
0367 const DTChamberId baseid(aDigi.whNum(), aDigi.stNum(), aDigi.scNum() + 1);
0368
0369
0370 std::unique_ptr<DTChamber> chamb(const_cast<DTChamber *>(_geodt->chamber(baseid)));
0371
0372 std::unique_ptr<DTTrigGeom> trig_geom(new DTTrigGeom(chamb.get(), false));
0373 chamb.release();
0374
0375
0376
0377
0378
0379
0380 const int NBTI_theta = ((baseid.station() != 4) ? trig_geom->nCell(2) : trig_geom->nCell(3));
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396 const L1MuDTChambThDigi *theta_segm =
0397 dtThDigis->chThetaSegm(aDigi.whNum(), aDigi.stNum(), aDigi.scNum(), aDigi.bxNum());
0398 int bti_group = -1;
0399 if (theta_segm) {
0400 for (unsigned int i = 0; i < 7; ++i)
0401 if (theta_segm->position(i) && bti_group < 0)
0402 bti_group = i;
0403 else if (theta_segm->position(i) && bti_group > -1)
0404 bti_group = 511;
0405 }
0406
0407 int iEta = 0;
0408 if (bti_group == 511)
0409 iEta = 95;
0410 else if (bti_group == -1 && aDigi.stNum() == 1)
0411 iEta = 92;
0412 else if (bti_group == -1 && aDigi.stNum() == 2)
0413 iEta = 79;
0414 else if (bti_group == -1 && aDigi.stNum() == 3)
0415 iEta = 75;
0416 else if (baseid.station() != 4 && bti_group >= 0) {
0417
0418 unsigned bti_actual = bti_group * NBTI_theta / 7 + NBTI_theta / 14 + 1;
0419 DTBtiId thetaBTI = DTBtiId(baseid, 2, bti_actual);
0420 GlobalPoint theta_gp = trig_geom->CMSPosition(thetaBTI);
0421 iEta = etaVal2Code(fabs(theta_gp.eta()));
0422 }
0423 int signEta = sgn(aDigi.whNum());
0424 iEta *= signEta;
0425 return iEta;
0426 }
0427
0428
0429
0430 int AngleConverter::getGlobalEta(unsigned int rawid, const CSCCorrelatedLCTDigi &aDigi) {
0431
0432
0433
0434
0435
0436
0437
0438
0439 const CSCDetId id(rawid);
0440
0441
0442 std::unique_ptr<const CSCChamber> chamb(_geocsc->chamber(id));
0443 std::unique_ptr<const CSCLayerGeometry> layer_geom(chamb->layer(CSCConstants::KEY_ALCT_LAYER)->geometry());
0444 std::unique_ptr<const CSCLayer> layer(chamb->layer(CSCConstants::KEY_ALCT_LAYER));
0445
0446 const uint16_t halfstrip = aDigi.getStrip();
0447 const uint16_t pattern = aDigi.getPattern();
0448 const uint16_t keyWG = aDigi.getKeyWG();
0449
0450
0451
0452
0453 double offset = 0.0;
0454 switch (1) {
0455 case 1:
0456 offset = CSCPatternBank::getLegacyPosition(pattern);
0457 }
0458 const unsigned halfstrip_offs = unsigned(0.5 + halfstrip + offset);
0459 const unsigned strip = halfstrip_offs / 2 + 1;
0460
0461
0462
0463 const LocalPoint coarse_lp = layer_geom->stripWireGroupIntersection(strip, keyWG);
0464 const GlobalPoint coarse_gp = layer->surface().toGlobal(coarse_lp);
0465
0466
0467
0468 const double hs_offset = layer_geom->stripPhiPitch() / 4.0;
0469
0470
0471 const bool ccw = isCSCCounterClockwise(layer);
0472
0473 const double phi_offset = ((halfstrip_offs % 2 ? 1 : -1) * (ccw ? -hs_offset : hs_offset));
0474
0475
0476
0477 const GlobalPoint final_gp(
0478 GlobalPoint::Polar(coarse_gp.theta(), (coarse_gp.phi().value() + phi_offset), coarse_gp.mag()));
0479
0480 chamb.release();
0481 layer_geom.release();
0482 layer.release();
0483
0484
0485
0486
0487
0488
0489 return etaKeyWG2Code(id, keyWG);
0490
0491
0492
0493
0494 }
0495
0496
0497 int AngleConverter::getGlobalEta(unsigned int rawid, const unsigned int &strip) {
0498 const RPCDetId id(rawid);
0499 std::unique_ptr<const RPCRoll> roll(_georpc->roll(id));
0500 const LocalPoint lp = roll->centreOfStrip((int)strip);
0501 const GlobalPoint gp = roll->toGlobal(lp);
0502 roll.release();
0503
0504 return etaVal2Code(gp.eta());
0505
0506
0507 }
0508
0509
0510 bool AngleConverter::isCSCCounterClockwise(const std::unique_ptr<const CSCLayer> &layer) const {
0511 const int nStrips = layer->geometry()->numberOfStrips();
0512 const double phi1 = layer->centerOfStrip(1).phi();
0513 const double phiN = layer->centerOfStrip(nStrips).phi();
0514 return ((std::abs(phi1 - phiN) < M_PI && phi1 >= phiN) || (std::abs(phi1 - phiN) >= M_PI && phi1 < phiN));
0515 }
0516
0517
0518 const int AngleConverter::findBTIgroup(const L1MuDTChambPhDigi &aDigi, const L1MuDTChambThContainer *dtThDigis) {
0519 int bti_group = -1;
0520
0521 const L1MuDTChambThDigi *theta_segm =
0522 dtThDigis->chThetaSegm(aDigi.whNum(), aDigi.stNum(), aDigi.scNum(), aDigi.bxNum());
0523 if (!theta_segm)
0524 return bti_group;
0525
0526 for (unsigned int i = 0; i < 7; ++i) {
0527 if (theta_segm->position(i) && bti_group < 0)
0528 bti_group = i;
0529
0530
0531
0532 else if (theta_segm->position(i) && bti_group > -1)
0533 return -1;
0534 }
0535
0536 return bti_group;
0537 }
0538
0539