File indexing completed on 2022-11-22 00:47:27
0001 #include "FWCore/Framework/interface/ModuleFactory.h"
0002 #include "FWCore/Framework/interface/ESProducer.h"
0003
0004 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0005 #include "DataFormats/TrackerCommon/interface/TrackerDetSide.h"
0006 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0007 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0008 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0009 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0010 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0011
0012 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0013 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0014
0015
0016 #include "RecoTracker/MkFit/interface/MkFitGeometry.h"
0017 #include "RecoTracker/MkFitCore/interface/ConfigWrapper.h"
0018 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0019 #include "RecoTracker/MkFitCore/interface/IterationConfig.h"
0020 #include "RecoTracker/MkFitCMS/interface/LayerNumberConverter.h"
0021
0022 #include <sstream>
0023
0024
0025
0026
0027
0028 class MkFitGeometryESProducer : public edm::ESProducer {
0029 public:
0030 MkFitGeometryESProducer(const edm::ParameterSet &iConfig);
0031
0032 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0033
0034 std::unique_ptr<MkFitGeometry> produce(const TrackerRecoGeometryRecord &iRecord);
0035
0036 private:
0037 struct GapCollector {
0038 struct Interval {
0039 float x, y;
0040 };
0041
0042 void reset_current() { m_current = {std::numeric_limits<float>::max(), -std::numeric_limits<float>::max()}; }
0043 void extend_current(float q) {
0044 m_current.x = std::min(m_current.x, q);
0045 m_current.y = std::max(m_current.y, q);
0046 }
0047 void add_current() { add_interval(m_current.x, m_current.y); }
0048
0049 void add_interval(float x, float y);
0050
0051 void sqrt_elements();
0052 bool find_gap(Interval &itvl, float eps);
0053 void print_gaps(std::ostream &ostr);
0054
0055 std::list<Interval> m_coverage;
0056 Interval m_current;
0057 };
0058 typedef std::unordered_map<int, GapCollector> layer_gap_map_t;
0059
0060 void considerPoint(const GlobalPoint &gp, mkfit::LayerInfo &lay_info);
0061 void fillShapeAndPlacement(const GeomDet *det, mkfit::TrackerInfo &trk_info, layer_gap_map_t *lgc_map = nullptr);
0062 void addPixBGeometry(mkfit::TrackerInfo &trk_info);
0063 void addPixEGeometry(mkfit::TrackerInfo &trk_info);
0064 void addTIBGeometry(mkfit::TrackerInfo &trk_info);
0065 void addTOBGeometry(mkfit::TrackerInfo &trk_info);
0066 void addTIDGeometry(mkfit::TrackerInfo &trk_info);
0067 void addTECGeometry(mkfit::TrackerInfo &trk_info);
0068
0069 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0070 edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> ttopoToken_;
0071 edm::ESGetToken<GeometricSearchTracker, TrackerRecoGeometryRecord> trackerToken_;
0072
0073 const TrackerTopology *trackerTopo_ = nullptr;
0074 const TrackerGeometry *trackerGeom_ = nullptr;
0075 mkfit::LayerNumberConverter layerNrConv_ = {mkfit::TkLayout::phase1};
0076 };
0077
0078 MkFitGeometryESProducer::MkFitGeometryESProducer(const edm::ParameterSet &iConfig) {
0079 auto cc = setWhatProduced(this);
0080 geomToken_ = cc.consumes();
0081 ttopoToken_ = cc.consumes();
0082 trackerToken_ = cc.consumes();
0083 }
0084
0085 void MkFitGeometryESProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0086 edm::ParameterSetDescription desc;
0087 descriptions.addWithDefaultLabel(desc);
0088 }
0089
0090
0091
0092 void MkFitGeometryESProducer::GapCollector::add_interval(float x, float y) {
0093 if (x > y)
0094 std::swap(x, y);
0095 bool handled = false;
0096 for (auto i = m_coverage.begin(); i != m_coverage.end(); ++i) {
0097 if (y < i->x) {
0098 m_coverage.insert(i, {x, y});
0099 handled = true;
0100 break;
0101 } else if (x > i->y) {
0102 continue;
0103 } else if (x < i->x) {
0104 i->x = x;
0105 handled = true;
0106 break;
0107 } else if (y > i->y) {
0108 i->y = y;
0109
0110 auto j = i;
0111 ++j;
0112 if (j != m_coverage.end() && i->y >= j->x) {
0113 i->y = j->y;
0114 m_coverage.erase(j);
0115 }
0116 handled = true;
0117 break;
0118 } else {
0119 handled = true;
0120 break;
0121 }
0122 }
0123 if (!handled) {
0124 m_coverage.push_back({x, y});
0125 }
0126 }
0127
0128 void MkFitGeometryESProducer::GapCollector::sqrt_elements() {
0129 for (auto &itvl : m_coverage) {
0130 itvl.x = std::sqrt(itvl.x);
0131 itvl.y = std::sqrt(itvl.y);
0132 }
0133 }
0134
0135 bool MkFitGeometryESProducer::GapCollector::find_gap(Interval &itvl, float eps) {
0136 auto i = m_coverage.begin();
0137 while (i != m_coverage.end()) {
0138 auto j = i;
0139 ++j;
0140 if (j != m_coverage.end()) {
0141 if (j->x - i->y > eps) {
0142 itvl = {i->y, j->x};
0143 return true;
0144 }
0145 i = j;
0146 } else {
0147 break;
0148 }
0149 }
0150 return false;
0151 }
0152
0153 void MkFitGeometryESProducer::GapCollector::print_gaps(std::ostream &ostr) {
0154 auto i = m_coverage.begin();
0155 while (i != m_coverage.end()) {
0156 auto j = i;
0157 ++j;
0158 if (j != m_coverage.end()) {
0159 ostr << "(" << i->y << ", " << j->x << ")->" << j->x - i->y;
0160 i = j;
0161 } else {
0162 break;
0163 }
0164 }
0165 }
0166
0167
0168
0169 void MkFitGeometryESProducer::considerPoint(const GlobalPoint &gp, mkfit::LayerInfo &li) {
0170
0171 float r = gp.perp2(), z = gp.z();
0172 li.extend_limits(r, z);
0173 }
0174
0175 void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det,
0176 mkfit::TrackerInfo &trk_info,
0177 layer_gap_map_t *lgc_map) {
0178 DetId detid = det->geographicalId();
0179
0180 float xy[4][2];
0181 float dz;
0182 const Bounds *b = &((det->surface()).bounds());
0183
0184 if (const TrapezoidalPlaneBounds *b2 = dynamic_cast<const TrapezoidalPlaneBounds *>(b)) {
0185
0186 std::array<const float, 4> const &par = b2->parameters();
0187 xy[0][0] = -par[0];
0188 xy[0][1] = -par[3];
0189 xy[1][0] = -par[1];
0190 xy[1][1] = par[3];
0191 xy[2][0] = par[1];
0192 xy[2][1] = par[3];
0193 xy[3][0] = par[0];
0194 xy[3][1] = -par[3];
0195 dz = par[2];
0196
0197 #ifdef DUMP_MKF_GEO
0198 printf("TRAP 0x%x %f %f %f %f ", detid.rawId(), par[0], par[1], par[2], par[3]);
0199 #endif
0200 } else if (const RectangularPlaneBounds *b2 = dynamic_cast<const RectangularPlaneBounds *>(b)) {
0201
0202 float dx = b2->width() * 0.5;
0203 float dy = b2->length() * 0.5;
0204 xy[0][0] = -dx;
0205 xy[0][1] = -dy;
0206 xy[1][0] = -dx;
0207 xy[1][1] = dy;
0208 xy[2][0] = dx;
0209 xy[2][1] = dy;
0210 xy[3][0] = dx;
0211 xy[3][1] = -dy;
0212 dz = b2->thickness() * 0.5;
0213
0214 #ifdef DUMP_MKF_GEO
0215 printf("RECT 0x%x %f %f %f ", detid.rawId(), dx, dy, dz);
0216 #endif
0217 } else {
0218 throw cms::Exception("UnimplementedFeature") << "unsupported Bounds class";
0219 }
0220
0221 const bool useMatched = false;
0222 int lay =
0223 layerNrConv_.convertLayerNumber(detid.subdetId(),
0224 trackerTopo_->layer(detid),
0225 useMatched,
0226 trackerTopo_->isStereo(detid),
0227 trackerTopo_->side(detid) == static_cast<unsigned>(TrackerDetSide::PosEndcap));
0228 #ifdef DUMP_MKF_GEO
0229 printf(" subdet=%d layer=%d side=%d is_stereo=%d --> mkflayer=%d\n",
0230 detid.subdetId(),
0231 trackerTopo_->layer(detid),
0232 trackerTopo_->side(detid),
0233 trackerTopo_->isStereo(detid),
0234 lay);
0235 #endif
0236
0237 mkfit::LayerInfo &layer_info = trk_info.layer_nc(lay);
0238 if (lgc_map) {
0239 (*lgc_map)[lay].reset_current();
0240 }
0241 for (int i = 0; i < 4; ++i) {
0242 Local3DPoint lp1(xy[i][0], xy[i][1], -dz);
0243 Local3DPoint lp2(xy[i][0], xy[i][1], dz);
0244 GlobalPoint gp1 = det->surface().toGlobal(lp1);
0245 GlobalPoint gp2 = det->surface().toGlobal(lp2);
0246 considerPoint(gp1, layer_info);
0247 considerPoint(gp2, layer_info);
0248 if (lgc_map) {
0249 (*lgc_map)[lay].extend_current(gp1.perp2());
0250 (*lgc_map)[lay].extend_current(gp2.perp2());
0251 }
0252 }
0253 if (lgc_map) {
0254 (*lgc_map)[lay].add_current();
0255 }
0256
0257 const auto &p = det->position();
0258 auto z = det->rotation().z();
0259 auto x = det->rotation().x();
0260 layer_info.register_module({{p.x(), p.y(), p.z()}, {z.x(), z.y(), z.z()}, {x.x(), x.y(), x.z()}, detid.rawId()});
0261
0262 layer_info.set_subdet(detid.subdetId());
0263 layer_info.set_is_pixel(detid.subdetId() <= 2);
0264 layer_info.set_is_stereo(trackerTopo_->isStereo(detid));
0265 }
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278 void MkFitGeometryESProducer::addPixBGeometry(mkfit::TrackerInfo &trk_info) {
0279 #ifdef DUMP_MKF_GEO
0280 printf("\n*** addPixBGeometry\n\n");
0281 #endif
0282 for (auto &det : trackerGeom_->detsPXB()) {
0283 fillShapeAndPlacement(det, trk_info);
0284 }
0285 }
0286
0287 void MkFitGeometryESProducer::addPixEGeometry(mkfit::TrackerInfo &trk_info) {
0288 #ifdef DUMP_MKF_GEO
0289 printf("\n*** addPixEGeometry\n\n");
0290 #endif
0291 for (auto &det : trackerGeom_->detsPXF()) {
0292 fillShapeAndPlacement(det, trk_info);
0293 }
0294 }
0295
0296 void MkFitGeometryESProducer::addTIBGeometry(mkfit::TrackerInfo &trk_info) {
0297 #ifdef DUMP_MKF_GEO
0298 printf("\n*** addTIBGeometry\n\n");
0299 #endif
0300 for (auto &det : trackerGeom_->detsTIB()) {
0301 fillShapeAndPlacement(det, trk_info);
0302 }
0303 }
0304
0305 void MkFitGeometryESProducer::addTOBGeometry(mkfit::TrackerInfo &trk_info) {
0306 #ifdef DUMP_MKF_GEO
0307 printf("\n*** addTOBGeometry\n\n");
0308 #endif
0309 for (auto &det : trackerGeom_->detsTOB()) {
0310 fillShapeAndPlacement(det, trk_info);
0311 }
0312 }
0313
0314 void MkFitGeometryESProducer::addTIDGeometry(mkfit::TrackerInfo &trk_info) {
0315 #ifdef DUMP_MKF_GEO
0316 printf("\n*** addTIDGeometry\n\n");
0317 #endif
0318 for (auto &det : trackerGeom_->detsTID()) {
0319 fillShapeAndPlacement(det, trk_info);
0320 }
0321 }
0322
0323 void MkFitGeometryESProducer::addTECGeometry(mkfit::TrackerInfo &trk_info) {
0324 #ifdef DUMP_MKF_GEO
0325 printf("\n*** addTECGeometry\n\n");
0326 #endif
0327
0328 layer_gap_map_t lgc_map;
0329 for (auto &det : trackerGeom_->detsTEC()) {
0330 fillShapeAndPlacement(det, trk_info, &lgc_map);
0331 }
0332
0333 std::ostringstream ostr;
0334 ostr << "addTECGeometry() gap report:\n";
0335 GapCollector::Interval itvl;
0336 for (auto &[layer, gcol] : lgc_map) {
0337 gcol.sqrt_elements();
0338 if (gcol.find_gap(itvl, 0.5)) {
0339 ostr << " layer: " << layer << ", gap: " << itvl.x << " -> " << itvl.y << " width = " << itvl.y - itvl.x << "\n";
0340 ostr << " all gaps: ";
0341 gcol.print_gaps(ostr);
0342 ostr << "\n";
0343 trk_info.layer_nc(layer).set_r_hole_range(itvl.x, itvl.y);
0344 }
0345 }
0346 edm::LogVerbatim("MkFitGeometryESProducer") << ostr.str();
0347 }
0348
0349
0350
0351 namespace {
0352 const float phase1QBins[] = {
0353
0354 2.0, 2.0, 2.0, 2.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
0355
0356 1.0, 1.0, 1.0, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5,
0357
0358 1.0, 1.0, 1.0, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5
0359 };
0360 const float phase2QBins[] = {
0361
0362
0363 2.0, 2.0, 2.0, 2.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0,
0364
0365 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6,
0366
0367 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6
0368 };
0369 }
0370
0371
0372
0373 std::unique_ptr<MkFitGeometry> MkFitGeometryESProducer::produce(const TrackerRecoGeometryRecord &iRecord) {
0374 auto trackerInfo = std::make_unique<mkfit::TrackerInfo>();
0375
0376 trackerGeom_ = &iRecord.get(geomToken_);
0377 trackerTopo_ = &iRecord.get(ttopoToken_);
0378
0379 const float *qBinDefaults = nullptr;
0380
0381
0382 if (trackerGeom_->isThere(GeomDetEnumerators::P1PXB) || trackerGeom_->isThere(GeomDetEnumerators::P1PXEC)) {
0383 edm::LogInfo("MkFitGeometryESProducer") << "Extracting PhaseI geometry";
0384 trackerInfo->create_layers(18, 27, 27);
0385 qBinDefaults = phase1QBins;
0386 } else if (trackerGeom_->isThere(GeomDetEnumerators::P2PXB) || trackerGeom_->isThere(GeomDetEnumerators::P2PXEC) ||
0387 trackerGeom_->isThere(GeomDetEnumerators::P2OTB) || trackerGeom_->isThere(GeomDetEnumerators::P2OTEC)) {
0388 edm::LogInfo("MkFitGeometryESProducer") << "Extracting PhaseII geometry";
0389 layerNrConv_.reset(mkfit::TkLayout::phase2);
0390 trackerInfo->create_layers(16, 22, 22);
0391 qBinDefaults = phase2QBins;
0392 } else {
0393 throw cms::Exception("UnimplementedFeature") << "unsupported / unknowen geometry version";
0394 }
0395
0396
0397 for (int i = 0; i < trackerInfo->n_layers(); ++i) {
0398 auto &li = trackerInfo->layer_nc(i);
0399 li.set_limits(
0400 std::numeric_limits<float>::max(), 0, std::numeric_limits<float>::max(), -std::numeric_limits<float>::max());
0401 li.reserve_modules(256);
0402 }
0403
0404 addPixBGeometry(*trackerInfo);
0405 addPixEGeometry(*trackerInfo);
0406 addTIBGeometry(*trackerInfo);
0407 addTIDGeometry(*trackerInfo);
0408 addTOBGeometry(*trackerInfo);
0409 addTECGeometry(*trackerInfo);
0410
0411
0412 unsigned int n_mod = 0;
0413 for (int i = 0; i < trackerInfo->n_layers(); ++i) {
0414 auto &li = trackerInfo->layer_nc(i);
0415 li.set_r_in_out(std::sqrt(li.rin()), std::sqrt(li.rout()));
0416 li.set_propagate_to(li.is_barrel() ? li.r_mean() : li.z_mean());
0417 li.set_q_bin(qBinDefaults[i]);
0418 unsigned int maxsid = li.shrink_modules();
0419
0420 n_mod += maxsid;
0421
0422
0423 assert(maxsid < 1u << 13);
0424 assert(n_mod > 0);
0425 }
0426 #ifdef DUMP_MKF_GEO
0427 printf("Total number of modules %u, 14-bits fit up to %u modules\n", n_mod, 1u << 13);
0428 #endif
0429
0430 return std::make_unique<MkFitGeometry>(iRecord.get(geomToken_),
0431 iRecord.get(trackerToken_),
0432 iRecord.get(ttopoToken_),
0433 std::move(trackerInfo),
0434 layerNrConv_);
0435 }
0436
0437 DEFINE_FWK_EVENTSETUP_MODULE(MkFitGeometryESProducer);