DetLayers/src/MuonRPCDetLayerGeometryBuilder.cc

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 }  // namespace rpcdetlayergeomsort
0028
0029 MuonRPCDetLayerGeometryBuilder::~MuonRPCDetLayerGeometryBuilder() {}
0030
0031 // Builds the forward (first) and backward (second) layers
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;  // +1: forward, -1: backward
0037
0038     // ME 1
0039     int firstStation = 1;
0040
0041     // ME 1/1
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     // ME 1/2 and ME1/3
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     // ME 2 and ME 3
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   //Sort in R
0191   precomputed_value_sort(geomDets.begin(), geomDets.end(), geomsort::DetR());
0192
0193   // Clusterize in phi - phi0
0194   float resolution(25);  // cm
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   // LogTrace(metname) << "     Found " << phiClust.size() << " clusters in Phi, ";
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     // LogTrace(metname) << "       cluster " << i
0224     // << " phisepar " << phiSepar <<endl;
0225     while (separ < last && float((*separ)->position().perp()) - r0 < rSepar) {
0226       // LogTrace(metname) << "         roll at dphi:  " << float((*separ)->position().phi())-phi0;
0227       separ++;
0228     }
0229
0230     if (int(separ - layerStart) > 0) {
0231       // we have a layer in R.  Now separate it into rods
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   //Sort in phi
0250   precomputed_value_sort(geomDets.begin(), geomDets.end(), geomsort::DetPhi());
0251
0252   // Clusterize in phi - phi0
0253   float resolution(0.01);  // rad
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   // LogTrace(metname) << "     Found " << phiClust.size() << " clusters in Phi, ";
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     // LogTrace(metname) << "       cluster " << i
0283     // << " phisepar " << phiSepar <<endl;
0284     while (separ < last && float((*separ)->position().phi()) - phi0 < phiSepar) {
0285       // LogTrace(metname) << "         roll at dphi:  " << float((*separ)->position().phi())-phi0;
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   // ME1/2 is always in back
0301   //  if(rpcId.station() == 1 && rpcId.ring() == 2)  return false;
0302
0303   bool result = false;
0304   int ring = rpcId.ring();
0305   int station = rpcId.station();
0306   // 20 degree rings are a little weird! not anymore from 17x
0307   if (ring == 1 && station > 1) {
0308     // RE2/1 RE3/1  Upscope Geometry
0309     /* goes (sector) (subsector)            1/3
0310     1 1 back   // front
0311     1 2 front  // back
0312     1 3 front  // front
0313     2 1 front  // back
0314     2 2 back   // from
0315     2 3 back   // back
0316
0317     */
0318     result = (rpcId.subsector() != 2);
0319     if (rpcId.sector() % 2 == 0)
0320       result = !result;
0321     return result;
0322   } else {
0323     // 10 degree rings have odd subsectors in front
0324     result = (rpcId.subsector() % 2 == 0);
0325   }
0326   return result;
0327 }