Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-27 14:07:46

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   const ME0DetId id(tp.detId<ME0DetId>());
0155   const ME0Chamber* chamber = _geome0->chamber(id);
0156   const ME0Layer* keylayer = chamber->layer(3);  // ME0 key layer is layer 3
0157   int partition = tp.getME0Data().partition;     // 'partition' is in half-roll unit
0158   int iroll = (partition >> 1) + 1;
0159   const ME0EtaPartition* roll = keylayer->etaPartition(iroll);
0160   assert(roll != nullptr);  // failed to get ME0 roll
0161   // See L1Trigger/ME0Trigger/src/ME0TriggerPseudoBuilder.cc
0162   int phiposition = tp.getME0Data().phiposition;  // 'phiposition' is in half-strip unit
0163   int istrip = (phiposition >> 1);
0164   int phiposition2 = (phiposition & 0x1);
0165   float centreOfStrip = istrip + 0.25 + phiposition2 * 0.5;
0166   const LocalPoint& lp = roll->centreOfStrip(centreOfStrip);
0167   const GlobalPoint& gp = roll->surface().toGlobal(lp);
0168   return gp;
0169 }
0170 
0171 double GeometryTranslator::calcME0SpecificEta(const TriggerPrimitive& tp) const {
0172   return getME0SpecificPoint(tp).eta();
0173 }
0174 
0175 double GeometryTranslator::calcME0SpecificPhi(const TriggerPrimitive& tp) const {
0176   return getME0SpecificPoint(tp).phi();
0177 }
0178 
0179 double GeometryTranslator::calcME0SpecificBend(const TriggerPrimitive& tp) const {
0180   return tp.getME0Data().deltaphi * (tp.getME0Data().bend == 0 ? 1 : -1);
0181 }
0182 
0183 // _____________________________________________________________________________
0184 // GEM
0185 GlobalPoint GeometryTranslator::getGEMSpecificPoint(const TriggerPrimitive& tp) const {
0186   const GEMDetId id(tp.detId<GEMDetId>());
0187   const GEMEtaPartition* roll = _geogem->etaPartition(id);
0188   assert(roll != nullptr);  // failed to get GEM roll
0189   //const uint16_t pad = tp.getGEMData().pad;
0190   // Use half-pad precision, + 0.5 at the end to get the center of the pad (pad starts from 0)
0191   const float pad = (0.5 * static_cast<float>(tp.getGEMData().pad_low + tp.getGEMData().pad_hi)) + 0.5f;
0192   const LocalPoint& lp = roll->centreOfPad(pad);
0193   const GlobalPoint& gp = roll->surface().toGlobal(lp);
0194   return gp;
0195 }
0196 
0197 double GeometryTranslator::calcGEMSpecificEta(const TriggerPrimitive& tp) const {
0198   return getGEMSpecificPoint(tp).eta();
0199 }
0200 
0201 double GeometryTranslator::calcGEMSpecificPhi(const TriggerPrimitive& tp) const {
0202   return getGEMSpecificPoint(tp).phi();
0203 }
0204 
0205 double GeometryTranslator::calcGEMSpecificBend(const TriggerPrimitive& tp) const { return 0.0; }
0206 
0207 // _____________________________________________________________________________
0208 // RPC
0209 GlobalPoint GeometryTranslator::getRPCSpecificPoint(const TriggerPrimitive& tp) const {
0210   const RPCDetId id(tp.detId<RPCDetId>());
0211   const RPCRoll* roll = _georpc->roll(id);
0212   assert(roll != nullptr);  // failed to get RPC roll
0213   //const int strip = static_cast<int>(tp.getRPCData().strip);
0214   // Use half-strip precision, - 0.5 at the end to get the center of the strip (strip starts from 1)
0215   const float strip = (0.5 * static_cast<float>(tp.getRPCData().strip_low + tp.getRPCData().strip_hi)) - 0.5f;
0216   const LocalPoint& lp = roll->centreOfStrip(strip);
0217   const GlobalPoint& gp = roll->surface().toGlobal(lp);
0218   return gp;
0219 }
0220 
0221 double GeometryTranslator::calcRPCSpecificEta(const TriggerPrimitive& tp) const {
0222   return getRPCSpecificPoint(tp).eta();
0223 }
0224 
0225 double GeometryTranslator::calcRPCSpecificPhi(const TriggerPrimitive& tp) const {
0226   return getRPCSpecificPoint(tp).phi();
0227 }
0228 
0229 // this function actually does nothing since RPC
0230 // hits are point-like objects
0231 double GeometryTranslator::calcRPCSpecificBend(const TriggerPrimitive& tp) const { return 0.0; }
0232 
0233 // _____________________________________________________________________________
0234 // CSC
0235 //
0236 // alot of this is transcription and consolidation of the CSC
0237 // global phi calculation code
0238 // this works directly with the geometry
0239 // rather than using the old phi luts
0240 GlobalPoint GeometryTranslator::getCSCSpecificPoint(const TriggerPrimitive& tp) const {
0241   const CSCDetId id(tp.detId<CSCDetId>());
0242   // we should change this to weak_ptrs at some point
0243   // requires introducing std::shared_ptrs to geometry
0244   std::unique_ptr<const CSCChamber> chamb(_geocsc->chamber(id));
0245   assert(chamb != nullptr);  // failed to get CSC chamber
0246   std::unique_ptr<const CSCLayerGeometry> layer_geom(chamb->layer(CSCConstants::KEY_ALCT_LAYER)->geometry());
0247   std::unique_ptr<const CSCLayer> layer(chamb->layer(CSCConstants::KEY_ALCT_LAYER));
0248 
0249   const uint16_t halfstrip = tp.getCSCData().strip;
0250   const uint16_t pattern = tp.getCSCData().pattern;
0251   const uint16_t keyWG = tp.getCSCData().keywire;
0252   //const unsigned maxStrips = layer_geom->numberOfStrips();
0253 
0254   // so we can extend this later
0255   // assume TMB2007 half-strips only as baseline
0256   double offset = 0.0;
0257   switch (1) {
0258     case 1:
0259       offset = CSCPatternBank::getLegacyPosition(pattern);
0260   }
0261   const unsigned halfstrip_offs = static_cast<unsigned>(0.5 + halfstrip + offset);
0262   const unsigned strip = halfstrip_offs / 2 + 1;  // geom starts from 1
0263 
0264   // the rough location of the hit at the ALCT key layer
0265   // we will refine this using the half strip information
0266   const LocalPoint& coarse_lp = layer_geom->stripWireGroupIntersection(strip, keyWG);
0267   const GlobalPoint& coarse_gp = layer->surface().toGlobal(coarse_lp);
0268 
0269   // the strip width/4.0 gives the offset of the half-strip
0270   // center with respect to the strip center
0271   const double hs_offset = layer_geom->stripPhiPitch() / 4.0;
0272 
0273   // determine handedness of the chamber
0274   const bool ccw = isCSCCounterClockwise(layer);
0275   // we need to subtract the offset of even half strips and add the odd ones
0276   const double phi_offset = ((halfstrip_offs % 2 ? 1 : -1) * (ccw ? -hs_offset : hs_offset));
0277 
0278   // the global eta calculation uses the middle of the strip
0279   // so no need to increment it
0280   const GlobalPoint final_gp(
0281       GlobalPoint::Polar(coarse_gp.theta(), (coarse_gp.phi().value() + phi_offset), coarse_gp.mag()));
0282 
0283   // We need to add in some notion of the 'error' on trigger primitives
0284   // like the width of the wire group by the width of the strip
0285   // or something similar
0286 
0287   // release ownership of the pointers
0288   chamb.release();
0289   layer_geom.release();
0290   layer.release();
0291 
0292   return final_gp;
0293 }
0294 
0295 double GeometryTranslator::calcCSCSpecificEta(const TriggerPrimitive& tp) const {
0296   return getCSCSpecificPoint(tp).eta();
0297 }
0298 
0299 double GeometryTranslator::calcCSCSpecificPhi(const TriggerPrimitive& tp) const {
0300   return getCSCSpecificPoint(tp).phi();
0301 }
0302 
0303 double GeometryTranslator::calcCSCSpecificBend(const TriggerPrimitive& tp) const { return tp.getCSCData().bend; }
0304 
0305 bool GeometryTranslator::isCSCCounterClockwise(const std::unique_ptr<const CSCLayer>& layer) const {
0306   const int nStrips = layer->geometry()->numberOfStrips();
0307   const double phi1 = layer->centerOfStrip(1).phi();
0308   const double phiN = layer->centerOfStrip(nStrips).phi();
0309   return ((std::abs(phi1 - phiN) < M_PI && phi1 >= phiN) || (std::abs(phi1 - phiN) >= M_PI && phi1 < phiN));
0310 }
0311 
0312 // _____________________________________________________________________________
0313 // DT
0314 GlobalPoint GeometryTranslator::calcDTSpecificPoint(const TriggerPrimitive& tp) const {
0315   const DTChamberId baseid(tp.detId<DTChamberId>());
0316   // do not use this pointer for anything other than creating a trig geom
0317   std::unique_ptr<DTChamber> chamb(const_cast<DTChamber*>(_geodt->chamber(baseid)));
0318   std::unique_ptr<DTTrigGeom> trig_geom(new DTTrigGeom(chamb.get(), false));
0319   chamb.release();  // release it here so no one gets funny ideas
0320   // super layer one is the theta superlayer in a DT chamber
0321   // station 4 does not have a theta super layer
0322   // the BTI index from the theta trigger is an OR of some BTI outputs
0323   // so, we choose the BTI that's in the middle of the group
0324   // as the BTI that we get theta from
0325   // TODO:::::>>> need to make sure this ordering doesn't flip under wheel sign
0326   const int NBTI_theta = ((baseid.station() != 4) ? trig_geom->nCell(2) : trig_geom->nCell(3));
0327   const int bti_group = tp.getDTData().theta_bti_group;
0328   const unsigned bti_actual = bti_group * NBTI_theta / 7 + NBTI_theta / 14 + 1;
0329   DTBtiId thetaBTI;
0330   if (baseid.station() != 4 && bti_group != -1) {
0331     thetaBTI = DTBtiId(baseid, 2, bti_actual);
0332   } else {
0333     // since this is phi oriented it'll give us theta in the middle
0334     // of the chamber
0335     thetaBTI = DTBtiId(baseid, 3, 1);
0336   }
0337   const GlobalPoint& theta_gp = trig_geom->CMSPosition(thetaBTI);
0338 
0339   // local phi in sector -> global phi
0340   double phi = static_cast<double>(tp.getDTData().radialAngle) / 4096.0;  // 12 bits for 1 radian
0341   phi += tp.getDTData().sector * M_PI / 6.0;                              // add sector offset, sector is [0,11]
0342 
0343   return GlobalPoint(GlobalPoint::Polar(theta_gp.theta(), phi, theta_gp.mag()));
0344 }
0345 
0346 double GeometryTranslator::calcDTSpecificEta(const TriggerPrimitive& tp) const { return calcDTSpecificPoint(tp).eta(); }
0347 
0348 double GeometryTranslator::calcDTSpecificPhi(const TriggerPrimitive& tp) const { return calcDTSpecificPoint(tp).phi(); }
0349 
0350 // we have the bend except for station 3
0351 double GeometryTranslator::calcDTSpecificBend(const TriggerPrimitive& tp) const {
0352   int bend = tp.getDTData().bendingAngle;
0353   double bendf = bend / 512.0;
0354   return bendf;
0355 }