L1TMuon/src/GeometryTranslator.cc

0001 #include "L1Trigger/L1TMuon/interface/GeometryTranslator.h"
0002
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005
0006 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0007 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0008 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0009 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0010 #include "Geometry/GEMGeometry/interface/GEMGeometry.h"
0011 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0012 #include "MagneticField/Engine/interface/MagneticField.h"
0013 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0014
0015 #include "L1Trigger/DTUtilities/interface/DTTrigGeom.h"
0016 #include "DataFormats/CSCDigi/interface/CSCConstants.h"
0017 #include "L1Trigger/CSCTriggerPrimitives/interface/CSCPatternBank.h"
0018
0019 #include "L1Trigger/L1TMuon/interface/MuonTriggerPrimitive.h"
0020
0021 #include <cmath>  // for pi
0022
0023 using namespace L1TMuon;
0024
0025 GeometryTranslator::GeometryTranslator(edm::ConsumesCollector iC)
0026     : _geom_cache_id(0ULL),
0027       geodtToken_(iC.esConsumes()),
0028       geocscToken_(iC.esConsumes()),
0029       georpcToken_(iC.esConsumes()),
0030       geogemToken_(iC.esConsumes()),
0031       geome0Token_(iC.esConsumes()),
0032       _magfield_cache_id(0ULL),
0033       magfieldToken_(iC.esConsumes()) {}
0034
0035 GeometryTranslator::~GeometryTranslator() {}
0036
0037 double GeometryTranslator::calculateGlobalEta(const TriggerPrimitive& tp) const {
0038   switch (tp.subsystem()) {
0039     case L1TMuon::kDT:
0040       return calcDTSpecificEta(tp);
0041       break;
0042     case L1TMuon::kCSC:
0043       return calcCSCSpecificEta(tp);
0044       break;
0045     case L1TMuon::kRPC:
0046       return calcRPCSpecificEta(tp);
0047       break;
0048     case L1TMuon::kGEM:
0049       return calcGEMSpecificEta(tp);
0050       break;
0051     case L1TMuon::kME0:
0052       return calcME0SpecificEta(tp);
0053       break;
0054     default:
0055       return std::nan("Invalid TP type!");
0056       break;
0057   }
0058 }
0059
0060 double GeometryTranslator::calculateGlobalPhi(const TriggerPrimitive& tp) const {
0061   switch (tp.subsystem()) {
0062     case L1TMuon::kDT:
0063       return calcDTSpecificPhi(tp);
0064       break;
0065     case L1TMuon::kCSC:
0066       return calcCSCSpecificPhi(tp);
0067       break;
0068     case L1TMuon::kRPC:
0069       return calcRPCSpecificPhi(tp);
0070       break;
0071     case L1TMuon::kGEM:
0072       return calcGEMSpecificPhi(tp);
0073       break;
0074     case L1TMuon::kME0:
0075       return calcME0SpecificPhi(tp);
0076       break;
0077     default:
0078       return std::nan("Invalid TP type!");
0079       break;
0080   }
0081 }
0082
0083 double GeometryTranslator::calculateBendAngle(const TriggerPrimitive& tp) const {
0084   switch (tp.subsystem()) {
0085     case L1TMuon::kDT:
0086       return calcDTSpecificBend(tp);
0087       break;
0088     case L1TMuon::kCSC:
0089       return calcCSCSpecificBend(tp);
0090       break;
0091     case L1TMuon::kRPC:
0092       return calcRPCSpecificBend(tp);
0093       break;
0094     case L1TMuon::kGEM:
0095       return calcGEMSpecificBend(tp);
0096       break;
0097     case L1TMuon::kME0:
0098       return calcME0SpecificBend(tp);
0099       break;
0100     default:
0101       return std::nan("Invalid TP type!");
0102       break;
0103   }
0104 }
0105
0106 GlobalPoint GeometryTranslator::getGlobalPoint(const TriggerPrimitive& tp) const {
0107   switch (tp.subsystem()) {
0108     case L1TMuon::kDT:
0109       return calcDTSpecificPoint(tp);
0110       break;
0111     case L1TMuon::kCSC:
0112       return getCSCSpecificPoint(tp);
0113       break;
0114     case L1TMuon::kRPC:
0115       return getRPCSpecificPoint(tp);
0116       break;
0117     case L1TMuon::kGEM:
0118       return getGEMSpecificPoint(tp);
0119       break;
0120     case L1TMuon::kME0:
0121       return getME0SpecificPoint(tp);
0122       break;
0123     default:
0124       GlobalPoint ret(
0125           GlobalPoint::Polar(std::nan("Invalid TP type!"), std::nan("Invalid TP type!"), std::nan("Invalid TP type!")));
0126       return ret;
0127       break;
0128   }
0129 }
0130
0131 void GeometryTranslator::checkAndUpdateGeometry(const edm::EventSetup& es) {
0132   const MuonGeometryRecord& geom = es.get<MuonGeometryRecord>();
0133   unsigned long long geomid = geom.cacheIdentifier();
0134   if (_geom_cache_id != geomid) {
0135     _geodt = geom.getHandle(geodtToken_);
0136     _geocsc = geom.getHandle(geocscToken_);
0137     _georpc = geom.getHandle(georpcToken_);
0138     _geogem = geom.getHandle(geogemToken_);
0139     _geome0 = geom.getHandle(geome0Token_);
0140     _geom_cache_id = geomid;
0141   }
0142
0143   const IdealMagneticFieldRecord& magfield = es.get<IdealMagneticFieldRecord>();
0144   unsigned long long magfieldid = magfield.cacheIdentifier();
0145   if (_magfield_cache_id != magfieldid) {
0146     _magfield = magfield.getHandle(magfieldToken_);
0147     _magfield_cache_id = magfieldid;
0148   }
0149 }
0150
0151 // _____________________________________________________________________________
0152 // ME0
0153 GlobalPoint GeometryTranslator::getME0SpecificPoint(const TriggerPrimitive& tp) const {
0154   if (tp.detId<DetId>().subdetId() == MuonSubdetId::GEM) {  // GE0
0155     const GEMDetId id(tp.detId<GEMDetId>());
0156     const GEMSuperChamber* chamber = _geogem->superChamber(id);
0157     const GEMChamber* keylayer = chamber->chamber(3);  // GEM key layer is layer 3
0158     int partition = tp.getME0Data().partition;         // 'partition' is in half-roll unit
0159     int iroll = (partition >> 1) + 1;
0160     const GEMEtaPartition* roll = keylayer->etaPartition(iroll);
0161
0162     // failed to get ME0 roll
0163     if (roll == nullptr) {
0164       throw cms::Exception("Invalid GE0 Roll") << "Failed to get GEM roll for GE0" << std::endl;
0165     }
0166
0167     // See L1Trigger/ME0Trigger/src/ME0TriggerPseudoBuilder.cc
0168     int phiposition = tp.getME0Data().phiposition;  // 'phiposition' is in half-strip unit
0169     int istrip = (phiposition >> 1);
0170     int phiposition2 = (phiposition & 0x1);
0171     float centreOfStrip = istrip + 0.25 + phiposition2 * 0.5;
0172     const LocalPoint& lp = roll->centreOfStrip(centreOfStrip);
0173     const GlobalPoint& gp = roll->surface().toGlobal(lp);
0174     return gp;
0175   } else {  // ME0
0176     const ME0DetId id(tp.detId<ME0DetId>());
0177     const ME0Chamber* chamber = _geome0->chamber(id);
0178     const ME0Layer* keylayer = chamber->layer(3);  // ME0 key layer is layer 3
0179     int partition = tp.getME0Data().partition;     // 'partition' is in half-roll unit
0180     int iroll = (partition >> 1) + 1;
0181     const ME0EtaPartition* roll = keylayer->etaPartition(iroll);
0182
0183     // failed to get ME0 roll
0184     if (roll == nullptr) {
0185       throw cms::Exception("Invalid ME0 Roll") << "Failed to get ME0 roll" << std::endl;
0186     }
0187
0188     // See L1Trigger/ME0Trigger/src/ME0TriggerPseudoBuilder.cc
0189     int phiposition = tp.getME0Data().phiposition;  // 'phiposition' is in half-strip unit
0190     int istrip = (phiposition >> 1);
0191     int phiposition2 = (phiposition & 0x1);
0192     float centreOfStrip = istrip + 0.25 + phiposition2 * 0.5;
0193     const LocalPoint& lp = roll->centreOfStrip(centreOfStrip);
0194     const GlobalPoint& gp = roll->surface().toGlobal(lp);
0195     return gp;
0196   }
0197 }
0198
0199 double GeometryTranslator::calcME0SpecificEta(const TriggerPrimitive& tp) const {
0200   return getME0SpecificPoint(tp).eta();
0201 }
0202
0203 double GeometryTranslator::calcME0SpecificPhi(const TriggerPrimitive& tp) const {
0204   return getME0SpecificPoint(tp).phi();
0205 }
0206
0207 double GeometryTranslator::calcME0SpecificBend(const TriggerPrimitive& tp) const {
0208   return tp.getME0Data().deltaphi * (tp.getME0Data().bend == 0 ? 1 : -1);
0209 }
0210
0211 // _____________________________________________________________________________
0212 // GEM
0213 GlobalPoint GeometryTranslator::getGEMSpecificPoint(const TriggerPrimitive& tp) const {
0214   const GEMDetId id(tp.detId<GEMDetId>());
0215   const GEMEtaPartition* roll = _geogem->etaPartition(id);
0216   assert(roll != nullptr);  // failed to get GEM roll
0217   //const uint16_t pad = tp.getGEMData().pad;
0218   // Use half-pad precision, + 0.5 at the end to get the center of the pad (pad starts from 0)
0219   const float pad = (0.5 * static_cast<float>(tp.getGEMData().pad_low + tp.getGEMData().pad_hi)) + 0.5f;
0220   const LocalPoint& lp = roll->centreOfPad(pad);
0221   const GlobalPoint& gp = roll->surface().toGlobal(lp);
0222   return gp;
0223 }
0224
0225 double GeometryTranslator::calcGEMSpecificEta(const TriggerPrimitive& tp) const {
0226   return getGEMSpecificPoint(tp).eta();
0227 }
0228
0229 double GeometryTranslator::calcGEMSpecificPhi(const TriggerPrimitive& tp) const {
0230   return getGEMSpecificPoint(tp).phi();
0231 }
0232
0233 double GeometryTranslator::calcGEMSpecificBend(const TriggerPrimitive& tp) const { return 0.0; }
0234
0235 // _____________________________________________________________________________
0236 // RPC
0237 GlobalPoint GeometryTranslator::getRPCSpecificPoint(const TriggerPrimitive& tp) const {
0238   const RPCDetId id(tp.detId<RPCDetId>());
0239   const RPCRoll* roll = _georpc->roll(id);
0240   assert(roll != nullptr);  // failed to get RPC roll
0241   //const int strip = static_cast<int>(tp.getRPCData().strip);
0242   // Use half-strip precision, - 0.5 at the end to get the center of the strip (strip starts from 1)
0243   const float strip = (0.5 * static_cast<float>(tp.getRPCData().strip_low + tp.getRPCData().strip_hi)) - 0.5f;
0244   const LocalPoint& lp = roll->centreOfStrip(strip);
0245   const GlobalPoint& gp = roll->surface().toGlobal(lp);
0246   return gp;
0247 }
0248
0249 double GeometryTranslator::calcRPCSpecificEta(const TriggerPrimitive& tp) const {
0250   return getRPCSpecificPoint(tp).eta();
0251 }
0252
0253 double GeometryTranslator::calcRPCSpecificPhi(const TriggerPrimitive& tp) const {
0254   return getRPCSpecificPoint(tp).phi();
0255 }
0256
0257 // this function actually does nothing since RPC
0258 // hits are point-like objects
0259 double GeometryTranslator::calcRPCSpecificBend(const TriggerPrimitive& tp) const { return 0.0; }
0260
0261 // _____________________________________________________________________________
0262 // CSC
0263 //
0264 // alot of this is transcription and consolidation of the CSC
0265 // global phi calculation code
0266 // this works directly with the geometry
0267 // rather than using the old phi luts
0268 GlobalPoint GeometryTranslator::getCSCSpecificPoint(const TriggerPrimitive& tp) const {
0269   const CSCDetId id(tp.detId<CSCDetId>());
0270   // we should change this to weak_ptrs at some point
0271   // requires introducing std::shared_ptrs to geometry
0272   std::unique_ptr<const CSCChamber> chamb(_geocsc->chamber(id));
0273   assert(chamb != nullptr);  // failed to get CSC chamber
0274   std::unique_ptr<const CSCLayerGeometry> layer_geom(chamb->layer(CSCConstants::KEY_ALCT_LAYER)->geometry());
0275   std::unique_ptr<const CSCLayer> layer(chamb->layer(CSCConstants::KEY_ALCT_LAYER));
0276
0277   const uint16_t halfstrip = tp.getCSCData().strip;
0278   const uint16_t pattern = tp.getCSCData().pattern;
0279   const uint16_t keyWG = tp.getCSCData().keywire;
0280   //const unsigned maxStrips = layer_geom->numberOfStrips();
0281
0282   // so we can extend this later
0283   // assume TMB2007 half-strips only as baseline
0284   double offset = 0.0;
0285   switch (1) {
0286     case 1:
0287       offset = CSCPatternBank::getLegacyPosition(pattern);
0288   }
0289   const unsigned halfstrip_offs = static_cast<unsigned>(0.5 + halfstrip + offset);
0290   const unsigned strip = halfstrip_offs / 2 + 1;  // geom starts from 1
0291
0292   // the rough location of the hit at the ALCT key layer
0293   // we will refine this using the half strip information
0294   const LocalPoint& coarse_lp = layer_geom->stripWireGroupIntersection(strip, keyWG);
0295   const GlobalPoint& coarse_gp = layer->surface().toGlobal(coarse_lp);
0296
0297   // the strip width/4.0 gives the offset of the half-strip
0298   // center with respect to the strip center
0299   const double hs_offset = layer_geom->stripPhiPitch() / 4.0;
0300
0301   // determine handedness of the chamber
0302   const bool ccw = isCSCCounterClockwise(layer);
0303   // we need to subtract the offset of even half strips and add the odd ones
0304   const double phi_offset = ((halfstrip_offs % 2 ? 1 : -1) * (ccw ? -hs_offset : hs_offset));
0305
0306   // the global eta calculation uses the middle of the strip
0307   // so no need to increment it
0308   const GlobalPoint final_gp(
0309       GlobalPoint::Polar(coarse_gp.theta(), (coarse_gp.phi().value() + phi_offset), coarse_gp.mag()));
0310
0311   // We need to add in some notion of the 'error' on trigger primitives
0312   // like the width of the wire group by the width of the strip
0313   // or something similar
0314
0315   // release ownership of the pointers
0316   chamb.release();
0317   layer_geom.release();
0318   layer.release();
0319
0320   return final_gp;
0321 }
0322
0323 double GeometryTranslator::calcCSCSpecificEta(const TriggerPrimitive& tp) const {
0324   return getCSCSpecificPoint(tp).eta();
0325 }
0326
0327 double GeometryTranslator::calcCSCSpecificPhi(const TriggerPrimitive& tp) const {
0328   return getCSCSpecificPoint(tp).phi();
0329 }
0330
0331 double GeometryTranslator::calcCSCSpecificBend(const TriggerPrimitive& tp) const { return tp.getCSCData().bend; }
0332
0333 bool GeometryTranslator::isCSCCounterClockwise(const std::unique_ptr<const CSCLayer>& layer) const {
0334   const int nStrips = layer->geometry()->numberOfStrips();
0335   const double phi1 = layer->centerOfStrip(1).phi();
0336   const double phiN = layer->centerOfStrip(nStrips).phi();
0337   return ((std::abs(phi1 - phiN) < M_PI && phi1 >= phiN) || (std::abs(phi1 - phiN) >= M_PI && phi1 < phiN));
0338 }
0339
0340 // _____________________________________________________________________________
0341 // DT
0342 GlobalPoint GeometryTranslator::calcDTSpecificPoint(const TriggerPrimitive& tp) const {
0343   const DTChamberId baseid(tp.detId<DTChamberId>());
0344   // do not use this pointer for anything other than creating a trig geom
0345   std::unique_ptr<DTChamber> chamb(const_cast<DTChamber*>(_geodt->chamber(baseid)));
0346   std::unique_ptr<DTTrigGeom> trig_geom(new DTTrigGeom(chamb.get(), false));
0347   chamb.release();  // release it here so no one gets funny ideas
0348   // super layer one is the theta superlayer in a DT chamber
0349   // station 4 does not have a theta super layer
0350   // the BTI index from the theta trigger is an OR of some BTI outputs
0351   // so, we choose the BTI that's in the middle of the group
0352   // as the BTI that we get theta from
0353   // TODO:::::>>> need to make sure this ordering doesn't flip under wheel sign
0354   const int NBTI_theta = ((baseid.station() != 4) ? trig_geom->nCell(2) : trig_geom->nCell(3));
0355   const int bti_group = tp.getDTData().theta_bti_group;
0356   const unsigned bti_actual = bti_group * NBTI_theta / 7 + NBTI_theta / 14 + 1;
0357   DTBtiId thetaBTI;
0358   if (baseid.station() != 4 && bti_group != -1) {
0359     thetaBTI = DTBtiId(baseid, 2, bti_actual);
0360   } else {
0361     // since this is phi oriented it'll give us theta in the middle
0362     // of the chamber
0363     thetaBTI = DTBtiId(baseid, 3, 1);
0364   }
0365   const GlobalPoint& theta_gp = trig_geom->CMSPosition(thetaBTI);
0366
0367   // local phi in sector -> global phi
0368   double phi = static_cast<double>(tp.getDTData().radialAngle) / 4096.0;  // 12 bits for 1 radian
0369   phi += tp.getDTData().sector * M_PI / 6.0;                              // add sector offset, sector is [0,11]
0370
0371   return GlobalPoint(GlobalPoint::Polar(theta_gp.theta(), phi, theta_gp.mag()));
0372 }
0373
0374 double GeometryTranslator::calcDTSpecificEta(const TriggerPrimitive& tp) const { return calcDTSpecificPoint(tp).eta(); }
0375
0376 double GeometryTranslator::calcDTSpecificPhi(const TriggerPrimitive& tp) const { return calcDTSpecificPoint(tp).phi(); }
0377
0378 // we have the bend except for station 3
0379 double GeometryTranslator::calcDTSpecificBend(const TriggerPrimitive& tp) const {
0380   int bend = tp.getDTData().bendingAngle;
0381   double bendf = bend / 512.0;
0382   return bendf;
0383 }