File indexing completed on 2024-04-06 12:26:53
0001 #include <RecoMuon/DetLayers/src/MuonRPCDetLayerGeometryBuilder.h>
0002
0003 #include <DataFormats/MuonDetId/interface/RPCDetId.h>
0004 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0005 #include <RecoMuon/DetLayers/interface/MuRingForwardDoubleLayer.h>
0006 #include <RecoMuon/DetLayers/interface/MuRodBarrelLayer.h>
0007 #include <RecoMuon/DetLayers/interface/MuDetRing.h>
0008 #include <RecoMuon/DetLayers/interface/MuDetRod.h>
0009
0010 #include <Utilities/General/interface/precomputed_value_sort.h>
0011 #include <Geometry/CommonDetUnit/interface/DetSorting.h>
0012 #include "Utilities/BinningTools/interface/ClusterizingHistogram.h"
0013
0014 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0015
0016 #include <iostream>
0017
0018 using namespace std;
0019
0020 namespace rpcdetlayergeomsort {
0021 template <class T, class Scalar = typename T::Scalar>
0022 struct ExtractInnerRadius {
0023 typedef Scalar result_type;
0024 Scalar operator()(const T* p) const { return fabs(p->specificSurface().innerRadius()); }
0025 Scalar operator()(const T& p) const { return fabs(p.specificSurface().innerRadius()); }
0026 };
0027 }
0028
0029 MuonRPCDetLayerGeometryBuilder::~MuonRPCDetLayerGeometryBuilder() {}
0030
0031
0032 pair<vector<DetLayer*>, vector<DetLayer*> > MuonRPCDetLayerGeometryBuilder::buildEndcapLayers(const RPCGeometry& geo) {
0033 vector<DetLayer*> result[2];
0034
0035 for (int endcap = -1; endcap <= 1; endcap += 2) {
0036 int iendcap = (endcap == 1) ? 0 : 1;
0037
0038
0039 int firstStation = 1;
0040
0041
0042 for (int layer = RPCDetId::minLayerId; layer <= RPCDetId::maxLayerId; ++layer) {
0043 vector<int> rolls;
0044 std::vector<int> rings;
0045 int FirstStationRing = 1;
0046 rings.push_back(FirstStationRing);
0047 for (int roll = RPCDetId::minRollId + 1; roll <= RPCDetId::maxRollId; ++roll) {
0048 rolls.push_back(roll);
0049 }
0050
0051 MuRingForwardDoubleLayer* ringLayer = buildLayer(endcap, rings, firstStation, layer, rolls, geo);
0052 if (ringLayer)
0053 result[iendcap].push_back(ringLayer);
0054 }
0055
0056
0057 for (int layer = RPCDetId::minLayerId; layer <= RPCDetId::maxLayerId; ++layer) {
0058 vector<int> rolls;
0059 std::vector<int> rings;
0060 for (int ring = 2; ring <= 3; ++ring) {
0061 rings.push_back(ring);
0062 }
0063 for (int roll = RPCDetId::minRollId + 1; roll <= RPCDetId::maxRollId; ++roll) {
0064 rolls.push_back(roll);
0065 }
0066
0067 MuRingForwardDoubleLayer* ringLayer = buildLayer(endcap, rings, firstStation, layer, rolls, geo);
0068 if (ringLayer)
0069 result[iendcap].push_back(ringLayer);
0070 }
0071
0072
0073 for (int station = 2; station <= RPCDetId::maxStationId; ++station) {
0074 for (int layer = RPCDetId::minLayerId; layer <= RPCDetId::maxLayerId; ++layer) {
0075 vector<int> rolls;
0076 std::vector<int> rings;
0077 for (int ring = RPCDetId::minRingForwardId; ring <= RPCDetId::maxRingForwardId; ++ring) {
0078 rings.push_back(ring);
0079 }
0080 for (int roll = RPCDetId::minRollId + 1; roll <= RPCDetId::maxRollId; ++roll) {
0081 rolls.push_back(roll);
0082 }
0083
0084 MuRingForwardDoubleLayer* ringLayer = buildLayer(endcap, rings, station, layer, rolls, geo);
0085 if (ringLayer)
0086 result[iendcap].push_back(ringLayer);
0087 }
0088 }
0089 }
0090 pair<vector<DetLayer*>, vector<DetLayer*> > res_pair(result[0], result[1]);
0091 return res_pair;
0092 }
0093
0094 MuRingForwardDoubleLayer* MuonRPCDetLayerGeometryBuilder::buildLayer(
0095 int endcap, const std::vector<int>& rings, int station, int layer, vector<int>& rolls, const RPCGeometry& geo) {
0096 const std::string metname = "Muon|RPC|RecoMuon|RecoMuonDetLayers|MuonRPCDetLayerGeometryBuilder";
0097
0098 vector<const ForwardDetRing*> frontRings, backRings;
0099
0100 for (std::vector<int>::const_iterator ring = rings.begin(); ring < rings.end(); ++ring) {
0101 for (vector<int>::iterator roll = rolls.begin(); roll != rolls.end(); ++roll) {
0102 vector<const GeomDet*> frontDets, backDets;
0103 for (int sector = RPCDetId::minSectorForwardId; sector <= RPCDetId::maxSectorForwardId; ++sector) {
0104 for (int subsector = RPCDetId::minSubSectorForwardId; subsector <= RPCDetId::maxSectorForwardId; ++subsector) {
0105 RPCDetId rpcId(endcap, *ring, station, sector, layer, subsector, (*roll));
0106 bool isInFront = isFront(rpcId);
0107 const GeomDet* geomDet = geo.idToDet(rpcId);
0108 if (geomDet) {
0109 if (isInFront) {
0110 frontDets.push_back(geomDet);
0111 } else {
0112 backDets.push_back(geomDet);
0113 }
0114 LogTrace(metname) << "get RPC Endcap roll " << rpcId << (isInFront ? "front" : "back ")
0115 << " at R=" << geomDet->position().perp() << ", phi=" << geomDet->position().phi()
0116 << ", Z=" << geomDet->position().z();
0117 }
0118 }
0119 }
0120 if (!frontDets.empty()) {
0121 precomputed_value_sort(frontDets.begin(), frontDets.end(), geomsort::DetPhi());
0122 frontRings.push_back(new MuDetRing(frontDets));
0123 LogTrace(metname) << "New front ring with " << frontDets.size()
0124 << " chambers at z=" << frontRings.back()->position().z();
0125 }
0126 if (!backDets.empty()) {
0127 precomputed_value_sort(backDets.begin(), backDets.end(), geomsort::DetPhi());
0128 backRings.push_back(new MuDetRing(backDets));
0129 LogTrace(metname) << "New back ring with " << backDets.size()
0130 << " chambers at z=" << backRings.back()->position().z();
0131 }
0132 }
0133 }
0134
0135 MuRingForwardDoubleLayer* result = nullptr;
0136
0137 if (!backRings.empty() || !frontRings.empty()) {
0138 typedef rpcdetlayergeomsort::ExtractInnerRadius<ForwardDetRing, float> SortRingByInnerR;
0139 precomputed_value_sort(frontRings.begin(), frontRings.end(), SortRingByInnerR());
0140 precomputed_value_sort(backRings.begin(), backRings.end(), SortRingByInnerR());
0141
0142 result = new MuRingForwardDoubleLayer(frontRings, backRings);
0143 LogTrace(metname) << "New layer with " << frontRings.size() << " front rings and " << backRings.size()
0144 << " back rings, at Z " << result->position().z();
0145 }
0146
0147 return result;
0148 }
0149
0150 vector<DetLayer*> MuonRPCDetLayerGeometryBuilder::buildBarrelLayers(const RPCGeometry& geo) {
0151 const std::string metname = "Muon|RPC|RecoMuon|RecoMuonDetLayers|MuonRPCDetLayerGeometryBuilder";
0152
0153 vector<DetLayer*> detlayers;
0154 vector<MuRodBarrelLayer*> result;
0155 int region = 0;
0156
0157 for (int station = RPCDetId::minStationId; station <= RPCDetId::maxStationId; station++) {
0158 vector<const GeomDet*> geomDets;
0159
0160 for (int layer = RPCDetId::minLayerId; layer <= RPCDetId::maxLayerId; ++layer) {
0161 for (int sector = RPCDetId::minSectorId; sector <= RPCDetId::maxSectorId; sector++) {
0162 for (int subsector = RPCDetId::minSubSectorId; subsector <= RPCDetId::maxSubSectorId; subsector++) {
0163 for (int wheel = RPCDetId::minRingBarrelId; wheel <= RPCDetId::maxRingBarrelId; wheel++) {
0164 for (int roll = RPCDetId::minRollId + 1; roll <= RPCDetId::maxRollId; roll++) {
0165 const GeomDet* geomDet = geo.idToDet(RPCDetId(region, wheel, station, sector, layer, subsector, roll));
0166 if (geomDet) {
0167 geomDets.push_back(geomDet);
0168 LogTrace(metname) << "get RPC Barrel roll "
0169 << RPCDetId(region, wheel, station, sector, layer, subsector, roll)
0170 << " at R=" << geomDet->position().perp() << ", phi=" << geomDet->position().phi();
0171 }
0172 }
0173 }
0174 }
0175 }
0176 }
0177 makeBarrelLayers(geomDets, result);
0178 }
0179
0180 for (vector<MuRodBarrelLayer*>::const_iterator it = result.begin(); it != result.end(); it++)
0181 detlayers.push_back((DetLayer*)(*it));
0182
0183 return detlayers;
0184 }
0185
0186 void MuonRPCDetLayerGeometryBuilder::makeBarrelLayers(vector<const GeomDet*>& geomDets,
0187 vector<MuRodBarrelLayer*>& result) {
0188 const std::string metname = "Muon|RPC|RecoMuon|RecoMuonDetLayers|MuonRPCDetLayerGeometryBuilder";
0189
0190
0191 precomputed_value_sort(geomDets.begin(), geomDets.end(), geomsort::DetR());
0192
0193
0194 float resolution(25);
0195 float r0 = float(geomDets.front()->position().perp());
0196 float rMin = -float(resolution);
0197 float rMax = float(geomDets.back()->position().perp()) - r0 + resolution;
0198
0199 ClusterizingHistogram hisR(int((rMax - rMin) / resolution) + 1, rMin, rMax);
0200
0201 vector<const GeomDet*>::iterator first = geomDets.begin();
0202 vector<const GeomDet*>::iterator last = geomDets.end();
0203
0204 for (vector<const GeomDet*>::iterator i = first; i != last; i++) {
0205 hisR.fill(float((*i)->position().perp()) - r0);
0206 LogTrace(metname) << "R " << float((*i)->position().perp()) - r0;
0207 }
0208 vector<float> rClust = hisR.clusterize(resolution);
0209
0210
0211
0212 vector<const GeomDet*>::iterator layerStart = first;
0213 vector<const GeomDet*>::iterator separ = first;
0214
0215 for (unsigned int i = 0; i < rClust.size(); i++) {
0216 float rSepar;
0217 if (i < rClust.size() - 1) {
0218 rSepar = (rClust[i] + rClust[i + 1]) / 2.f;
0219 } else {
0220 rSepar = rMax;
0221 }
0222
0223
0224
0225 while (separ < last && float((*separ)->position().perp()) - r0 < rSepar) {
0226
0227 separ++;
0228 }
0229
0230 if (int(separ - layerStart) > 0) {
0231
0232 vector<const DetRod*> rods;
0233 vector<const GeomDet*> layerDets(layerStart, separ);
0234 makeBarrelRods(layerDets, rods);
0235
0236 if (!rods.empty()) {
0237 result.push_back(new MuRodBarrelLayer(rods));
0238 LogTrace(metname) << " New MuRodBarrelLayer with " << rods.size() << " rods, at R "
0239 << result.back()->specificSurface().radius();
0240 }
0241 }
0242 layerStart = separ;
0243 }
0244 }
0245
0246 void MuonRPCDetLayerGeometryBuilder::makeBarrelRods(vector<const GeomDet*>& geomDets, vector<const DetRod*>& result) {
0247 const std::string metname = "Muon|RPC|RecoMuon|RecoMuonDetLayers|MuonRPCDetLayerGeometryBuilder";
0248
0249
0250 precomputed_value_sort(geomDets.begin(), geomDets.end(), geomsort::DetPhi());
0251
0252
0253 float resolution(0.01);
0254 float phi0 = float(geomDets.front()->position().phi());
0255 float phiMin = -float(resolution);
0256 float phiMax = float(geomDets.back()->position().phi()) - phi0 + resolution;
0257
0258 ClusterizingHistogram hisPhi(int((phiMax - phiMin) / resolution) + 1, phiMin, phiMax);
0259
0260 vector<const GeomDet*>::iterator first = geomDets.begin();
0261 vector<const GeomDet*>::iterator last = geomDets.end();
0262
0263 for (vector<const GeomDet*>::iterator i = first; i != last; i++) {
0264 hisPhi.fill(float((*i)->position().phi()) - phi0);
0265 LogTrace(metname) << "C " << float((*i)->position().phi()) - phi0;
0266 }
0267 vector<float> phiClust = hisPhi.clusterize(resolution);
0268
0269
0270
0271 vector<const GeomDet*>::iterator rodStart = first;
0272 vector<const GeomDet*>::iterator separ = first;
0273
0274 for (unsigned int i = 0; i < phiClust.size(); i++) {
0275 float phiSepar;
0276 if (i < phiClust.size() - 1) {
0277 phiSepar = (phiClust[i] + phiClust[i + 1]) / 2.f;
0278 } else {
0279 phiSepar = phiMax;
0280 }
0281
0282
0283
0284 while (separ < last && float((*separ)->position().phi()) - phi0 < phiSepar) {
0285
0286 separ++;
0287 }
0288
0289 if (int(separ - rodStart) > 0) {
0290 result.push_back(new MuDetRod(rodStart, separ));
0291 LogTrace(metname) << " New MuDetRod with " << int(separ - rodStart)
0292 << " rolls at R=" << (*rodStart)->position().perp()
0293 << ", phi=" << float((*rodStart)->position().phi());
0294 }
0295 rodStart = separ;
0296 }
0297 }
0298
0299 bool MuonRPCDetLayerGeometryBuilder::isFront(const RPCDetId& rpcId) {
0300
0301
0302
0303 bool result = false;
0304 int ring = rpcId.ring();
0305 int station = rpcId.station();
0306
0307 if (ring == 1 && station > 1) {
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318 result = (rpcId.subsector() != 2);
0319 if (rpcId.sector() % 2 == 0)
0320 result = !result;
0321 return result;
0322 } else {
0323
0324 result = (rpcId.subsector() % 2 == 0);
0325 }
0326 return result;
0327 }