Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:28:13

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   DetId detid = det->geographicalId();
0199 
0200   float xy[4][2];
0201   float half_length, dz;
0202   const Bounds *b = &((det->surface()).bounds());
0203 
0204   if (const TrapezoidalPlaneBounds *b2 = dynamic_cast<const TrapezoidalPlaneBounds *>(b)) {
0205     // See sec. "TrapezoidalPlaneBounds parameters" in doc/reco-geom-notes.txt
0206     std::array<const float, 4> const &par = b2->parameters();
0207     xy[0][0] = -par[0];
0208     xy[0][1] = -par[3];
0209     xy[1][0] = -par[1];
0210     xy[1][1] = par[3];
0211     xy[2][0] = par[1];
0212     xy[2][1] = par[3];
0213     xy[3][0] = par[0];
0214     xy[3][1] = -par[3];
0215     half_length = par[3];
0216     dz = par[2];
0217 
0218 #ifdef DUMP_MKF_GEO
0219     printf("TRAP 0x%x %f %f %f %f  ", detid.rawId(), par[0], par[1], par[2], par[3]);
0220 #endif
0221   } else if (const RectangularPlaneBounds *b2 = dynamic_cast<const RectangularPlaneBounds *>(b)) {
0222     // Rectangular
0223     float dx = b2->width() * 0.5;   // half width
0224     float dy = b2->length() * 0.5;  // half length
0225     xy[0][0] = -dx;
0226     xy[0][1] = -dy;
0227     xy[1][0] = -dx;
0228     xy[1][1] = dy;
0229     xy[2][0] = dx;
0230     xy[2][1] = dy;
0231     xy[3][0] = dx;
0232     xy[3][1] = -dy;
0233     half_length = dy;
0234     dz = b2->thickness() * 0.5;  // half thickness
0235 
0236 #ifdef DUMP_MKF_GEO
0237     printf("RECT 0x%x %f %f %f  ", detid.rawId(), dx, dy, dz);
0238 #endif
0239   } else {
0240     throw cms::Exception("UnimplementedFeature") << "unsupported Bounds class";
0241   }
0242 
0243   const bool useMatched = false;
0244   int lay =
0245       layerNrConv_.convertLayerNumber(detid.subdetId(),
0246                                       trackerTopo_->layer(detid),
0247                                       useMatched,
0248                                       trackerTopo_->isStereo(detid),
0249                                       trackerTopo_->side(detid) == static_cast<unsigned>(TrackerDetSide::PosEndcap));
0250 #ifdef DUMP_MKF_GEO
0251   printf("  subdet=%d layer=%d side=%d is_stereo=%d --> mkflayer=%d\n",
0252          detid.subdetId(),
0253          trackerTopo_->layer(detid),
0254          trackerTopo_->side(detid),
0255          trackerTopo_->isStereo(detid),
0256          lay);
0257 #endif
0258 
0259   mkfit::LayerInfo &layer_info = trk_info.layer_nc(lay);
0260   if (lgc_map) {
0261     (*lgc_map)[lay].reset_current();
0262   }
0263   float zbox_min = 1000, zbox_max = 0, rbox_min = 1000, rbox_max = 0;
0264   for (int i = 0; i < 4; ++i) {
0265     Local3DPoint lp1(xy[i][0], xy[i][1], -dz);
0266     Local3DPoint lp2(xy[i][0], xy[i][1], dz);
0267     GlobalPoint gp1 = det->surface().toGlobal(lp1);
0268     GlobalPoint gp2 = det->surface().toGlobal(lp2);
0269     considerPoint(gp1, layer_info);
0270     considerPoint(gp2, layer_info);
0271     findRZBox(gp1, rbox_min, rbox_max, zbox_min, zbox_max);
0272     findRZBox(gp2, rbox_min, rbox_max, zbox_min, zbox_max);
0273     if (lgc_map) {
0274       (*lgc_map)[lay].extend_current(gp1.perp2());
0275       (*lgc_map)[lay].extend_current(gp2.perp2());
0276     }
0277   }
0278   if (lgc_map) {
0279     (*lgc_map)[lay].add_current();
0280   }
0281   // Module information
0282   const auto &p = det->position();
0283   auto z = det->rotation().z();
0284   auto x = det->rotation().x();
0285   layer_info.register_module(
0286       {{p.x(), p.y(), p.z()}, {z.x(), z.y(), z.z()}, {x.x(), x.y(), x.z()}, half_length, detid.rawId()});
0287   // Set some layer parameters (repeatedly, would require hard-coding otherwise)
0288   layer_info.set_subdet(detid.subdetId());
0289   layer_info.set_is_pixel(detid.subdetId() <= 2);
0290   layer_info.set_is_stereo(trackerTopo_->isStereo(detid));
0291 
0292   bool doubleSide = false;  //double modules have double material
0293   if (detid.subdetId() == SiStripSubdetector::TIB)
0294     doubleSide = trackerTopo_->tibIsDoubleSide(detid);
0295   else if (detid.subdetId() == SiStripSubdetector::TID)
0296     doubleSide = trackerTopo_->tidIsDoubleSide(detid);
0297   else if (detid.subdetId() == SiStripSubdetector::TOB)
0298     doubleSide = trackerTopo_->tobIsDoubleSide(detid);
0299   else if (detid.subdetId() == SiStripSubdetector::TEC)
0300     doubleSide = trackerTopo_->tecIsDoubleSide(detid);
0301 
0302   if (!doubleSide)  //fill material
0303   {
0304     //module material
0305     const float bbxi = det->surface().mediumProperties().xi();
0306     const float radL = det->surface().mediumProperties().radLen();
0307     //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
0308     const float iBin = trk_info.mat_range_z() / trk_info.mat_nbins_z();
0309     const float jBin = trk_info.mat_range_r() / trk_info.mat_nbins_r();
0310     for (int i = std::floor(zbox_min / iBin); i < std::ceil(zbox_max / iBin); i++) {
0311       for (int j = std::floor(rbox_min / jBin); j < std::ceil(rbox_max / jBin); j++) {
0312         const float iF = i * iBin;
0313         const float jF = j * jBin;
0314         float overlap = std::max(0.f, std::min(jF + jBin, rbox_max) - std::max(jF, rbox_min)) *
0315                         std::max(0.f, std::min(iF + iBin, zbox_max) - std::max(iF, zbox_min));
0316         if (overlap > 0)
0317           material_histogram(i, j).add(overlap, bbxi, radL);
0318       }
0319     }
0320   }
0321 }
0322 
0323 //==============================================================================
0324 
0325 // These functions do the following:
0326 // 0. Detect bounding cylinder of each layer.
0327 // 1. Setup LayerInfo data.
0328 // 2. Establish short module ids.
0329 // 3. Store module normal and strip direction vectors.
0330 // 4. Extract stereo coverage holes where they exist (TEC, all but last 3 double-layers).
0331 //
0332 // See python/dumpMkFitGeometry.py and dumpMkFitGeometryPhase2.py
0333 
0334 void MkFitGeometryESProducer::addPixBGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0335 #ifdef DUMP_MKF_GEO
0336   printf("\n*** addPixBGeometry\n\n");
0337 #endif
0338   for (auto &det : trackerGeom_->detsPXB()) {
0339     fillShapeAndPlacement(det, trk_info, material_histogram);
0340   }
0341 }
0342 
0343 void MkFitGeometryESProducer::addPixEGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0344 #ifdef DUMP_MKF_GEO
0345   printf("\n*** addPixEGeometry\n\n");
0346 #endif
0347   for (auto &det : trackerGeom_->detsPXF()) {
0348     fillShapeAndPlacement(det, trk_info, material_histogram);
0349   }
0350 }
0351 
0352 void MkFitGeometryESProducer::addTIBGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0353 #ifdef DUMP_MKF_GEO
0354   printf("\n*** addTIBGeometry\n\n");
0355 #endif
0356   for (auto &det : trackerGeom_->detsTIB()) {
0357     fillShapeAndPlacement(det, trk_info, material_histogram);
0358   }
0359 }
0360 
0361 void MkFitGeometryESProducer::addTOBGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0362 #ifdef DUMP_MKF_GEO
0363   printf("\n*** addTOBGeometry\n\n");
0364 #endif
0365   for (auto &det : trackerGeom_->detsTOB()) {
0366     fillShapeAndPlacement(det, trk_info, material_histogram);
0367   }
0368 }
0369 
0370 void MkFitGeometryESProducer::addTIDGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0371 #ifdef DUMP_MKF_GEO
0372   printf("\n*** addTIDGeometry\n\n");
0373 #endif
0374   for (auto &det : trackerGeom_->detsTID()) {
0375     fillShapeAndPlacement(det, trk_info, material_histogram);
0376   }
0377 }
0378 
0379 void MkFitGeometryESProducer::addTECGeometry(mkfit::TrackerInfo &trk_info, MaterialHistogram &material_histogram) {
0380 #ifdef DUMP_MKF_GEO
0381   printf("\n*** addTECGeometry\n\n");
0382 #endif
0383   // For TEC we also need to discover hole in radial extents.
0384   layer_gap_map_t lgc_map;
0385   for (auto &det : trackerGeom_->detsTEC()) {
0386     fillShapeAndPlacement(det, trk_info, material_histogram, &lgc_map);
0387   }
0388   // Now loop over the GapCollectors and see if there is a coverage gap.
0389   std::ostringstream ostr;
0390   ostr << "addTECGeometry() gap report:\n";
0391   GapCollector::Interval itvl;
0392   for (auto &[layer, gcol] : lgc_map) {
0393     gcol.sqrt_elements();
0394     if (gcol.find_gap(itvl, 0.5)) {
0395       ostr << "  layer: " << layer << ", gap: " << itvl.x << " -> " << itvl.y << " width = " << itvl.y - itvl.x << "\n";
0396       ostr << "    all gaps: ";
0397       gcol.print_gaps(ostr);
0398       ostr << "\n";
0399       trk_info.layer_nc(layer).set_r_hole_range(itvl.x, itvl.y);
0400     }
0401   }
0402   edm::LogVerbatim("MkFitGeometryESProducer") << ostr.str();
0403 }
0404 
0405 void MkFitGeometryESProducer::findRZBox(const GlobalPoint &gp, float &rmin, float &rmax, float &zmin, float &zmax) {
0406   float r = gp.perp(), z = std::abs(gp.z());
0407   rmax = std::max(r, rmax);
0408   rmin = std::min(r, rmin);
0409   zmax = std::max(z, zmax);
0410   zmin = std::min(z, zmin);
0411 }
0412 
0413 void MkFitGeometryESProducer::aggregateMaterialInfo(mkfit::TrackerInfo &trk_info,
0414                                                     MaterialHistogram &material_histogram) {
0415   //from histogram (vector of tuples) to grid
0416   for (int i = 0; i < trk_info.mat_nbins_z(); i++) {
0417     for (int j = 0; j < trk_info.mat_nbins_r(); j++) {
0418       const MatHistBin &mhb = material_histogram(i, j);
0419       if (mhb.weight > 0) {
0420         trk_info.material_bbxi(i, j) = mhb.xi / mhb.weight;
0421         trk_info.material_radl(i, j) = mhb.rl / mhb.weight;
0422       }
0423     }
0424   }
0425 }
0426 
0427 void MkFitGeometryESProducer::fillLayers(mkfit::TrackerInfo &trk_info) {
0428   mkfit::rectvec<int> rneighbor_map(trk_info.mat_nbins_z(), trk_info.mat_nbins_r());
0429   mkfit::rectvec<int> zneighbor_map(trk_info.mat_nbins_z(), trk_info.mat_nbins_r());
0430 
0431   for (int im = 0; im < trk_info.n_layers(); ++im) {
0432     const mkfit::LayerInfo &li = trk_info.layer(im);
0433     if (!li.is_barrel() && li.zmax() < 0)
0434       continue;  // neg endcap covered by pos
0435     int rin, rout, zmin, zmax;
0436     rin = trk_info.mat_bin_r(li.rin());
0437     rout = trk_info.mat_bin_r(li.rout()) + 1;
0438     if (li.is_barrel()) {
0439       zmin = 0;
0440       zmax = trk_info.mat_bin_z(std::max(std::abs(li.zmax()), std::abs(li.zmin()))) + 1;
0441     } else {
0442       zmin = trk_info.mat_bin_z(li.zmin());
0443       zmax = trk_info.mat_bin_z(li.zmax()) + 1;
0444     }
0445     for (int i = zmin; i < zmax; i++) {
0446       for (int j = rin; j < rout; j++) {
0447         if (trk_info.material_bbxi(i, j) == 0) {
0448           float distancesqmin = 100000;
0449           for (int i2 = zmin; i2 < zmax; i2++) {
0450             for (int j2 = rin; j2 < rout; j2++) {
0451               if (j == j2 && i == i2)
0452                 continue;
0453               auto mydistsq = (i - i2) * (i - i2) + (j - j2) * (j - j2);
0454               if (mydistsq < distancesqmin && trk_info.material_radl(i2, j2) > 0) {
0455                 distancesqmin = mydistsq;
0456                 zneighbor_map(i, j) = i2;
0457                 rneighbor_map(i, j) = j2;
0458               }
0459             }
0460           }  // can work on speedup here
0461         }
0462       }
0463     }
0464     for (int i = zmin; i < zmax; i++) {
0465       for (int j = rin; j < rout; j++) {
0466         if (trk_info.material_bbxi(i, j) == 0) {
0467           int iN = zneighbor_map(i, j);
0468           int jN = rneighbor_map(i, j);
0469           trk_info.material_bbxi(i, j) = trk_info.material_bbxi(iN, jN);
0470           trk_info.material_radl(i, j) = trk_info.material_radl(iN, jN);
0471         }
0472       }
0473     }
0474   }  //module loop
0475 }
0476 
0477 //------------------------------------------------------------------------------
0478 // clang-format off
0479 namespace {
0480   const float phase1QBins[] = {
0481     // PIXB, TIB, TOB
0482     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,
0483     // PIXE+, TID+, TEC+
0484     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,
0485     // PIXE-, TID-, TEC-
0486     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
0487   };
0488   const float phase2QBins[] = {
0489     // TODO: Review these numbers.
0490     // PIXB, TOB
0491     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,
0492     // PIXE+, TEC+
0493     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,
0494     // PIXE-, TEC-
0495     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
0496   };
0497 }
0498 // clang-format on
0499 //------------------------------------------------------------------------------
0500 
0501 std::unique_ptr<MkFitGeometry> MkFitGeometryESProducer::produce(const TrackerRecoGeometryRecord &iRecord) {
0502   auto trackerInfo = std::make_unique<mkfit::TrackerInfo>();
0503 
0504   trackerGeom_ = &iRecord.get(geomToken_);
0505   trackerTopo_ = &iRecord.get(ttopoToken_);
0506 
0507   const float *qBinDefaults = nullptr;
0508 
0509   // std::string path = "Geometry/TrackerCommonData/data/";
0510   if (trackerGeom_->isThere(GeomDetEnumerators::P1PXB) || trackerGeom_->isThere(GeomDetEnumerators::P1PXEC)) {
0511     edm::LogInfo("MkFitGeometryESProducer") << "Extracting PhaseI geometry";
0512     trackerInfo->create_layers(18, 27, 27);
0513     qBinDefaults = phase1QBins;
0514 
0515     trackerInfo->create_material(300, 300.0f, 120, 120.0f);
0516   } else if (trackerGeom_->isThere(GeomDetEnumerators::P2PXB) || trackerGeom_->isThere(GeomDetEnumerators::P2PXEC) ||
0517              trackerGeom_->isThere(GeomDetEnumerators::P2OTB) || trackerGeom_->isThere(GeomDetEnumerators::P2OTEC)) {
0518     edm::LogInfo("MkFitGeometryESProducer") << "Extracting PhaseII geometry";
0519     layerNrConv_.reset(mkfit::TkLayout::phase2);
0520     trackerInfo->create_layers(16, 22, 22);
0521     qBinDefaults = phase2QBins;
0522     trackerInfo->create_material(300, 300.0f, 120, 120.0f);
0523   } else {
0524     throw cms::Exception("UnimplementedFeature") << "unsupported / unknowen geometry version";
0525   }
0526 
0527   // Prepare layer boundaries for bounding-box search
0528   for (int i = 0; i < trackerInfo->n_layers(); ++i) {
0529     auto &li = trackerInfo->layer_nc(i);
0530     li.set_limits(
0531         std::numeric_limits<float>::max(), 0, std::numeric_limits<float>::max(), -std::numeric_limits<float>::max());
0532     li.reserve_modules(256);
0533   }
0534 
0535   MaterialHistogram material_histogram(trackerInfo->mat_nbins_z(), trackerInfo->mat_nbins_r());
0536 
0537   // This is sort of CMS-phase1 specific ... but fireworks code uses it for PhaseII as well.
0538   addPixBGeometry(*trackerInfo, material_histogram);
0539   addPixEGeometry(*trackerInfo, material_histogram);
0540   addTIBGeometry(*trackerInfo, material_histogram);
0541   addTIDGeometry(*trackerInfo, material_histogram);
0542   addTOBGeometry(*trackerInfo, material_histogram);
0543   addTECGeometry(*trackerInfo, material_histogram);
0544 
0545   // r_in/out kept as squares until here, root them
0546   unsigned int n_mod = 0;
0547   for (int i = 0; i < trackerInfo->n_layers(); ++i) {
0548     auto &li = trackerInfo->layer_nc(i);
0549     li.set_r_in_out(std::sqrt(li.rin()), std::sqrt(li.rout()));
0550     li.set_propagate_to(li.is_barrel() ? li.r_mean() : li.z_mean());
0551     li.set_q_bin(qBinDefaults[i]);
0552     unsigned int maxsid = li.shrink_modules();
0553 
0554     n_mod += maxsid;
0555 
0556     // Make sure the short id fits in the 14 bits...
0557     assert(maxsid < 1u << 13);
0558     assert(n_mod > 0);
0559   }
0560 
0561   // Material grid
0562   aggregateMaterialInfo(*trackerInfo, material_histogram);
0563   fillLayers(*trackerInfo);
0564 
0565   // Propagation configuration
0566   {
0567     using namespace mkfit;
0568     PropagationConfig &pconf = trackerInfo->prop_config_nc();
0569     pconf.backward_fit_to_pca = false;
0570     pconf.finding_requires_propagation_to_hit_pos = true;
0571     pconf.finding_inter_layer_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
0572     if (Config::usePropToPlane)
0573       pconf.finding_intra_layer_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
0574     else
0575       pconf.finding_intra_layer_pflags = PropagationFlags(PF_none);
0576     pconf.backward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
0577     pconf.forward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
0578     pconf.seed_fit_pflags = PropagationFlags(PF_none);
0579     pconf.pca_prop_pflags = PropagationFlags(PF_none);
0580     pconf.apply_tracker_info(trackerInfo.get());
0581   }
0582 
0583 #ifdef DUMP_MKF_GEO
0584   printf("Total number of modules %u, 14-bits fit up to %u modules\n", n_mod, 1u << 13);
0585 #endif
0586 
0587   return std::make_unique<MkFitGeometry>(iRecord.get(geomToken_),
0588                                          iRecord.get(trackerToken_),
0589                                          iRecord.get(ttopoToken_),
0590                                          std::move(trackerInfo),
0591                                          layerNrConv_);
0592 }
0593 
0594 DEFINE_FWK_EVENTSETUP_MODULE(MkFitGeometryESProducer);