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 #include "DataFormats/SiStripDetId/interface/SiStripEnums.h"
0016
0017 // mkFit includes
0018 #include "RecoTracker/MkFit/interface/MkFitGeometry.h"
0019 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0020 #include "RecoTracker/MkFitCore/interface/IterationConfig.h"
0021 #include "RecoTracker/MkFitCMS/interface/LayerNumberConverter.h"
0022 #include "RecoTracker/MkFitCore/interface/Config.h"
0023
0024 #include <list>
0025 #include <vector>
0026 #include <unordered_map>
0027 #include <sstream>
0028
0029 // #define DUMP_MKF_GEO
0030
0031 //------------------------------------------------------------------------------
0032
0033 class MkFitGeometryESProducer : public edm::ESProducer {
0034 public:
0035   MkFitGeometryESProducer(const edm::ParameterSet &iConfig);
0036
0037   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0038
0039   std::unique_ptr<MkFitGeometry> produce(const TrackerRecoGeometryRecord &iRecord);
0040
0041 private:
0042   struct GapCollector {
0043     struct Interval {
0044       float x, y;
0045     };
0046
0047     void reset_current() { m_current = {std::numeric_limits<float>::max(), -std::numeric_limits<float>::max()}; }
0048     void extend_current(float q) {
0049       m_current.x = std::min(m_current.x, q);
0050       m_current.y = std::max(m_current.y, q);
0051     }
0052     void add_current() { add_interval(m_current.x, m_current.y); }
0053
0054     void add_interval(float x, float y);
0055
0056     void sqrt_elements();
0057     bool find_gap(Interval &itvl, float eps);
0058     void print_gaps(std::ostream &ostr);
0059
0060     std::list<Interval> m_coverage;
0061     Interval m_current;
0062   };
0063   using layer_gap_map_t = std::unordered_map<int, GapCollector>;
0064
0065   struct ModuleShape_hash {
0066     std::size_t operator()(const mkfit::ModuleShape &s) const noexcept {
0067       return std::hash<float>{}(s.dx1 + s.dx2 + s.dy + s.dz);
0068     }
0069   };
0070   using module_shape_hmap_t = std::unordered_map<mkfit::ModuleShape, unsigned short, ModuleShape_hash>;
0071   using layer_module_shape_vec_t = std::vector<module_shape_hmap_t>;
0072
0073   struct MatHistBin {
0074     float weight{0}, xi{0}, rl{0};
0075     void add(float w, float x, float r) {
0076       weight += w;
0077       xi += w * x;
0078       rl += w * r;
0079     }
0080   };
0081   using MaterialHistogram = mkfit::rectvec<MatHistBin>;
0082
0083   void considerPoint(const GlobalPoint &gp, mkfit::LayerInfo &lay_info);
0084   void fillShapeAndPlacement(const GeomDet *det,
0085                              mkfit::TrackerInfo &trk_info,
0086                              MaterialHistogram &material_histogram,
0087                              layer_gap_map_t *lgc_map = nullptr);
0088   void addPixBGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram);
0089   void addPixEGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram);
0090   void addTIBGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram);
0091   void addTOBGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram);
0092   void addTIDGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram);
0093   void addTECGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram);
0094
0095   void findRZBox(const GlobalPoint &gp, float &rmin, float &rmax, float &zmin, float &zmax);
0096   void aggregateMaterialInfo(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram);
0097   void fillLayers(mkfit::TrackerInfo &trk_info);
0098
0099   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0100   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> ttopoToken_;
0101   edm::ESGetToken<GeometricSearchTracker, TrackerRecoGeometryRecord> trackerToken_;
0102
0103   const TrackerTopology *trackerTopo_ = nullptr;
0104   const TrackerGeometry *trackerGeom_ = nullptr;
0105   mkfit::LayerNumberConverter layerNrConv_ = {mkfit::TkLayout::phase1};
0106   layer_module_shape_vec_t layerModuleShapeVec_;
0107 };
0108
0109 MkFitGeometryESProducer::MkFitGeometryESProducer(const edm::ParameterSet &iConfig) {
0110   auto cc = setWhatProduced(this);
0111   geomToken_ = cc.consumes();
0112   ttopoToken_ = cc.consumes();
0113   trackerToken_ = cc.consumes();
0114 }
0115
0116 void MkFitGeometryESProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0117   edm::ParameterSetDescription desc;
0118   descriptions.addWithDefaultLabel(desc);
0119 }
0120
0121 //------------------------------------------------------------------------------
0122
0123 void MkFitGeometryESProducer::GapCollector::add_interval(float x, float y) {
0124   if (x > y)
0125     std::swap(x, y);
0126   bool handled = false;
0127   for (auto i = m_coverage.begin(); i != m_coverage.end(); ++i) {
0128     if (y < i->x) {  // fully on 'left'
0129       m_coverage.insert(i, {x, y});
0130       handled = true;
0131       break;
0132     } else if (x > i->y) {  // fully on 'right'
0133       continue;
0134     } else if (x < i->x) {  // sticking out on 'left'
0135       i->x = x;
0136       handled = true;
0137       break;
0138     } else if (y > i->y) {  // sticking out on 'right'
0139       i->y = y;
0140       // check for overlap with the next interval, potentially merge
0141       auto j = i;
0142       ++j;
0143       if (j != m_coverage.end() && i->y >= j->x) {
0144         i->y = j->y;
0145         m_coverage.erase(j);
0146       }
0147       handled = true;
0148       break;
0149     } else {  // contained in current interval
0150       handled = true;
0151       break;
0152     }
0153   }
0154   if (!handled) {
0155     m_coverage.push_back({x, y});
0156   }
0157 }
0158
0159 void MkFitGeometryESProducer::GapCollector::sqrt_elements() {
0160   for (auto &itvl : m_coverage) {
0161     itvl.x = std::sqrt(itvl.x);
0162     itvl.y = std::sqrt(itvl.y);
0163   }
0164 }
0165
0166 bool MkFitGeometryESProducer::GapCollector::find_gap(Interval &itvl, float eps) {
0167   auto i = m_coverage.begin();
0168   while (i != m_coverage.end()) {
0169     auto j = i;
0170     ++j;
0171     if (j != m_coverage.end()) {
0172       if (j->x - i->y > eps) {
0173         itvl = {i->y, j->x};
0174         return true;
0175       }
0176       i = j;
0177     } else {
0178       break;
0179     }
0180   }
0181   return false;
0182 }
0183
0184 void MkFitGeometryESProducer::GapCollector::print_gaps(std::ostream &ostr) {
0185   auto i = m_coverage.begin();
0186   while (i != m_coverage.end()) {
0187     auto j = i;
0188     ++j;
0189     if (j != m_coverage.end()) {
0190       ostr << "(" << i->y << ", " << j->x << ")->" << j->x - i->y;
0191       i = j;
0192     } else {
0193       break;
0194     }
0195   }
0196 }
0197
0198 //------------------------------------------------------------------------------
0199
0200 void MkFitGeometryESProducer::considerPoint(const GlobalPoint &gp, mkfit::LayerInfo &li) {
0201   // Use radius squared during bounding-region search.
0202   float r = gp.perp2(), z = gp.z();
0203   li.extend_limits(r, z);
0204 }
0205
0206 void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det,
0207                                                     mkfit::TrackerInfo &trk_info,
0208                                                     MaterialHistogram &material_histogram,
0209                                                     layer_gap_map_t *lgc_map) {
0210   const DetId detid = det->geographicalId();
0211
0212   bool doubleSide = false;  //double modules have double material
0213   if (detid.subdetId() == SiStripSubdetector::TIB)
0214     doubleSide = trackerTopo_->tibIsDoubleSide(detid);
0215   else if (detid.subdetId() == SiStripSubdetector::TID)
0216     doubleSide = trackerTopo_->tidIsDoubleSide(detid);
0217   else if (detid.subdetId() == SiStripSubdetector::TOB)
0218     doubleSide = trackerTopo_->tobIsDoubleSide(detid);
0219   else if (detid.subdetId() == SiStripSubdetector::TEC)
0220     doubleSide = trackerTopo_->tecIsDoubleSide(detid);
0221
0222   float xy[4][2];
0223   float dz;
0224   const Bounds *b = &((det->surface()).bounds());
0225   mkfit::ModuleShape ms;
0226
0227   if (const TrapezoidalPlaneBounds *b2 = dynamic_cast<const TrapezoidalPlaneBounds *>(b)) {
0228     // See sec. "TrapezoidalPlaneBounds parameters" in doc/reco-geom-notes.txt
0229     std::array<const float, 4> const &par = b2->parameters();
0230     xy[0][0] = -par[0];
0231     xy[0][1] = -par[3];
0232     xy[1][0] = -par[1];
0233     xy[1][1] = par[3];
0234     xy[2][0] = par[1];
0235     xy[2][1] = par[3];
0236     xy[3][0] = par[0];
0237     xy[3][1] = -par[3];
0238     dz = par[2];
0239     ms.round_assign(par[0], par[1], par[3], par[2]);
0240
0241 #ifdef DUMP_MKF_GEO
0242     printf("TRAP 0x%x %f %f %f %f  ", detid.rawId(), par[0], par[1], par[2], par[3]);
0243 #endif
0244   } else if (const RectangularPlaneBounds *b2 = dynamic_cast<const RectangularPlaneBounds *>(b)) {
0245     // Rectangular
0246     float dx = b2->width() * 0.5;   // half width
0247     float dy = b2->length() * 0.5;  // half length
0248     xy[0][0] = -dx;
0249     xy[0][1] = -dy;
0250     xy[1][0] = -dx;
0251     xy[1][1] = dy;
0252     xy[2][0] = dx;
0253     xy[2][1] = dy;
0254     xy[3][0] = dx;
0255     xy[3][1] = -dy;
0256     dz = b2->thickness() * 0.5;  // half thickness
0257     ms.round_assign(dx, 0.0f, dy, dz);
0258
0259 #ifdef DUMP_MKF_GEO
0260     printf("RECT 0x%x %f %f %f           ", detid.rawId(), dx, dy, dz);
0261 #endif
0262   } else {
0263     throw cms::Exception("UnimplementedFeature") << "unsupported Bounds class";
0264   }
0265
0266   const bool useMatched = false;
0267   int lay =
0268       layerNrConv_.convertLayerNumber(detid.subdetId(),
0269                                       trackerTopo_->layer(detid),
0270                                       useMatched,
0271                                       trackerTopo_->isStereo(detid),
0272                                       trackerTopo_->side(detid) == static_cast<unsigned>(TrackerDetSide::PosEndcap));
0273
0274   unsigned short shape_id = 9999;
0275   if (!doubleSide) {
0276     module_shape_hmap_t &bhm = layerModuleShapeVec_[lay];
0277     auto bhmi = bhm.find(ms);
0278     if (bhmi == bhm.end()) {
0279       bhmi = bhm.insert({ms, (unsigned short)bhm.size()}).first;
0280     }
0281     shape_id = bhmi->second;
0282   }
0283
0284 #ifdef DUMP_MKF_GEO
0285   printf("  subdet=%d layer=%d side=%d is_stereo=%d is_double_side=%d --> mkflayer=%d; unique shape id=%hu\n",
0286          detid.subdetId(),
0287          trackerTopo_->layer(detid),
0288          trackerTopo_->side(detid),
0289          trackerTopo_->isStereo(detid),
0290          doubleSide,
0291          lay,
0292          shape_id);
0293 #endif
0294
0295   mkfit::LayerInfo &layer_info = trk_info.layer_nc(lay);
0296   if (lgc_map) {
0297     (*lgc_map)[lay].reset_current();
0298   }
0299   float zbox_min = 1000, zbox_max = 0, rbox_min = 1000, rbox_max = 0;
0300   for (int i = 0; i < 4; ++i) {
0301     Local3DPoint lp1(xy[i][0], xy[i][1], -dz);
0302     Local3DPoint lp2(xy[i][0], xy[i][1], dz);
0303     GlobalPoint gp1 = det->surface().toGlobal(lp1);
0304     GlobalPoint gp2 = det->surface().toGlobal(lp2);
0305     considerPoint(gp1, layer_info);
0306     considerPoint(gp2, layer_info);
0307     findRZBox(gp1, rbox_min, rbox_max, zbox_min, zbox_max);
0308     findRZBox(gp2, rbox_min, rbox_max, zbox_min, zbox_max);
0309     if (lgc_map) {
0310       (*lgc_map)[lay].extend_current(gp1.perp2());
0311       (*lgc_map)[lay].extend_current(gp2.perp2());
0312     }
0313   }
0314   if (lgc_map) {
0315     (*lgc_map)[lay].add_current();
0316   }
0317
0318   // Double-sided module (join of two modules) information is not used in mkFit and
0319   // also not needed for the material calculation.
0320   // NOTE: This check should actually be performed even before the bounding box calculation
0321   // as double-side module bounding box is a box-shaped union of bounding shapes of the
0322   // constituent modules -- thus extending the layer bounding box unnecessarily.
0323   // The 'if' is here as mkFit layer hit/miss logic has been tuned with this extended
0324   // bounding boxes. To be moved higher up once this detection is improved through usage of
0325   // the thin-thick layer concept (expected towards the end of 2024).
0326   if (doubleSide)
0327     return;
0328
0329   // Module information
0330   const auto &p = det->position();
0331   auto z = det->rotation().z();
0332   auto x = det->rotation().x();
0333   layer_info.register_module(
0334       {{p.x(), p.y(), p.z()}, {z.x(), z.y(), z.z()}, {x.x(), x.y(), x.z()}, detid.rawId(), shape_id});
0335   // Set some layer parameters (repeatedly, would require hard-coding otherwise)
0336   layer_info.set_subdet(detid.subdetId());
0337   layer_info.set_is_pixel(detid.subdetId() <= 2);
0338   layer_info.set_is_stereo(trackerTopo_->isStereo(detid));
0339   if (layerNrConv_.isPhase2() && !layer_info.is_pixel())
0340     layer_info.set_has_charge(false);
0341
0342   // Fill material
0343   {
0344     // module material
0345     const float bbxi = det->surface().mediumProperties().xi();
0346     const float radL = det->surface().mediumProperties().radLen();
0347     // loop over bins to fill histogram with bbxi, radL and their weight, which the overlap surface in r-z with the cmsquare of a bin
0348     const float iBin = trk_info.mat_range_z() / trk_info.mat_nbins_z();
0349     const float jBin = trk_info.mat_range_r() / trk_info.mat_nbins_r();
0350     for (int i = std::floor(zbox_min / iBin); i < std::ceil(zbox_max / iBin); i++) {
0351       for (int j = std::floor(rbox_min / jBin); j < std::ceil(rbox_max / jBin); j++) {
0352         const float iF = i * iBin;
0353         const float jF = j * jBin;
0354         float overlap = std::max(0.f, std::min(jF + jBin, rbox_max) - std::max(jF, rbox_min)) *
0355                         std::max(0.f, std::min(iF + iBin, zbox_max) - std::max(iF, zbox_min));
0356         if (overlap > 0)
0357           material_histogram(i, j).add(overlap, bbxi, radL);
0358       }
0359     }
0360   }
0361 }
0362
0363 //==============================================================================
0364
0365 // These functions do the following:
0366 // 0. Detect bounding cylinder of each layer.
0367 // 1. Setup LayerInfo data.
0368 // 2. Establish short module ids.
0369 // 3. Store module normal and strip direction vectors.
0370 // 4. Extract stereo coverage holes where they exist (TEC, all but last 3 double-layers).
0371 //
0372 // See python/dumpMkFitGeometry.py and dumpMkFitGeometryPhase2.py
0373
0374 void MkFitGeometryESProducer::addPixBGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0375 #ifdef DUMP_MKF_GEO
0376   printf("\n*** addPixBGeometry\n\n");
0377 #endif
0378   for (auto &det : trackerGeom_->detsPXB()) {
0379     fillShapeAndPlacement(det, trk_info, material_histogram);
0380   }
0381 }
0382
0383 void MkFitGeometryESProducer::addPixEGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0384 #ifdef DUMP_MKF_GEO
0385   printf("\n*** addPixEGeometry\n\n");
0386 #endif
0387   for (auto &det : trackerGeom_->detsPXF()) {
0388     fillShapeAndPlacement(det, trk_info, material_histogram);
0389   }
0390 }
0391
0392 void MkFitGeometryESProducer::addTIBGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0393 #ifdef DUMP_MKF_GEO
0394   printf("\n*** addTIBGeometry\n\n");
0395 #endif
0396   for (auto &det : trackerGeom_->detsTIB()) {
0397     fillShapeAndPlacement(det, trk_info, material_histogram);
0398   }
0399 }
0400
0401 void MkFitGeometryESProducer::addTOBGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0402 #ifdef DUMP_MKF_GEO
0403   printf("\n*** addTOBGeometry\n\n");
0404 #endif
0405   for (auto &det : trackerGeom_->detsTOB()) {
0406     fillShapeAndPlacement(det, trk_info, material_histogram);
0407   }
0408 }
0409
0410 void MkFitGeometryESProducer::addTIDGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0411 #ifdef DUMP_MKF_GEO
0412   printf("\n*** addTIDGeometry\n\n");
0413 #endif
0414   for (auto &det : trackerGeom_->detsTID()) {
0415     fillShapeAndPlacement(det, trk_info, material_histogram);
0416   }
0417 }
0418
0419 void MkFitGeometryESProducer::addTECGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0420 #ifdef DUMP_MKF_GEO
0421   printf("\n*** addTECGeometry\n\n");
0422 #endif
0423   // For TEC we also need to discover hole in radial extents.
0424   layer_gap_map_t lgc_map;
0425   for (auto &det : trackerGeom_->detsTEC()) {
0426     fillShapeAndPlacement(det, trk_info, material_histogram, &lgc_map);
0427   }
0428   // Now loop over the GapCollectors and see if there is a coverage gap.
0429   std::ostringstream ostr;
0430   ostr << "addTECGeometry() gap report:\n";
0431   GapCollector::Interval itvl;
0432   for (auto &[layer, gcol] : lgc_map) {
0433     gcol.sqrt_elements();
0434     if (gcol.find_gap(itvl, 0.5)) {
0435       ostr << "  layer: " << layer << ", gap: " << itvl.x << " -> " << itvl.y << " width = " << itvl.y - itvl.x << "\n";
0436       ostr << "    all gaps: ";
0437       gcol.print_gaps(ostr);
0438       ostr << "\n";
0439       trk_info.layer_nc(layer).set_r_hole_range(itvl.x, itvl.y);
0440     }
0441   }
0442   edm::LogVerbatim("MkFitGeometryESProducer") << ostr.str();
0443 }
0444
0445 void MkFitGeometryESProducer::findRZBox(const GlobalPoint &gp, float &rmin, float &rmax, float &zmin, float &zmax) {
0446   float r = gp.perp(), z = std::abs(gp.z());
0447   rmax = std::max(r, rmax);
0448   rmin = std::min(r, rmin);
0449   zmax = std::max(z, zmax);
0450   zmin = std::min(z, zmin);
0451 }
0452
0453 void MkFitGeometryESProducer::aggregateMaterialInfo(mkfit::TrackerInfo &trk_info,
0454                                                     MaterialHistogram &material_histogram) {
0455   //from histogram (vector of tuples) to grid
0456   for (int i = 0; i < trk_info.mat_nbins_z(); i++) {
0457     for (int j = 0; j < trk_info.mat_nbins_r(); j++) {
0458       const MatHistBin &mhb = material_histogram(i, j);
0459       if (mhb.weight > 0) {
0460         trk_info.material_bbxi(i, j) = mhb.xi / mhb.weight;
0461         trk_info.material_radl(i, j) = mhb.rl / mhb.weight;
0462       }
0463     }
0464   }
0465 }
0466
0467 void MkFitGeometryESProducer::fillLayers(mkfit::TrackerInfo &trk_info) {
0468   mkfit::rectvec<int> rneighbor_map(trk_info.mat_nbins_z(), trk_info.mat_nbins_r());
0469   mkfit::rectvec<int> zneighbor_map(trk_info.mat_nbins_z(), trk_info.mat_nbins_r());
0470
0471   for (int im = 0; im < trk_info.n_layers(); ++im) {
0472     const mkfit::LayerInfo &li = trk_info.layer(im);
0473     if (!li.is_barrel() && li.zmax() < 0)
0474       continue;  // neg endcap covered by pos
0475     int rin, rout, zmin, zmax;
0476     rin = trk_info.mat_bin_r(li.rin());
0477     rout = trk_info.mat_bin_r(li.rout()) + 1;
0478     if (li.is_barrel()) {
0479       zmin = 0;
0480       zmax = trk_info.mat_bin_z(std::max(std::abs(li.zmax()), std::abs(li.zmin()))) + 1;
0481     } else {
0482       zmin = trk_info.mat_bin_z(li.zmin());
0483       zmax = trk_info.mat_bin_z(li.zmax()) + 1;
0484     }
0485     for (int i = zmin; i < zmax; i++) {
0486       for (int j = rin; j < rout; j++) {
0487         if (trk_info.material_bbxi(i, j) == 0) {
0488           float distancesqmin = 100000;
0489           for (int i2 = zmin; i2 < zmax; i2++) {
0490             for (int j2 = rin; j2 < rout; j2++) {
0491               if (j == j2 && i == i2)
0492                 continue;
0493               auto mydistsq = (i - i2) * (i - i2) + (j - j2) * (j - j2);
0494               if (mydistsq < distancesqmin && trk_info.material_radl(i2, j2) > 0) {
0495                 distancesqmin = mydistsq;
0496                 zneighbor_map(i, j) = i2;
0497                 rneighbor_map(i, j) = j2;
0498               }
0499             }
0500           }  // can work on speedup here
0501         }
0502       }
0503     }
0504     for (int i = zmin; i < zmax; i++) {
0505       for (int j = rin; j < rout; j++) {
0506         if (trk_info.material_bbxi(i, j) == 0) {
0507           int iN = zneighbor_map(i, j);
0508           int jN = rneighbor_map(i, j);
0509           trk_info.material_bbxi(i, j) = trk_info.material_bbxi(iN, jN);
0510           trk_info.material_radl(i, j) = trk_info.material_radl(iN, jN);
0511         }
0512       }
0513     }
0514   }  //module loop
0515 }
0516
0517 //------------------------------------------------------------------------------
0518 // clang-format off
0519 namespace {
0520   const float phase1QBins[] = {
0521     // PIXB, TIB, TOB
0522     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,
0523     // PIXE+, TID+, TEC+
0524     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,
0525     // PIXE-, TID-, TEC-
0526     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
0527   };
0528   const float phase2QBins[] = {
0529     // TODO: Review these numbers.
0530     // PIXB, TOB
0531     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,
0532     // PIXE+, TEC+
0533     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,
0534     // PIXE-, TEC-
0535     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
0536   };
0537 }
0538 // clang-format on
0539 //------------------------------------------------------------------------------
0540
0541 std::unique_ptr<MkFitGeometry> MkFitGeometryESProducer::produce(const TrackerRecoGeometryRecord &iRecord) {
0542   auto trackerInfo = std::make_unique<mkfit::TrackerInfo>();
0543
0544   trackerGeom_ = &iRecord.get(geomToken_);
0545   trackerTopo_ = &iRecord.get(ttopoToken_);
0546
0547   const float *qBinDefaults = nullptr;
0548
0549   // std::string path = "Geometry/TrackerCommonData/data/";
0550   if (trackerGeom_->isThere(GeomDetEnumerators::P1PXB) || trackerGeom_->isThere(GeomDetEnumerators::P1PXEC)) {
0551     edm::LogInfo("MkFitGeometryESProducer") << "Extracting PhaseI geometry";
0552     trackerInfo->create_layers(18, 27, 27);
0553     qBinDefaults = phase1QBins;
0554     trackerInfo->create_material(300, 300.0f, 120, 120.0f);
0555   } else if (trackerGeom_->isThere(GeomDetEnumerators::P2PXB) || trackerGeom_->isThere(GeomDetEnumerators::P2PXEC) ||
0556              trackerGeom_->isThere(GeomDetEnumerators::P2OTB) || trackerGeom_->isThere(GeomDetEnumerators::P2OTEC)) {
0557     edm::LogInfo("MkFitGeometryESProducer") << "Extracting PhaseII geometry";
0558 #if !defined(MKFIT_PHASE2CUSTOMFLAGS)
0559     // In Phase2, by default use prop-to-plane and pT-dependent MS.
0560     // Option is kept to use custom flags, for R&D/test purposes.
0561     using namespace mkfit;
0562     Config::usePropToPlane = true;
0563     Config::usePtMultScat = true;
0564 #endif
0565     layerNrConv_.reset(mkfit::TkLayout::phase2);
0566     trackerInfo->create_layers(16, 22, 22);
0567     qBinDefaults = phase2QBins;
0568     trackerInfo->create_material(300, 300.0f, 120, 120.0f);
0569   } else {
0570     throw cms::Exception("UnimplementedFeature") << "unsupported / unknowen geometry version";
0571   }
0572
0573   // Prepare layer boundaries for bounding-rz-box search
0574   for (int i = 0; i < trackerInfo->n_layers(); ++i) {
0575     auto &li = trackerInfo->layer_nc(i);
0576     li.set_limits(
0577         std::numeric_limits<float>::max(), 0, std::numeric_limits<float>::max(), -std::numeric_limits<float>::max());
0578     li.reserve_modules(256);
0579   }
0580   // Resize the module bounds collector vector
0581   layerModuleShapeVec_.resize(trackerInfo->n_layers());
0582
0583   MaterialHistogram material_histogram(trackerInfo->mat_nbins_z(), trackerInfo->mat_nbins_r());
0584
0585   // This works for both Phase1 and Phase2.
0586   // Phase2 TrackerGeometry returns empty det-vectors for TIB and TEC.
0587   addPixBGeometry(*trackerInfo, material_histogram);
0588   addPixEGeometry(*trackerInfo, material_histogram);
0589   addTIBGeometry(*trackerInfo, material_histogram);
0590   addTIDGeometry(*trackerInfo, material_histogram);
0591   addTOBGeometry(*trackerInfo, material_histogram);
0592   addTECGeometry(*trackerInfo, material_histogram);
0593
0594   // r_in/out kept as squares until here, root them; fill ModuleShape vectors of layers
0595   unsigned int n_mod = 0;
0596   for (int i = 0; i < trackerInfo->n_layers(); ++i) {
0597     auto &li = trackerInfo->layer_nc(i);
0598     li.set_r_in_out(std::sqrt(li.rin()), std::sqrt(li.rout()));
0599     li.set_propagate_to(li.is_barrel() ? li.r_mean() : li.z_mean());
0600     li.set_q_bin(qBinDefaults[i]);
0601     unsigned int maxsid = li.shrink_modules();
0602
0603     n_mod += maxsid;
0604
0605     // Make sure the short id fits in the 14 bits...
0606     assert(maxsid < 1u << 13);
0607     assert(n_mod > 0);
0608
0609     // Fill ModuleShape vectors
0610     int n_shapes = layerModuleShapeVec_[i].size();
0611     li.resize_shapes(n_shapes);
0612     for (auto &[shape, id] : layerModuleShapeVec_[i]) {
0613       li.register_shape(shape, id);
0614     }
0615   }
0616
0617   // Material grid
0618   aggregateMaterialInfo(*trackerInfo, material_histogram);
0619   fillLayers(*trackerInfo);
0620
0621   // Propagation configuration
0622   {
0623     using namespace mkfit;
0624     PropagationConfig &pconf = trackerInfo->prop_config_nc();
0625     pconf.backward_fit_to_pca = false;
0626     pconf.finding_requires_propagation_to_hit_pos = true;
0627     pconf.finding_inter_layer_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
0628     if (Config::usePropToPlane)
0629       pconf.finding_intra_layer_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
0630     else
0631       pconf.finding_intra_layer_pflags = PropagationFlags(PF_none);
0632     pconf.backward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
0633     pconf.forward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
0634     pconf.seed_fit_pflags = PropagationFlags(PF_none);
0635     pconf.pca_prop_pflags = PropagationFlags(PF_none);
0636     pconf.apply_tracker_info(trackerInfo.get());
0637   }
0638
0639 #ifdef DUMP_MKF_GEO
0640   printf("Total number of modules %u, 14-bits fit up to %u modules\n", n_mod, 1u << 13);
0641 #endif
0642
0643   return std::make_unique<MkFitGeometry>(iRecord.get(geomToken_),
0644                                          iRecord.get(trackerToken_),
0645                                          iRecord.get(ttopoToken_),
0646                                          std::move(trackerInfo),
0647                                          layerNrConv_);
0648 }
0649
0650 DEFINE_FWK_EVENTSETUP_MODULE(MkFitGeometryESProducer);