MkFit/plugins/MkFitGeometryESProducer.cc

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 // mkFit includes
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 class MkFitGeometryESProducer : public edm::ESProducer {
0027 public:
0028   MkFitGeometryESProducer(const edm::ParameterSet &iConfig);
0029
0030   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0031
0032   std::unique_ptr<MkFitGeometry> produce(const TrackerRecoGeometryRecord &iRecord);
0033
0034 private:
0035   struct GapCollector {
0036     struct Interval {
0037       float x, y;
0038     };
0039
0040     void reset_current() { m_current = {std::numeric_limits<float>::max(), -std::numeric_limits<float>::max()}; }
0041     void extend_current(float q) {
0042       m_current.x = std::min(m_current.x, q);
0043       m_current.y = std::max(m_current.y, q);
0044     }
0045     void add_current() { add_interval(m_current.x, m_current.y); }
0046
0047     void add_interval(float x, float y);
0048
0049     void sqrt_elements();
0050     bool find_gap(Interval &itvl, float eps);
0051     void print_gaps(std::ostream &ostr);
0052
0053     std::list<Interval> m_coverage;
0054     Interval m_current;
0055   };
0056   typedef std::unordered_map<int, GapCollector> layer_gap_map_t;
0057
0058   void considerPoint(const GlobalPoint &gp, mkfit::LayerInfo &lay_info);
0059   void fillShapeAndPlacement(const GeomDet *det, mkfit::TrackerInfo &trk_info, layer_gap_map_t *lgc_map = nullptr);
0060   void addPixBGeometry(mkfit::TrackerInfo &trk_info);
0061   void addPixEGeometry(mkfit::TrackerInfo &trk_info);
0062   void addTIBGeometry(mkfit::TrackerInfo &trk_info);
0063   void addTOBGeometry(mkfit::TrackerInfo &trk_info);
0064   void addTIDGeometry(mkfit::TrackerInfo &trk_info);
0065   void addTECGeometry(mkfit::TrackerInfo &trk_info);
0066
0067   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0068   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> ttopoToken_;
0069   edm::ESGetToken<GeometricSearchTracker, TrackerRecoGeometryRecord> trackerToken_;
0070
0071   const TrackerTopology *trackerTopo_ = nullptr;
0072   const TrackerGeometry *trackerGeom_ = nullptr;
0073   mkfit::LayerNumberConverter layerNrConv_ = {mkfit::TkLayout::phase1};
0074 };
0075
0076 MkFitGeometryESProducer::MkFitGeometryESProducer(const edm::ParameterSet &iConfig) {
0077   auto cc = setWhatProduced(this);
0078   geomToken_ = cc.consumes();
0079   ttopoToken_ = cc.consumes();
0080   trackerToken_ = cc.consumes();
0081 }
0082
0083 void MkFitGeometryESProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0084   edm::ParameterSetDescription desc;
0085   descriptions.addWithDefaultLabel(desc);
0086 }
0087
0088 //------------------------------------------------------------------------------
0089
0090 void MkFitGeometryESProducer::GapCollector::add_interval(float x, float y) {
0091   if (x > y)
0092     std::swap(x, y);
0093   bool handled = false;
0094   for (auto i = m_coverage.begin(); i != m_coverage.end(); ++i) {
0095     if (y < i->x) {  // fully on 'left'
0096       m_coverage.insert(i, {x, y});
0097       handled = true;
0098       break;
0099     } else if (x > i->y) {  // fully on 'right'
0100       continue;
0101     } else if (x < i->x) {  // sticking out on 'left'
0102       i->x = x;
0103       handled = true;
0104       break;
0105     } else if (y > i->y) {  // sticking out on 'right'
0106       i->y = y;
0107       // check for overlap with the next interval, potentially merge
0108       auto j = i;
0109       ++j;
0110       if (j != m_coverage.end() && i->y >= j->x) {
0111         i->y = j->y;
0112         m_coverage.erase(j);
0113       }
0114       handled = true;
0115       break;
0116     } else {  // contained in current interval
0117       handled = true;
0118       break;
0119     }
0120   }
0121   if (!handled) {
0122     m_coverage.push_back({x, y});
0123   }
0124 }
0125
0126 void MkFitGeometryESProducer::GapCollector::sqrt_elements() {
0127   for (auto &itvl : m_coverage) {
0128     itvl.x = std::sqrt(itvl.x);
0129     itvl.y = std::sqrt(itvl.y);
0130   }
0131 }
0132
0133 bool MkFitGeometryESProducer::GapCollector::find_gap(Interval &itvl, float eps) {
0134   auto i = m_coverage.begin();
0135   while (i != m_coverage.end()) {
0136     auto j = i;
0137     ++j;
0138     if (j != m_coverage.end()) {
0139       if (j->x - i->y > eps) {
0140         itvl = {i->y, j->x};
0141         return true;
0142       }
0143       i = j;
0144     } else {
0145       break;
0146     }
0147   }
0148   return false;
0149 }
0150
0151 void MkFitGeometryESProducer::GapCollector::print_gaps(std::ostream &ostr) {
0152   auto i = m_coverage.begin();
0153   while (i != m_coverage.end()) {
0154     auto j = i;
0155     ++j;
0156     if (j != m_coverage.end()) {
0157       ostr << "(" << i->y << ", " << j->x << ")->" << j->x - i->y;
0158       i = j;
0159     } else {
0160       break;
0161     }
0162   }
0163 }
0164
0165 //------------------------------------------------------------------------------
0166
0167 void MkFitGeometryESProducer::considerPoint(const GlobalPoint &gp, mkfit::LayerInfo &li) {
0168   // Use radius squared during bounding-region search.
0169   float r = gp.perp2(), z = gp.z();
0170   li.extend_limits(r, z);
0171 }
0172
0173 void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det,
0174                                                     mkfit::TrackerInfo &trk_info,
0175                                                     layer_gap_map_t *lgc_map) {
0176   DetId detid = det->geographicalId();
0177
0178   float xy[4][2];
0179   float dz;
0180   const Bounds *b = &((det->surface()).bounds());
0181
0182   if (const TrapezoidalPlaneBounds *b2 = dynamic_cast<const TrapezoidalPlaneBounds *>(b)) {
0183     // See sec. "TrapezoidalPlaneBounds parameters" in doc/reco-geom-notes.txt
0184     std::array<const float, 4> const &par = b2->parameters();
0185     xy[0][0] = -par[0];
0186     xy[0][1] = -par[3];
0187     xy[1][0] = -par[1];
0188     xy[1][1] = par[3];
0189     xy[2][0] = par[1];
0190     xy[2][1] = par[3];
0191     xy[3][0] = par[0];
0192     xy[3][1] = -par[3];
0193     dz = par[2];
0194
0195     // printf("TRAP 0x%x %f %f %f %f\n", detid.rawId(), par[0], par[1], par[2], par[3]);
0196   } else if (const RectangularPlaneBounds *b2 = dynamic_cast<const RectangularPlaneBounds *>(b)) {
0197     // Rectangular
0198     float dx = b2->width() * 0.5;   // half width
0199     float dy = b2->length() * 0.5;  // half length
0200     xy[0][0] = -dx;
0201     xy[0][1] = -dy;
0202     xy[1][0] = -dx;
0203     xy[1][1] = dy;
0204     xy[2][0] = dx;
0205     xy[2][1] = dy;
0206     xy[3][0] = dx;
0207     xy[3][1] = -dy;
0208     dz = b2->thickness() * 0.5;  // half thickness
0209
0210     // printf("RECT 0x%x %f %f %f\n", detid.rawId(), dx, dy, dz);
0211   } else {
0212     throw cms::Exception("UnimplementedFeature") << "unsupported Bounds class";
0213   }
0214
0215   const bool useMatched = false;
0216   int lay =
0217       layerNrConv_.convertLayerNumber(detid.subdetId(),
0218                                       trackerTopo_->layer(detid),
0219                                       useMatched,
0220                                       trackerTopo_->isStereo(detid),
0221                                       trackerTopo_->side(detid) == static_cast<unsigned>(TrackerDetSide::PosEndcap));
0222
0223   mkfit::LayerInfo &layer_info = trk_info.layer_nc(lay);
0224   if (lgc_map) {
0225     (*lgc_map)[lay].reset_current();
0226   }
0227   for (int i = 0; i < 4; ++i) {
0228     Local3DPoint lp1(xy[i][0], xy[i][1], -dz);
0229     Local3DPoint lp2(xy[i][0], xy[i][1], dz);
0230     GlobalPoint gp1 = det->surface().toGlobal(lp1);
0231     GlobalPoint gp2 = det->surface().toGlobal(lp2);
0232     considerPoint(gp1, layer_info);
0233     considerPoint(gp2, layer_info);
0234     if (lgc_map) {
0235       (*lgc_map)[lay].extend_current(gp1.perp2());
0236       (*lgc_map)[lay].extend_current(gp2.perp2());
0237     }
0238   }
0239   if (lgc_map) {
0240     (*lgc_map)[lay].add_current();
0241   }
0242   // Module information
0243   const auto &p = det->position();
0244   auto z = det->rotation().z();
0245   auto x = det->rotation().x();
0246   layer_info.register_module({{p.x(), p.y(), p.z()}, {z.x(), z.y(), z.z()}, {x.x(), x.y(), x.z()}, detid.rawId()});
0247   // Set some layer parameters (repeatedly, would require hard-coding otherwise)
0248   layer_info.set_subdet(detid.subdetId());
0249   layer_info.set_is_pixel(detid.subdetId() <= 2);
0250   layer_info.set_is_stereo(trackerTopo_->isStereo(detid));
0251 }
0252
0253 //==============================================================================
0254
0255 // Ideally these functions would also:
0256 // 0. Setup LayerInfo data (which is now done in auto-generated code).
0257 //    Some data-members are a bit over specific, esp/ bools for CMS sub-detectors.
0258 // 1. Establish short module ids (now done in MkFitGeometry constructor).
0259 // 2. Store module normal and strip direction vectors
0260 // 3. ? Any other information ?
0261 // 4. Extract stereo coverage holes where they exist (TEC, all but last 3 double-layers).
0262 //
0263 // Plugin DumpMkFitGeometry.cc can then be used to export this for stand-alone.
0264 // Would also need to be picked up with tk-ntuple converter (to get module ids as
0265 // they will now be used as indices into module info vectors).
0266 //
0267 // An attempt at export cmsRun config is in python/dumpMkFitGeometry.py
0268
0269 void MkFitGeometryESProducer::addPixBGeometry(mkfit::TrackerInfo &trk_info) {
0270   for (auto &det : trackerGeom_->detsPXB()) {
0271     fillShapeAndPlacement(det, trk_info);
0272   }
0273 }
0274
0275 void MkFitGeometryESProducer::addPixEGeometry(mkfit::TrackerInfo &trk_info) {
0276   for (auto &det : trackerGeom_->detsPXF()) {
0277     fillShapeAndPlacement(det, trk_info);
0278   }
0279 }
0280
0281 void MkFitGeometryESProducer::addTIBGeometry(mkfit::TrackerInfo &trk_info) {
0282   for (auto &det : trackerGeom_->detsTIB()) {
0283     fillShapeAndPlacement(det, trk_info);
0284   }
0285 }
0286
0287 void MkFitGeometryESProducer::addTOBGeometry(mkfit::TrackerInfo &trk_info) {
0288   for (auto &det : trackerGeom_->detsTOB()) {
0289     fillShapeAndPlacement(det, trk_info);
0290   }
0291 }
0292
0293 void MkFitGeometryESProducer::addTIDGeometry(mkfit::TrackerInfo &trk_info) {
0294   for (auto &det : trackerGeom_->detsTID()) {
0295     fillShapeAndPlacement(det, trk_info);
0296   }
0297 }
0298
0299 void MkFitGeometryESProducer::addTECGeometry(mkfit::TrackerInfo &trk_info) {
0300   // For TEC we also need to discover hole in radial extents.
0301   layer_gap_map_t lgc_map;
0302   for (auto &det : trackerGeom_->detsTEC()) {
0303     fillShapeAndPlacement(det, trk_info, &lgc_map);
0304   }
0305   // Now loop over the GapCollectors and see if there is a coverage gap.
0306   std::ostringstream ostr;
0307   ostr << "addTECGeometry() gap report:\n";
0308   GapCollector::Interval itvl;
0309   for (auto &[layer, gcol] : lgc_map) {
0310     gcol.sqrt_elements();
0311     if (gcol.find_gap(itvl, 0.5)) {
0312       ostr << "  layer: " << layer << ", gap: " << itvl.x << " -> " << itvl.y << " width = " << itvl.y - itvl.x << "\n";
0313       ostr << "    all gaps: ";
0314       gcol.print_gaps(ostr);
0315       ostr << "\n";
0316       trk_info.layer_nc(layer).set_r_hole_range(itvl.x, itvl.y);
0317     }
0318   }
0319   edm::LogVerbatim("MkFitGeometryESProducer") << ostr.str();
0320 }
0321
0322 //------------------------------------------------------------------------------
0323 // clang-format off
0324 namespace {
0325   const float phase1QBins[] = {
0326     // PIXB, TIB, TOB
0327     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,
0328     // PIXE+, TID+, TEC+
0329     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,
0330     // PIXE-, TID-, TEC-
0331     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
0332   };
0333 }
0334 // clang-format on
0335 //------------------------------------------------------------------------------
0336
0337 std::unique_ptr<MkFitGeometry> MkFitGeometryESProducer::produce(const TrackerRecoGeometryRecord &iRecord) {
0338   auto trackerInfo = std::make_unique<mkfit::TrackerInfo>();
0339
0340   trackerGeom_ = &iRecord.get(geomToken_);
0341   trackerTopo_ = &iRecord.get(ttopoToken_);
0342
0343   // std::string path = "Geometry/TrackerCommonData/data/";
0344   if (trackerGeom_->isThere(GeomDetEnumerators::P1PXB) || trackerGeom_->isThere(GeomDetEnumerators::P1PXEC)) {
0345     edm::LogInfo("MkFitGeometryESProducer") << "Extracting PhaseI geometry";
0346     trackerInfo->create_layers(18, 27, 27);
0347   } else if (trackerGeom_->isThere(GeomDetEnumerators::P2PXB) || trackerGeom_->isThere(GeomDetEnumerators::P2PXEC) ||
0348              trackerGeom_->isThere(GeomDetEnumerators::P2OTB) || trackerGeom_->isThere(GeomDetEnumerators::P2OTEC)) {
0349     throw cms::Exception("UnimplementedFeature") << "PhaseII geometry extraction";
0350   } else {
0351     throw cms::Exception("UnimplementedFeature") << "unsupported / unknowen geometry version";
0352   }
0353
0354   // Prepare layer boundaries for bounding-box search
0355   for (int i = 0; i < trackerInfo->n_layers(); ++i) {
0356     auto &li = trackerInfo->layer_nc(i);
0357     li.set_limits(
0358         std::numeric_limits<float>::max(), 0, std::numeric_limits<float>::max(), -std::numeric_limits<float>::max());
0359     li.reserve_modules(256);
0360   }
0361   // This is sort of CMS-2017 specific ... but fireworks code uses it for PhaseII as well.
0362   addPixBGeometry(*trackerInfo);
0363   addPixEGeometry(*trackerInfo);
0364   addTIBGeometry(*trackerInfo);
0365   addTIDGeometry(*trackerInfo);
0366   addTOBGeometry(*trackerInfo);
0367   addTECGeometry(*trackerInfo);
0368
0369   // r_in/out kept as squres until here, root them
0370   for (int i = 0; i < trackerInfo->n_layers(); ++i) {
0371     auto &li = trackerInfo->layer_nc(i);
0372     li.set_r_in_out(std::sqrt(li.rin()), std::sqrt(li.rout()));
0373     li.set_propagate_to(li.is_barrel() ? li.r_mean() : li.z_mean());
0374     li.set_q_bin(phase1QBins[i]);
0375     unsigned int maxsid = li.shrink_modules();
0376     // Make sure the short id fits in the 12 bits...
0377     assert(maxsid < 1u << 11);
0378   }
0379
0380   return std::make_unique<MkFitGeometry>(
0381       iRecord.get(geomToken_), iRecord.get(trackerToken_), iRecord.get(ttopoToken_), std::move(trackerInfo));
0382 }
0383
0384 DEFINE_FWK_EVENTSETUP_MODULE(MkFitGeometryESProducer);