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