File indexing completed on 2022-02-02 23:43:35
0001 #include "L1Trigger/L1TMuonOverlapPhase1/interface/AngleConverterBase.h"
0002 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Omtf/OMTFConfiguration.h"
0003
0004 #include "DataFormats/CSCDigi/interface/CSCConstants.h"
0005 #include "L1Trigger/DTUtilities/interface/DTTrigGeom.h"
0006
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0009 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0010 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0011 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0012 #include "DataFormats/CSCDigi/interface/CSCCorrelatedLCTDigi.h"
0013 #include "DataFormats/L1DTTrackFinder/interface/L1MuDTChambPhDigi.h"
0014 #include "DataFormats/L1DTTrackFinder/interface/L1MuDTChambThContainer.h"
0015 #include "DataFormats/RPCDigi/interface/RPCDigi.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 #include <cmath>
0019
0020 namespace {
0021 template <typename T>
0022 int sgn(T val) {
0023 return (T(0) < val) - (val < T(0));
0024 }
0025
0026 int fixCscOffsetGeom(int offsetLoc) {
0027
0028
0029
0030 const std::vector<int> offCSC = {-154, -133, -17, -4, 4, 17, 133, 146, 154, 167, 283, 296, 304, 317,
0031 433, 446, 454, 467, 583, 596, 604, 617, 733, 746, 754, 767, 883, 904};
0032 auto gep = std::lower_bound(offCSC.begin(), offCSC.end(), offsetLoc);
0033 int fixOff = (gep != offCSC.end()) ? *gep : *(gep - 1);
0034 if (gep != offCSC.begin() && std::abs(*(gep - 1) - offsetLoc) < std::abs(fixOff - offsetLoc))
0035 fixOff = *(gep - 1);
0036 return fixOff;
0037 }
0038
0039 }
0040
0041 AngleConverterBase::AngleConverterBase() : _geom_cache_id(0ULL) {}
0042
0043
0044 AngleConverterBase::~AngleConverterBase() {}
0045
0046
0047 void AngleConverterBase::checkAndUpdateGeometry(const edm::EventSetup& es,
0048 const ProcConfigurationBase* config,
0049 const MuonGeometryTokens& muonGeometryTokens) {
0050 if (muonGeometryRecordWatcher.check(es)) {
0051 _georpc = es.getHandle(muonGeometryTokens.rpcGeometryEsToken);
0052 _geocsc = es.getHandle(muonGeometryTokens.cscGeometryEsToken);
0053 _geodt = es.getHandle(muonGeometryTokens.dtGeometryEsToken);
0054 }
0055 this->config = config;
0056 nPhiBins = config->nPhiBins();
0057 }
0058
0059
0060 int AngleConverterBase::getProcessorPhi(int phiZero, l1t::tftype part, int dtScNum, int dtPhi) const {
0061 int dtPhiBins = 4096;
0062
0063 double hsPhiPitch = 2 * M_PI / nPhiBins;
0064
0065 int sector = dtScNum + 1;
0066
0067 double scale = 1. / dtPhiBins / hsPhiPitch;
0068 int scale_coeff = lround(scale * pow(2, 11));
0069
0070 int ichamber = sector - 1;
0071 if (ichamber > 6)
0072 ichamber = ichamber - 12;
0073
0074 int offsetGlobal = (int)nPhiBins * ichamber / 12;
0075
0076 int phiConverted = floor(dtPhi * scale_coeff / pow(2, 11)) + offsetGlobal - phiZero;
0077
0078
0079 return config->foldPhi(phiConverted);
0080 }
0081
0082
0083 int AngleConverterBase::getProcessorPhi(int phiZero,
0084 l1t::tftype part,
0085 const CSCDetId& csc,
0086 const CSCCorrelatedLCTDigi& digi) const {
0087 const double hsPhiPitch = 2 * M_PI / nPhiBins;
0088
0089
0090
0091
0092
0093 int halfStrip = digi.getStrip();
0094
0095 const CSCChamber* chamber = _geocsc->chamber(csc);
0096
0097
0098
0099 if (csc.station() == 1 && csc.ring() == 1 && halfStrip > 128) {
0100 CSCDetId cscME11 = CSCDetId(csc.endcap(), csc.station(), 4, csc.chamber());
0101 chamber = _geocsc->chamber(cscME11);
0102 }
0103
0104 const CSCChamberSpecs* cspec = chamber->specs();
0105 const CSCLayer* layer = chamber->layer(3);
0106 int order = (layer->centerOfStrip(2).phi() - layer->centerOfStrip(1).phi() > 0) ? 1 : -1;
0107 double stripPhiPitch = cspec->stripPhiPitch();
0108 double scale = fabs(stripPhiPitch / hsPhiPitch / 2.);
0109 if (fabs(scale - 1.) < 0.0002)
0110 scale = 1.;
0111
0112 double phiHalfStrip0 = layer->centerOfStrip(1).phi() - order * stripPhiPitch / 4.;
0113
0114 int offsetLoc = lround((phiHalfStrip0) / hsPhiPitch - phiZero);
0115 offsetLoc = config->foldPhi(offsetLoc);
0116
0117 if (csc.station() == 1 && csc.ring() == 1 && halfStrip > 128) {
0118
0119
0120
0121
0122 halfStrip -= 128;
0123 }
0124
0125
0126
0127 int fixOff = offsetLoc;
0128
0129
0130
0131 if (config->getFixCscGeometryOffset())
0132 fixOff = fixCscOffsetGeom(offsetLoc);
0133
0134 int phi = fixOff + order * scale * halfStrip;
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155 return config->foldPhi(phi);
0156 }
0157
0158
0159
0160 int AngleConverterBase::getProcessorPhi(
0161 int phiZero, l1t::tftype part, const RPCDetId& rollId, const unsigned int& digi1, const unsigned int& digi2) const {
0162 const double hsPhiPitch = 2 * M_PI / nPhiBins;
0163 const int dummy = nPhiBins;
0164 const RPCRoll* roll = _georpc->roll(rollId);
0165 if (!roll)
0166 return dummy;
0167
0168 double stripPhi1 = (roll->toGlobal(roll->centreOfStrip((int)digi1))).phi();
0169 double stripPhi2 = (roll->toGlobal(roll->centreOfStrip((int)digi2))).phi();
0170
0171
0172 if (std::signbit(stripPhi1) != std::signbit(stripPhi2) && abs(stripPhi1) > M_PI / 2.) {
0173 if (std::signbit(stripPhi1)) {
0174 stripPhi1 += 2 * M_PI;
0175 } else
0176 stripPhi2 += 2 * M_PI;
0177 }
0178 int halfStrip = lround(((stripPhi1 + stripPhi2) / 2.) / hsPhiPitch);
0179 halfStrip = config->foldPhi(halfStrip);
0180
0181 LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << 185 << " roll " << rollId.rawId() << " " << rollId
0182 << " cluster: firstStrip " << digi1 << " stripPhi1Global " << stripPhi1
0183 << " stripPhi1LocalPhi " << roll->centreOfStrip((int)digi1).x() << " y "
0184 << roll->centreOfStrip((int)digi1).y() << " lastStrip " << digi2 << " stripPhi2Global "
0185 << stripPhi2 << " stripPhi2LocalPhi x " << roll->centreOfStrip((int)digi2).x() << " y "
0186 << roll->centreOfStrip((int)digi2).y() << " halfStrip " << halfStrip << std::endl;
0187
0188 return config->foldPhi(halfStrip - phiZero);
0189 }
0190
0191 int AngleConverterBase::getProcessorPhi(unsigned int iProcessor,
0192 l1t::tftype part,
0193 const RPCDetId& rollId,
0194 const unsigned int& digi) const {
0195 const double hsPhiPitch = 2 * M_PI / nPhiBins;
0196 const int dummy = nPhiBins;
0197 int processor = iProcessor + 1;
0198 const RPCRoll* roll = _georpc->roll(rollId);
0199 if (!roll)
0200 return dummy;
0201
0202 double phi15deg = M_PI / 3. * (processor - 1) + M_PI / 12.;
0203
0204
0205 double stripPhi = (roll->toGlobal(roll->centreOfStrip((int)digi))).phi();
0206
0207
0208 switch (processor) {
0209 case 1:
0210 break;
0211 case 6: {
0212 phi15deg -= 2 * M_PI;
0213 break;
0214 }
0215 default: {
0216 if (stripPhi < 0)
0217 stripPhi += 2 * M_PI;
0218 break;
0219 }
0220 }
0221
0222
0223 int halfStrip = lround((stripPhi - phi15deg) / hsPhiPitch);
0224
0225 return halfStrip;
0226 }
0227
0228
0229 EtaValue AngleConverterBase::getGlobalEtaDt(const DTChamberId& detId) const {
0230 Local2DPoint chamberMiddleLP(0, 0);
0231 GlobalPoint chamberMiddleGP = _geodt->chamber(detId)->toGlobal(chamberMiddleLP);
0232
0233 const DTChamberId baseidNeigh(detId.wheel() + (detId.wheel() >= 0 ? -1 : +1), detId.station(), detId.sector());
0234 GlobalPoint chambNeighMiddleGP = _geodt->chamber(baseidNeigh)->toGlobal(chamberMiddleLP);
0235
0236 EtaValue etaValue = {
0237 config->etaToHwEta(chamberMiddleGP.eta()),
0238 config->etaToHwEta(fabs(chamberMiddleGP.eta() - chambNeighMiddleGP.eta())) / 2,
0239 0,
0240 0,
0241 0
0242 };
0243
0244
0245 return etaValue;
0246 }
0247
0248
0249
0250 void AngleConverterBase::getGlobalEta(const L1MuDTChambThDigi& thetaDigi, std::vector<EtaValue>& etaSegments) const {
0251 const DTChamberId baseid(thetaDigi.whNum(), thetaDigi.stNum(), thetaDigi.scNum() + 1);
0252 DTTrigGeom trig_geom(_geodt->chamber(baseid), false);
0253
0254
0255
0256
0257
0258
0259
0260 const int NBTI_theta = trig_geom.nCell(2);
0261 for (unsigned int btiGroup = 0; btiGroup < 7; ++btiGroup) {
0262 if (thetaDigi.position(btiGroup)) {
0263 unsigned btiActual = btiGroup * NBTI_theta / 7 + NBTI_theta / 14 + 1;
0264 DTBtiId thetaBTI = DTBtiId(baseid, 2, btiActual);
0265 GlobalPoint theta_gp = trig_geom.CMSPosition(thetaBTI);
0266
0267 EtaValue etaValue = {
0268 config->etaToHwEta(theta_gp.eta()),
0269 0,
0270 thetaDigi.quality(btiGroup),
0271 thetaDigi.bxNum(),
0272 0
0273 };
0274 etaSegments.emplace_back(etaValue);
0275
0276
0277 }
0278 }
0279 }
0280
0281 std::vector<EtaValue> AngleConverterBase::getGlobalEta(const L1MuDTChambThContainer* dtThDigis,
0282 int bxFrom,
0283 int bxTo) const {
0284
0285
0286 std::vector<EtaValue> etaSegments;
0287
0288 for (auto& thetaDigi : (*(dtThDigis->getContainer()))) {
0289 if (thetaDigi.bxNum() >= bxFrom && thetaDigi.bxNum() <= bxTo) {
0290 getGlobalEta(thetaDigi, etaSegments);
0291 }
0292 }
0293 return etaSegments;
0294 }
0295
0296
0297 float AngleConverterBase::cscChamberEtaSize(const CSCDetId& detId) const {
0298 if (detId.station() == 1) {
0299 if (detId.ring() == 1)
0300 return (2.5 - 1.6) / 2.;
0301
0302 if (detId.ring() == 2)
0303 return (1.7 - 1.2) / 2.;
0304 if (detId.ring() == 3)
0305 return (1.12 - 0.9) / 2.;
0306 if (detId.ring() == 4)
0307 return (2.5 - 1.6) / 2.;
0308 } else if (detId.station() == 2) {
0309 if (detId.ring() == 1)
0310 return (2.5 - 1.6) / 2.;
0311 if (detId.ring() == 2)
0312 return (1.6 - 1.0) / 2.;
0313 } else if (detId.station() == 3) {
0314 if (detId.ring() == 1)
0315 return (2.5 - 1.7) / 2.;
0316 if (detId.ring() == 2)
0317 return (1.7 - 1.1) / 2.;
0318 } else if (detId.station() == 4) {
0319 if (detId.ring() == 1)
0320 return (2.45 - 1.8) / 2.;
0321 if (detId.ring() == 2)
0322 return (1.8 - 1.2) / 2.;
0323 }
0324 return 0;
0325 }
0326
0327 EtaValue AngleConverterBase::getGlobalEta(const CSCDetId& detId, const CSCCorrelatedLCTDigi& aDigi) const {
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337 auto chamb = _geocsc->chamber(detId);
0338 auto layer_geom = chamb->layer(CSCConstants::KEY_ALCT_LAYER)->geometry();
0339 auto layer = chamb->layer(CSCConstants::KEY_ALCT_LAYER);
0340
0341 const uint16_t keyWG = aDigi.getKeyWG();
0342
0343 const LocalPoint lpWg = layer_geom->localCenterOfWireGroup(keyWG);
0344 const GlobalPoint gpWg = layer->surface().toGlobal(lpWg);
0345
0346 EtaValue etaSegment = {
0347 config->etaToHwEta(gpWg.eta()),
0348 0,
0349 0,
0350 aDigi.getBX(),
0351 0
0352 };
0353
0354
0355 return etaSegment;
0356 }
0357
0358
0359
0360 EtaValue AngleConverterBase::getGlobalEtaCsc(const CSCDetId& detId) const {
0361 auto chamb = _geocsc->chamber(detId);
0362
0363 Local2DPoint chamberMiddleLP(0, 0);
0364 GlobalPoint chamberMiddleGP = chamb->toGlobal(chamberMiddleLP);
0365
0366 EtaValue etaValue = {
0367 config->etaToHwEta(chamberMiddleGP.eta()),
0368 config->etaToHwEta(cscChamberEtaSize(detId)),
0369 0,
0370 0,
0371 0
0372 };
0373
0374
0375 return etaValue;
0376 }
0377
0378
0379
0380 EtaValue AngleConverterBase::getGlobalEta(unsigned int rawid, const unsigned int& strip) const {
0381 const RPCDetId id(rawid);
0382
0383 auto roll = _georpc->roll(id);
0384 const LocalPoint lp = roll->centreOfStrip((int)strip);
0385 const GlobalPoint gp = roll->toGlobal(lp);
0386
0387 int neighbRoll = 1;
0388
0389 if (id.region() == 0) {
0390 if (id.station() == 2 &&
0391 ((abs(id.ring()) == 2 && id.layer() == 2) || (abs(id.ring()) != 2 && id.layer() == 1))) {
0392 if (id.roll() == 2)
0393 neighbRoll = 1;
0394 else {
0395 neighbRoll = 2;
0396 }
0397 } else
0398 neighbRoll = (id.roll() == 1 ? 3 : 1);
0399 } else {
0400 neighbRoll = id.roll() + (id.roll() == 1 ? +1 : -1);
0401 }
0402
0403 const RPCDetId idNeigh =
0404 RPCDetId(id.region(), id.ring(), id.station(), id.sector(), id.layer(), id.subsector(), neighbRoll);
0405
0406
0407
0408 auto rollNeigh = _georpc->roll(idNeigh);
0409 const LocalPoint lpNeigh = rollNeigh->centreOfStrip((int)strip);
0410 const GlobalPoint gpNeigh = rollNeigh->toGlobal(lpNeigh);
0411
0412 EtaValue etaValue = {config->etaToHwEta(gp.eta()),
0413 config->etaToHwEta(abs(gp.eta() - gpNeigh.eta())) /
0414 2,
0415 0};
0416
0417
0418 return etaValue;
0419 }
0420
0421
0422 bool AngleConverterBase::isCSCCounterClockwise(const CSCLayer* layer) const {
0423 const int nStrips = layer->geometry()->numberOfStrips();
0424 const double phi1 = layer->centerOfStrip(1).phi();
0425 const double phiN = layer->centerOfStrip(nStrips).phi();
0426 return ((std::abs(phi1 - phiN) < M_PI && phi1 >= phiN) || (std::abs(phi1 - phiN) >= M_PI && phi1 < phiN));
0427 }
0428
0429
0430 const int AngleConverterBase::findBTIgroup(const L1MuDTChambPhDigi& aDigi, const L1MuDTChambThContainer* dtThDigis) {
0431 int bti_group = -1;
0432
0433 const L1MuDTChambThDigi* theta_segm =
0434 dtThDigis->chThetaSegm(aDigi.whNum(), aDigi.stNum(), aDigi.scNum(), aDigi.bxNum());
0435 if (!theta_segm)
0436 return bti_group;
0437
0438 for (unsigned int i = 0; i < 7; ++i) {
0439 if (theta_segm->position(i) && bti_group < 0)
0440 bti_group = i;
0441
0442
0443
0444 else if (theta_segm->position(i) && bti_group > -1)
0445 return -1;
0446 }
0447
0448 return bti_group;
0449 }
0450
0451