Back to home page

Project CMSSW displayed by LXR
 
 

    

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 // 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 // #define DUMP_MKF_GEO
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) {  // fully on 'left'
0098       m_coverage.insert(i, {x, y});
0099       handled = true;
0100       break;
0101     } else if (x > i->y) {  // fully on 'right'
0102       continue;
0103     } else if (x < i->x) {  // sticking out on 'left'
0104       i->x = x;
0105       handled = true;
0106       break;
0107     } else if (y > i->y) {  // sticking out on 'right'
0108       i->y = y;
0109       // check for overlap with the next interval, potentially merge
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 {  // contained in current interval
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   // Use radius squared during bounding-region search.
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     // See sec. "TrapezoidalPlaneBounds parameters" in doc/reco-geom-notes.txt
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     // Rectangular
0202     float dx = b2->width() * 0.5;   // half width
0203     float dy = b2->length() * 0.5;  // half length
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;  // half thickness
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   // Module information
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   // Set some layer parameters (repeatedly, would require hard-coding otherwise)
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 // These functions do the following:
0270 // 0. Detect bounding cylinder of each layer.
0271 // 1. Setup LayerInfo data.
0272 // 2. Establish short module ids.
0273 // 3. Store module normal and strip direction vectors.
0274 // 4. Extract stereo coverage holes where they exist (TEC, all but last 3 double-layers).
0275 //
0276 // See python/dumpMkFitGeometry.py and dumpMkFitGeometryPhase2.py
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   // For TEC we also need to discover hole in radial extents.
0328   layer_gap_map_t lgc_map;
0329   for (auto &det : trackerGeom_->detsTEC()) {
0330     fillShapeAndPlacement(det, trk_info, &lgc_map);
0331   }
0332   // Now loop over the GapCollectors and see if there is a coverage gap.
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 // clang-format off
0351 namespace {
0352   const float phase1QBins[] = {
0353     // PIXB, TIB, TOB
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     // PIXE+, TID+, TEC+
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     // PIXE-, TID-, TEC-
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     // TODO: Review these numbers.
0362     // PIXB, TOB
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     // PIXE+, TEC+
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     // PIXE-, TEC-
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 // clang-format on
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   // std::string path = "Geometry/TrackerCommonData/data/";
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   // Prepare layer boundaries for bounding-box search
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   // This is sort of CMS-phase1 specific ... but fireworks code uses it for PhaseII as well.
0404   addPixBGeometry(*trackerInfo);
0405   addPixEGeometry(*trackerInfo);
0406   addTIBGeometry(*trackerInfo);
0407   addTIDGeometry(*trackerInfo);
0408   addTOBGeometry(*trackerInfo);
0409   addTECGeometry(*trackerInfo);
0410 
0411   // r_in/out kept as squares until here, root them
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     // Make sure the short id fits in the 14 bits...
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!