TkTrackingRegions/plugins/PixelInactiveAreaFinder.cc

0001 #include "PixelInactiveAreaFinder.h"
0002
0003 #include "FWCore/Utilities/interface/VecArray.h"
0004 #include "FWCore/Utilities/interface/transform.h"
0005 #include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
0006 #include "DataFormats/Math/interface/deltaPhi.h"
0007 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0008 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0009
0010 #include "TrackingTools/TransientTrackingRecHit/interface/SeedingLayerSetsLooper.h"
0011
0012 #include <fstream>
0013 #include <queue>
0014 #include <algorithm>
0015
0016 std::ostream& operator<<(std::ostream& os, SeedingLayerSetsBuilder::SeedingLayerId layer) {
0017   if (std::get<0>(layer) == GeomDetEnumerators::PixelBarrel) {
0018     os << "BPix";
0019   } else {
0020     os << "FPix";
0021   }
0022   os << std::get<2>(layer);
0023   if (std::get<1>(layer) == TrackerDetSide::PosEndcap) {
0024     os << "_pos";
0025   } else if (std::get<1>(layer) == TrackerDetSide::NegEndcap) {
0026     os << "_neg";
0027   }
0028   return os;
0029 }
0030
0031 namespace {
0032   using LayerPair = std::pair<SeedingLayerSetsBuilder::SeedingLayerId, SeedingLayerSetsBuilder::SeedingLayerId>;
0033   using ActiveLayerSetToInactiveSetsMap = std::map<LayerPair, edm::VecArray<LayerPair, 5>>;
0034   using Stream = std::stringstream;
0035   using Span_t = std::pair<float, float>;
0036
0037   ActiveLayerSetToInactiveSetsMap createActiveToInactiveMap() {
0038     ActiveLayerSetToInactiveSetsMap map;
0039
0040     auto bpix = [](int layer) {
0041       return SeedingLayerSetsBuilder::SeedingLayerId(GeomDetEnumerators::PixelBarrel, TrackerDetSide::Barrel, layer);
0042     };
0043     auto fpix_pos = [](int disk) {
0044       return SeedingLayerSetsBuilder::SeedingLayerId(GeomDetEnumerators::PixelEndcap, TrackerDetSide::PosEndcap, disk);
0045     };
0046     auto fpix_neg = [](int disk) {
0047       return SeedingLayerSetsBuilder::SeedingLayerId(GeomDetEnumerators::PixelEndcap, TrackerDetSide::NegEndcap, disk);
0048     };
0049
0050     auto add_permutations = [&](std::array<SeedingLayerSetsBuilder::SeedingLayerId, 4> quads) {
0051       do {
0052         // skip permutations like BPix2+BPix1 or FPix1+BPix1
0053         // operator> works automatically
0054         if (quads[0] > quads[1] || quads[2] > quads[3])
0055           continue;
0056
0057         map[std::make_pair(quads[0], quads[1])].emplace_back(quads[2], quads[3]);
0058       } while (std::next_permutation(quads.begin(), quads.end()));
0059     };
0060
0061     // 4 barrel
0062     add_permutations({{bpix(1), bpix(2), bpix(3), bpix(4)}});
0063
0064     // 3 barrel, 1 forward
0065     add_permutations({{bpix(1), bpix(2), bpix(3), fpix_pos(1)}});
0066     add_permutations({{bpix(1), bpix(2), bpix(3), fpix_neg(1)}});
0067
0068     // 2 barrel, 2 forward
0069     add_permutations({{bpix(1), bpix(2), fpix_pos(1), fpix_pos(2)}});
0070     add_permutations({{bpix(1), bpix(2), fpix_neg(1), fpix_neg(2)}});
0071
0072     // 1 barrel, 3 forward
0073     add_permutations({{bpix(1), fpix_pos(1), fpix_pos(2), fpix_pos(3)}});
0074     add_permutations({{bpix(1), fpix_neg(1), fpix_neg(2), fpix_neg(3)}});
0075
0076 #ifdef EDM_ML_DEBUG
0077     LogDebug("PixelInactiveAreaFinder") << "Active to inactive mapping";
0078     for (const auto& elem : map) {
0079       std::stringstream ss;
0080       for (const auto& layerPair : elem.second) {
0081         ss << layerPair.first << "+" << layerPair.second << ",";
0082       }
0083       LogTrace("PixelInactiveAreaFinder") << " " << elem.first.first << "+" << elem.first.second << " => " << ss.str();
0084     }
0085 #endif
0086
0087     return map;
0088   }
0089
0090   void detGroupSpanInfo(const PixelInactiveAreaFinder::DetGroupSpan& cspan, Stream& ss) {
0091     using std::fixed;
0092     using std::left;
0093     using std::noshowpos;
0094     using std::right;
0095     using std::setfill;
0096     using std::setprecision;
0097     using std::setw;
0098     using std::showpos;
0099     std::string deli = "; ";
0100     ss << "subdetid:[" << cspan.subdetId << "]" << deli;
0101     if (cspan.subdetId == PixelSubdetector::PixelBarrel) {
0102       ss << "layer:[" << cspan.layer << "]" << deli;
0103     } else {
0104       ss << "disk:[" << cspan.disk << "]" << deli;
0105     }
0106     ss
0107         //<< setfill(' ') << setw(36) << " "
0108         << setprecision(16) << showpos << "phi:<" << right << setw(12) << cspan.phiSpan.first << "," << left << setw(12)
0109         << cspan.phiSpan.second << ">" << deli << "z:<" << right << setw(7) << cspan.zSpan.first << "," << left
0110         << setw(7) << cspan.zSpan.second << ">" << deli << noshowpos << "r:<" << right << setw(10) << cspan.rSpan.first
0111         << "," << left << setw(10) << cspan.rSpan.second << ">" << deli;
0112   }
0113   void printOverlapSpans(const PixelInactiveAreaFinder::InactiveAreas& areasLayers) {
0114     auto spansLayerSets = areasLayers.spansAndLayerSets(GlobalPoint(0, 0, 0), std::numeric_limits<float>::infinity());
0115
0116     Stream ss;
0117     for (auto const& spansLayers : spansLayerSets) {
0118       ss << "Overlapping detGroups:\n";
0119       for (auto const& cspan : spansLayers.first) {
0120         detGroupSpanInfo(cspan, ss);
0121         ss << std::endl;
0122       }
0123     }
0124     edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0125   }
0126
0127   // Functions for finding bad detGroups
0128   bool phiRangesOverlap(const Span_t& phiSpanA, const Span_t& phiSpanB) {
0129     float x1, x2, y1, y2;
0130     std::tie(x1, x2) = phiSpanA;
0131     std::tie(y1, y2) = phiSpanB;
0132     // assuming phi ranges are [x1,x2] and [y1,y2] and xi,yi in [-pi,pi]
0133     if (x1 <= x2 && y1 <= y2) {
0134       return x1 <= y2 && y1 <= x2;
0135     } else if ((x1 > x2 && y1 <= y2) || (y1 > y2 && x1 <= x2)) {
0136       return y1 <= x2 || x1 <= y2;
0137     } else if (x1 > x2 && y1 > y2) {
0138       return true;
0139     } else {
0140       return false;
0141     }
0142   }
0143
0144   // Functions for finding ranges that detGroups cover
0145   bool phiMoreClockwise(float phiA, float phiB) {
0146     // return true if a is more clockwise than b
0147     return reco::deltaPhi(phiA, phiB) <= 0.f;
0148   }
0149   bool phiMoreCounterclockwise(float phiA, float phiB) {
0150     // return true if a is more counterclockwise than b
0151     return reco::deltaPhi(phiA, phiB) >= 0.f;
0152   }
0153
0154   // Functions for findind overlapping functions
0155   float zAxisIntersection(const float zrPointA[2], const float zrPointB[2]) {
0156     return (zrPointB[0] - zrPointA[0]) / (zrPointB[1] - zrPointA[1]) * (-zrPointA[1]) + zrPointA[0];
0157   }
0158   bool getZAxisOverlapRangeBarrel(const PixelInactiveAreaFinder::DetGroupSpan& cspanA,
0159                                   const PixelInactiveAreaFinder::DetGroupSpan& cspanB,
0160                                   std::pair<float, float>& range) {
0161     PixelInactiveAreaFinder::DetGroupSpan cspanUpper;
0162     PixelInactiveAreaFinder::DetGroupSpan cspanLower;
0163     if (cspanA.rSpan.second < cspanB.rSpan.first) {
0164       cspanLower = cspanA;
0165       cspanUpper = cspanB;
0166     } else if (cspanA.rSpan.first > cspanB.rSpan.second) {
0167       cspanUpper = cspanA;
0168       cspanLower = cspanB;
0169     } else {
0170       return false;
0171     }
0172     float lower = 0;
0173     float upper = 0;
0174     if (cspanUpper.zSpan.second < cspanLower.zSpan.first) {
0175       // lower intersectionpoint, point = {z,r} in cylindrical coordinates
0176       const float pointUpperDetGroupL[2] = {cspanUpper.zSpan.second, cspanUpper.rSpan.second};
0177       const float pointLowerDetGroupL[2] = {cspanLower.zSpan.first, cspanLower.rSpan.first};
0178       lower = zAxisIntersection(pointUpperDetGroupL, pointLowerDetGroupL);
0179       // upper intersectionpoint
0180       const float pointUpperDetGroupU[2] = {cspanUpper.zSpan.first, cspanUpper.rSpan.first};
0181       const float pointLowerDetGroupU[2] = {cspanLower.zSpan.second, cspanLower.rSpan.second};
0182       upper = zAxisIntersection(pointUpperDetGroupU, pointLowerDetGroupU);
0183     } else if (cspanUpper.zSpan.first <= cspanLower.zSpan.second && cspanLower.zSpan.first <= cspanUpper.zSpan.second) {
0184       // lower intersectionpoint, point = {z,r} in cylindrical coordinates
0185       const float pointUpperDetGroupL[2] = {cspanUpper.zSpan.second, cspanUpper.rSpan.first};
0186       const float pointLowerDetGroupL[2] = {cspanLower.zSpan.first, cspanLower.rSpan.second};
0187       lower = zAxisIntersection(pointUpperDetGroupL, pointLowerDetGroupL);
0188       // upper intersectionpoint
0189       const float pointUpperDetGroupU[2] = {cspanUpper.zSpan.first, cspanUpper.rSpan.first};
0190       const float pointLowerDetGroupU[2] = {cspanLower.zSpan.second, cspanLower.rSpan.second};
0191       upper = zAxisIntersection(pointUpperDetGroupU, pointLowerDetGroupU);
0192     } else if (cspanUpper.zSpan.first > cspanLower.zSpan.second) {
0193       // lower intersectionpoint, point = {z,r} in cylindrical coordinates
0194       const float pointUpperDetGroupL[2] = {cspanUpper.zSpan.second, cspanUpper.rSpan.first};
0195       const float pointLowerDetGroupL[2] = {cspanLower.zSpan.first, cspanLower.rSpan.second};
0196       lower = zAxisIntersection(pointUpperDetGroupL, pointLowerDetGroupL);
0197       // upper intersectionpoint
0198       const float pointUpperDetGroupU[2] = {cspanUpper.zSpan.first, cspanUpper.rSpan.second};
0199       const float pointLowerDetGroupU[2] = {cspanLower.zSpan.second, cspanLower.rSpan.first};
0200       upper = zAxisIntersection(pointUpperDetGroupU, pointLowerDetGroupU);
0201     } else {
0202       //something wrong
0203       return false;
0204     }
0205     range = std::pair<float, float>(lower, upper);
0206     return true;
0207   }
0208   bool getZAxisOverlapRangeEndcap(const PixelInactiveAreaFinder::DetGroupSpan& cspanA,
0209                                   const PixelInactiveAreaFinder::DetGroupSpan& cspanB,
0210                                   std::pair<float, float>& range) {
0211     // While on left hand side of pixel detector
0212     PixelInactiveAreaFinder::DetGroupSpan cspanNearer;
0213     PixelInactiveAreaFinder::DetGroupSpan cspanFurther;
0214     float lower = 0;
0215     float upper = 0;
0216     if (cspanA.zSpan.first < 0 && cspanB.zSpan.first < 0) {
0217       if (cspanA.zSpan.second < cspanB.zSpan.first) {
0218         cspanFurther = cspanA;
0219         cspanNearer = cspanB;
0220       } else if (cspanB.zSpan.second < cspanA.zSpan.first) {
0221         cspanFurther = cspanB;
0222         cspanNearer = cspanA;
0223       } else {
0224 #ifdef EDM_ML_DEBUG
0225         LogTrace("PixelInactiveAreaFinder") << "No overlap, same disk propably. Spans:";
0226         Stream ss;
0227         detGroupSpanInfo(cspanA, ss);
0228         ss << std::endl;
0229         detGroupSpanInfo(cspanB, ss);
0230         ss << std::endl;
0231         LogTrace("PixelInactiveAreaFinder") << ss.str();
0232         ss.str(std::string());
0233         LogTrace("PixelInactiveAreaFinder") << "**";
0234 #endif
0235         return false;
0236       }
0237       if (cspanFurther.rSpan.second > cspanNearer.rSpan.first) {
0238         const float pointA[2] = {cspanFurther.zSpan.second, cspanFurther.rSpan.second};
0239         const float pointB[2] = {cspanNearer.zSpan.first, cspanNearer.rSpan.first};
0240         lower = zAxisIntersection(pointA, pointB);
0241         if (cspanFurther.rSpan.first > cspanNearer.rSpan.second) {
0242           const float pointC[2] = {cspanFurther.zSpan.first, cspanFurther.rSpan.first};
0243           const float pointD[2] = {cspanNearer.zSpan.second, cspanFurther.rSpan.second};
0244           upper = zAxisIntersection(pointC, pointD);
0245         } else {
0246           upper = std::numeric_limits<float>::infinity();
0247         }
0248       } else {
0249 #ifdef EDM_ML_DEBUG
0250         LogTrace("PixelInactiveAreaFinder") << "No overlap, further detGroup is lower. Spans:";
0251         Stream ss;
0252         detGroupSpanInfo(cspanA, ss);
0253         ss << std::endl;
0254         detGroupSpanInfo(cspanB, ss);
0255         ss << std::endl;
0256         LogTrace("PixelInactiveAreaFinder") << ss.str();
0257         ss.str(std::string());
0258         LogTrace("PixelInactiveAreaFinder") << "**";
0259 #endif
0260         return false;
0261       }
0262     } else if (cspanA.zSpan.first > 0 && cspanB.zSpan.first > 0) {
0263       if (cspanA.zSpan.first > cspanB.zSpan.second) {
0264         cspanFurther = cspanA;
0265         cspanNearer = cspanB;
0266       } else if (cspanB.zSpan.first > cspanA.zSpan.second) {
0267         cspanFurther = cspanB;
0268         cspanNearer = cspanA;
0269       } else {
0270 #ifdef EDM_ML_DEBUG
0271         LogTrace("PixelInactiveAreaFinder") << "No overlap, same disk propably. Spans:";
0272         Stream ss;
0273         detGroupSpanInfo(cspanA, ss);
0274         ss << std::endl;
0275         detGroupSpanInfo(cspanB, ss);
0276         ss << std::endl;
0277         LogTrace("PixelInactiveAreaFinder") << ss.str();
0278         ss.str(std::string());
0279         LogTrace("PixelInactiveAreaFinder") << "**";
0280 #endif
0281         return false;
0282       }
0283       if (cspanFurther.rSpan.second > cspanNearer.rSpan.first) {
0284         const float pointA[2] = {cspanFurther.zSpan.first, cspanFurther.rSpan.second};
0285         const float pointB[2] = {cspanNearer.zSpan.second, cspanNearer.rSpan.first};
0286         upper = zAxisIntersection(pointA, pointB);
0287         if (cspanFurther.rSpan.first > cspanNearer.rSpan.second) {
0288           const float pointC[2] = {cspanFurther.zSpan.second, cspanFurther.rSpan.first};
0289           const float pointD[2] = {cspanNearer.zSpan.first, cspanFurther.rSpan.second};
0290           lower = zAxisIntersection(pointC, pointD);
0291         } else {
0292           lower = -std::numeric_limits<float>::infinity();
0293         }
0294       } else {
0295 #ifdef EDM_ML_DEBUG
0296         LogTrace("PixelInactiveAreaFinder") << "No overlap, further detGroup lower. Spans:";
0297         Stream ss;
0298         detGroupSpanInfo(cspanA, ss);
0299         ss << std::endl;
0300         detGroupSpanInfo(cspanB, ss);
0301         ss << std::endl;
0302         LogTrace("PixelInactiveAreaFinder") << ss.str();
0303         ss.str(std::string());
0304         LogTrace("PixelInactiveAreaFinder") << "**";
0305 #endif
0306         return false;
0307       }
0308     } else {
0309 #ifdef EDM_ML_DEBUG
0310       LogTrace("PixelInactiveAreaFinder") << "No overlap, different sides of z axis. Spans:";
0311       Stream ss;
0312       detGroupSpanInfo(cspanA, ss);
0313       ss << std::endl;
0314       detGroupSpanInfo(cspanB, ss);
0315       ss << std::endl;
0316       LogTrace("PixelInactiveAreaFinder") << ss.str();
0317       ss.str(std::string());
0318       LogTrace("PixelInactiveAreaFinder") << "**";
0319 #endif
0320       return false;
0321     }
0322     range = std::pair<float, float>(lower, upper);
0323     return true;
0324   }
0325
0326   bool getZAxisOverlapRangeBarrelEndcap(const PixelInactiveAreaFinder::DetGroupSpan& cspanBar,
0327                                         const PixelInactiveAreaFinder::DetGroupSpan& cspanEnd,
0328                                         std::pair<float, float>& range) {
0329     float lower = 0;
0330     float upper = 0;
0331     if (cspanEnd.rSpan.second > cspanBar.rSpan.first) {
0332       if (cspanEnd.zSpan.second < cspanBar.zSpan.first) {
0333         // if we are on the left hand side of pixel detector
0334         const float pointA[2] = {cspanEnd.zSpan.second, cspanEnd.rSpan.second};
0335         const float pointB[2] = {cspanBar.zSpan.first, cspanBar.rSpan.first};
0336         lower = zAxisIntersection(pointA, pointB);
0337         if (cspanEnd.rSpan.first > cspanBar.rSpan.second) {
0338           // if does not overlap, then there is also upper limit
0339           const float pointC[2] = {cspanEnd.zSpan.first, cspanEnd.rSpan.first};
0340           const float pointD[2] = {cspanBar.zSpan.second, cspanBar.rSpan.second};
0341           upper = zAxisIntersection(pointC, pointD);
0342         } else {
0343           upper = std::numeric_limits<float>::infinity();
0344         }
0345       } else if (cspanEnd.zSpan.first > cspanBar.zSpan.second) {
0346         // if we are on the right hand side of pixel detector
0347         const float pointA[2] = {cspanEnd.zSpan.first, cspanEnd.rSpan.second};
0348         const float pointB[2] = {cspanBar.zSpan.second, cspanBar.rSpan.first};
0349         upper = zAxisIntersection(pointA, pointB);
0350         if (cspanEnd.rSpan.first > cspanBar.rSpan.second) {
0351           const float pointC[2] = {cspanEnd.zSpan.second, cspanEnd.rSpan.first};
0352           const float pointD[2] = {cspanBar.zSpan.first, cspanBar.rSpan.second};
0353           lower = zAxisIntersection(pointC, pointD);
0354         } else {
0355           lower = -std::numeric_limits<float>::infinity();
0356         }
0357       } else {
0358         return false;
0359       }
0360     } else {
0361       return false;
0362     }
0363     range = std::pair<float, float>(lower, upper);
0364     return true;
0365   }
0366 }  // namespace
0367
0368 std::vector<
0369     std::pair<edm::VecArray<PixelInactiveAreaFinder::Area, 2>, std::vector<PixelInactiveAreaFinder::LayerSetIndex>>>
0370 PixelInactiveAreaFinder::InactiveAreas::areasAndLayerSets(const GlobalPoint& point, float zwidth) const {
0371   auto spansLayerSets = spansAndLayerSets(point, zwidth);
0372
0373   // TODO: try to remove this conversion...
0374   std::vector<std::pair<VecArray2<Area>, std::vector<LayerSetIndex>>> ret;
0375   for (auto& item : spansLayerSets) {
0376     auto& innerSpan = item.first[0];
0377     auto& outerSpan = item.first[1];
0378     VecArray2<Area> areas;
0379     areas.emplace_back(innerSpan.rSpan.first,
0380                        innerSpan.rSpan.second,
0381                        innerSpan.phiSpan.first,
0382                        innerSpan.phiSpan.second,
0383                        innerSpan.zSpan.first,
0384                        innerSpan.zSpan.second);
0385     areas.emplace_back(outerSpan.rSpan.first,
0386                        outerSpan.rSpan.second,
0387                        outerSpan.phiSpan.first,
0388                        outerSpan.phiSpan.second,
0389                        outerSpan.zSpan.first,
0390                        outerSpan.zSpan.second);
0391     ret.emplace_back(areas, std::move(item.second));
0392   }
0393
0394   return ret;
0395 }
0396
0397 std::vector<std::pair<edm::VecArray<PixelInactiveAreaFinder::DetGroupSpan, 2>,
0398                       std::vector<PixelInactiveAreaFinder::LayerSetIndex>>>
0399 PixelInactiveAreaFinder::InactiveAreas::spansAndLayerSets(const GlobalPoint& point, float zwidth) const {
0400   // TODO: in the future use 2D-r for the origin for the phi overlap check
0401   const float zmin = point.z() - zwidth;
0402   const float zmax = point.z() + zwidth;
0403
0404   std::vector<std::pair<VecArray2<DetGroupSpan>, std::vector<LayerSetIndex>>> ret;
0405
0406   LogDebug("PixelInactiveAreaFinder") << "Origin at " << point.x() << "," << point.y() << "," << point.z()
0407                                       << " z half width " << zwidth;
0408
0409   for (LayerSetIndex i = 0, end = inactiveLayerPairIndices_->size(); i < end; ++i) {
0410     const auto& layerIdxPair = (*inactiveLayerPairIndices_)[i];
0411     const auto& innerSpans = inactiveSpans_[layerIdxPair.first];
0412     const auto& outerSpans = inactiveSpans_[layerIdxPair.second];
0413
0414     for (const auto& innerSpan : innerSpans) {
0415       for (const auto& outerSpan : outerSpans) {
0416         if (phiRangesOverlap(innerSpan.phiSpan, outerSpan.phiSpan)) {
0417           std::pair<float, float> range(0, 0);
0418
0419           bool zOverlap = false;
0420           const auto innerDet = std::get<0>((*inactiveLayers_)[layerIdxPair.first]);
0421           const auto outerDet = std::get<0>((*inactiveLayers_)[layerIdxPair.second]);
0422           if (innerDet == GeomDetEnumerators::PixelBarrel) {
0423             if (outerDet == GeomDetEnumerators::PixelBarrel)
0424               zOverlap = getZAxisOverlapRangeBarrel(innerSpan, outerSpan, range);
0425             else
0426               zOverlap = getZAxisOverlapRangeBarrelEndcap(innerSpan, outerSpan, range);
0427           } else {
0428             if (outerDet == GeomDetEnumerators::PixelEndcap)
0429               zOverlap = getZAxisOverlapRangeEndcap(innerSpan, outerSpan, range);
0430             else
0431               throw cms::Exception("LogicError") << "Forward->barrel transition is not supported";
0432           }
0433
0434           if (zOverlap && zmin <= range.second && range.first <= zmax) {
0435 #ifdef EDM_ML_DEBUG
0436             Stream ss;
0437             for (auto ind : (*layerSetIndexInactiveToActive_)[i]) {
0438               ss << ind << ",";
0439             }
0440             ss << "\n  ";
0441             detGroupSpanInfo(innerSpan, ss);
0442             ss << "\n  ";
0443             detGroupSpanInfo(outerSpan, ss);
0444             LogTrace("PixelInactiveAreaFinder") << " adding areas for active layer sets " << ss.str();
0445 #endif
0446             VecArray2<DetGroupSpan> vec;
0447             vec.emplace_back(innerSpan);
0448             vec.emplace_back(outerSpan);
0449             ret.emplace_back(std::move(vec), (*layerSetIndexInactiveToActive_)[i]);
0450           }
0451         }
0452       }
0453     }
0454   }
0455
0456   return ret;
0457 }
0458
0459 PixelInactiveAreaFinder::PixelInactiveAreaFinder(
0460     const edm::ParameterSet& iConfig,
0461     const std::vector<SeedingLayerSetsBuilder::SeedingLayerId>& seedingLayers,
0462     const SeedingLayerSetsLooper& seedingLayerSetsLooper,
0463     edm::ConsumesCollector&& iC)
0464     : debug_(iConfig.getUntrackedParameter<bool>("debug")),
0465       createPlottingFiles_(iConfig.getUntrackedParameter<bool>("createPlottingFiles")),
0466       ignoreSingleFPixPanelModules_(iConfig.getParameter<bool>("ignoreSingleFPixPanelModules")),
0467       inactivePixelDetectorTokens_(
0468           edm::vector_transform(iConfig.getParameter<std::vector<edm::InputTag>>("inactivePixelDetectorLabels"),
0469                                 [&](const auto& tag) { return iC.consumes<DetIdCollection>(tag); })),
0470       badPixelFEDChannelsTokens_(
0471           edm::vector_transform(iConfig.getParameter<std::vector<edm::InputTag>>("badPixelFEDChannelCollectionLabels"),
0472                                 [&](const auto& tag) { return iC.consumes<PixelFEDChannelCollection>(tag); })),
0473       trackerGeometryToken_(iC.esConsumes()),
0474       trackerTopologyToken_(iC.esConsumes()),
0475       pixelQualityToken_(iC.esConsumes()) {
0476 #ifdef EDM_ML_DEBUG
0477   for (const auto& layer : seedingLayers) {
0478     LogTrace("PixelInactiveAreaFinder") << "Input layer subdet " << std::get<0>(layer) << " side "
0479                                         << static_cast<unsigned int>(std::get<1>(layer)) << " layer "
0480                                         << std::get<2>(layer);
0481   }
0482 #endif
0483
0484   auto findOrAdd = [&](SeedingLayerId layer) -> unsigned short {
0485     auto found = std::find(inactiveLayers_.cbegin(), inactiveLayers_.cend(), layer);
0486     if (found == inactiveLayers_.cend()) {
0487       auto ret = inactiveLayers_.size();
0488       inactiveLayers_.push_back(layer);
0489       return ret;
0490     }
0491     return std::distance(inactiveLayers_.cbegin(), found);
0492   };
0493
0494   // mapping from active layer pairs to inactive layer pairs
0495   const auto activeToInactiveMap = createActiveToInactiveMap();
0496
0497   // convert input layer pairs (that are for active layers) to layer
0498   // pairs to look for inactive areas
0499   LayerSetIndex i = 0;
0500   for (const auto& layerSet : seedingLayerSetsLooper.makeRange(seedingLayers)) {
0501     assert(layerSet.size() == 2);
0502     auto found = activeToInactiveMap.find(std::make_pair(layerSet[0], layerSet[1]));
0503     if (found == activeToInactiveMap.end()) {
0504       throw cms::Exception("Configuration")
0505           << "Encountered layer pair " << layerSet[0] << "+" << layerSet[1]
0506           << " not found from the internal 'active layer pairs' to 'inactive layer pairs' mapping; either fix the "
0507              "input or the mapping (in PixelInactiveAreaFinder.cc)";
0508     }
0509
0510     LogTrace("PixelInactiveAreaFinder") << "Input layer set " << layerSet[0] << "+" << layerSet[1];
0511     for (const auto& inactiveLayerSet : found->second) {
0512       auto innerInd = findOrAdd(inactiveLayerSet.first);
0513       auto outerInd = findOrAdd(inactiveLayerSet.second);
0514
0515       auto found = std::find(
0516           inactiveLayerSetIndices_.cbegin(), inactiveLayerSetIndices_.cend(), std::make_pair(innerInd, outerInd));
0517       if (found == inactiveLayerSetIndices_.end()) {
0518         inactiveLayerSetIndices_.emplace_back(innerInd, outerInd);
0519         layerSetIndexInactiveToActive_.push_back(std::vector<LayerSetIndex>{i});
0520       } else {
0521         layerSetIndexInactiveToActive_.at(std::distance(inactiveLayerSetIndices_.cbegin(), found))
0522             .push_back(i);  // TODO: move to operator[] once finished
0523       }
0524
0525       LogTrace("PixelInactiveAreaFinder")
0526           << " inactive layer set " << inactiveLayerSet.first << "+" << inactiveLayerSet.second;
0527     }
0528
0529     ++i;
0530   }
0531
0532 #ifdef EDM_ML_DEBUG
0533   LogDebug("PixelInactiveAreaFinder") << "All inactive layer sets";
0534   for (const auto& idxPair : inactiveLayerSetIndices_) {
0535     LogTrace("PixelInactiveAreaFinder") << " " << inactiveLayers_[idxPair.first] << "+"
0536                                         << inactiveLayers_[idxPair.second];
0537   }
0538 #endif
0539 }
0540
0541 void PixelInactiveAreaFinder::fillDescriptions(edm::ParameterSetDescription& desc) {
0542   desc.add<std::vector<edm::InputTag>>("inactivePixelDetectorLabels",
0543                                        std::vector<edm::InputTag>{{edm::InputTag("siPixelDigis")}})
0544       ->setComment("One or more DetIdCollections of modules to mask on the fly for a given event");
0545   desc.add<std::vector<edm::InputTag>>("badPixelFEDChannelCollectionLabels",
0546                                        std::vector<edm::InputTag>{{edm::InputTag("siPixelDigis")}})
0547       ->setComment("One or more PixelFEDChannelCollections of modules+ROCs to mask on the fly for a given event");
0548   desc.add<bool>("ignoreSingleFPixPanelModules", false);
0549
0550   desc.addUntracked<bool>("debug", false);
0551   desc.addUntracked<bool>("createPlottingFiles", false);
0552 }
0553
0554 PixelInactiveAreaFinder::InactiveAreas PixelInactiveAreaFinder::inactiveAreas(const edm::Event& iEvent,
0555                                                                               const edm::EventSetup& iSetup) {
0556   // Set data to handles
0557   trackerGeometry_ = &iSetup.getData(trackerGeometryToken_);
0558   trackerTopology_ = &iSetup.getData(trackerTopologyToken_);
0559
0560   // assign data to instance variables
0561   updatePixelDets(iSetup);
0562
0563   // clear the list of bad pixel modules at each event!
0564   badPixelDetsBarrel_.clear();
0565   badPixelDetsEndcap_.clear();
0566   getBadPixelDets(iEvent, iSetup);
0567
0568   //write files for plotting
0569   if (createPlottingFiles_) {
0570     createPlottingFiles();
0571   }
0572
0573   // find detGroupSpans ie phi,r,z limits for detector detGroups that are not working
0574   // returns pair where first is barrel spans and second endcap spans
0575   DetGroupSpanContainerPair cspans = detGroupSpans();
0576
0577   // map spans to a vector with consistent indexing with inactiveLayers_ and inactiveLayerSetIndices_
0578   // TODO: try to move the inner logic towards this direction as well
0579   std::vector<DetGroupSpanContainer> spans(inactiveLayers_.size());
0580
0581   auto doWork = [&](const DetGroupSpanContainer& container) {
0582     for (const auto& span : container) {
0583       const auto subdet = span.subdetId == PixelSubdetector::PixelBarrel ? GeomDetEnumerators::PixelBarrel
0584                                                                          : GeomDetEnumerators::PixelEndcap;
0585       const auto side = (subdet == GeomDetEnumerators::PixelBarrel
0586                              ? TrackerDetSide::Barrel
0587                              : (span.zSpan.first < 0 ? TrackerDetSide::NegEndcap : TrackerDetSide::PosEndcap));
0588       const auto layer = subdet == GeomDetEnumerators::PixelBarrel ? span.layer : span.disk;
0589       auto found = std::find(inactiveLayers_.begin(), inactiveLayers_.end(), SeedingLayerId(subdet, side, layer));
0590       if (found != inactiveLayers_.end()) {  // it is possible that this layer is ignored by the configuration
0591         spans[std::distance(inactiveLayers_.begin(), found)].push_back(span);
0592       }
0593     }
0594   };
0595   doWork(cspans.first);
0596   doWork(cspans.second);
0597
0598   auto ret =
0599       InactiveAreas(&inactiveLayers_, std::move(spans), &inactiveLayerSetIndices_, &layerSetIndexInactiveToActive_);
0600
0601   if (debug_) {
0602     printOverlapSpans(ret);
0603   }
0604
0605   return ret;
0606 }
0607
0608 // Functions for fetching date from handles
0609 void PixelInactiveAreaFinder::updatePixelDets(const edm::EventSetup& iSetup) {
0610   if (!geometryWatcher_.check(iSetup))
0611     return;
0612
0613   // deduce the number of ladders per layer and the number of modules per ladder
0614   {
0615     // sanity checks
0616     if (trackerGeometry_->numberOfLayers(PixelSubdetector::PixelBarrel) != 4) {
0617       throw cms::Exception("NotImplemented")
0618           << "This module supports only a detector with 4 pixel barrel layers, the current geometry has "
0619           << trackerGeometry_->numberOfLayers(PixelSubdetector::PixelBarrel);
0620     }
0621     if (trackerGeometry_->numberOfLayers(PixelSubdetector::PixelEndcap) != 3) {
0622       throw cms::Exception("NotImplemented")
0623           << "This module supports only a detector with 3 pixel forward disks, the current geometry has "
0624           << trackerGeometry_->numberOfLayers(PixelSubdetector::PixelEndcap);
0625     }
0626
0627     std::array<std::array<unsigned short, 100>, 4> counts = {};  // assume at most 100 ladders per layer
0628     for (const auto& det : trackerGeometry_->detsPXB()) {
0629       const auto layer = trackerTopology_->layer(det->geographicalId());
0630       const auto ladder = trackerTopology_->pxbLadder(det->geographicalId());
0631       if (ladder >= 100) {
0632         throw cms::Exception("LogicError")
0633             << "Got a ladder with number " << ladder
0634             << " while the expected maximum was 100; either something is wrong or the maximum has to be increased.";
0635       }
0636       counts[layer - 1][ladder - 1] += 1;  // numbering of layer and ladder starts at 1
0637     }
0638
0639     // take number of modules per ladder from the first ladder of the first layer (such better exist)
0640     // other ladders better have the same number
0641     nModulesPerLadder = counts[0][0];
0642     if (nModulesPerLadder == 0) {
0643       throw cms::Exception("LogicError") << "Ladder 1 of layer 1 has 0 modules, something fishy is going on.";
0644     }
0645
0646     LogDebug("PixelInactiveAreaFinder") << "Number of modules per ladder " << nModulesPerLadder
0647                                         << "; below are number of ladders per layer";
0648
0649     // number of ladders
0650     for (unsigned layer = 0; layer < 4; ++layer) {
0651       nBPixLadders[layer] =
0652           std::count_if(counts[layer].begin(), counts[layer].end(), [](unsigned short val) { return val > 0; });
0653       LogTrace("PixelInactiveAreaFinder")
0654           << "BPix layer " << (layer + 1) << " has " << nBPixLadders[layer] << " ladders";
0655
0656       auto fail = std::count_if(counts[layer].begin(), counts[layer].end(), [&](unsigned short val) {
0657         return val != nModulesPerLadder && val > 0;
0658       });
0659       if (fail != 0) {
0660         throw cms::Exception("LogicError")
0661             << "Layer " << (layer + 1) << " had " << fail
0662             << " ladders whose number of modules/ladder differed from the ladder 1 of layer 1 (" << nModulesPerLadder
0663             << "). Something fishy is going on.";
0664       }
0665     }
0666   }
0667
0668   // don't bother with the rest if not needed
0669   if (!createPlottingFiles_)
0670     return;
0671
0672   pixelDetsBarrel_.clear();
0673   pixelDetsEndcap_.clear();
0674
0675   for (auto const& geomDetPtr : trackerGeometry_->detsPXB()) {
0676     if (geomDetPtr->geographicalId().subdetId() == PixelSubdetector::PixelBarrel) {
0677       pixelDetsBarrel_.push_back(geomDetPtr->geographicalId().rawId());
0678     }
0679   }
0680   for (auto const& geomDetPtr : trackerGeometry_->detsPXF()) {
0681     if (geomDetPtr->geographicalId().subdetId() == PixelSubdetector::PixelEndcap) {
0682       pixelDetsEndcap_.push_back(geomDetPtr->geographicalId().rawId());
0683     }
0684   }
0685   std::sort(pixelDetsBarrel_.begin(), pixelDetsBarrel_.end());
0686   std::sort(pixelDetsEndcap_.begin(), pixelDetsEndcap_.end());
0687 }
0688 void PixelInactiveAreaFinder::getBadPixelDets(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0689   auto addDetId = [&](const auto id) {
0690     const auto detid = DetId(id);
0691     const auto subdet = detid.subdetId();
0692     if (subdet == PixelSubdetector::PixelBarrel) {
0693       badPixelDetsBarrel_.push_back(detid.rawId());
0694     } else if (subdet == PixelSubdetector::PixelEndcap) {
0695       badPixelDetsEndcap_.push_back(detid.rawId());
0696     }
0697   };
0698
0699   // SiPixelQuality
0700   auto const& pixelQuality = iSetup.getData(pixelQualityToken_);
0701
0702   for (auto const& disabledModule : pixelQuality.getBadComponentList()) {
0703     addDetId(disabledModule.DetID);
0704   }
0705
0706   // dynamic bad modules
0707   for (const auto& token : inactivePixelDetectorTokens_) {
0708     edm::Handle<DetIdCollection> detIds;
0709     iEvent.getByToken(token, detIds);
0710     for (const auto& id : *detIds) {
0711       addDetId(id);
0712     }
0713   }
0714
0715   // dynamic bad ROCs ("Fed25")
0716   // TODO: consider moving to finer-grained areas for inactive ROCs
0717   for (const auto& token : badPixelFEDChannelsTokens_) {
0718     edm::Handle<PixelFEDChannelCollection> pixelFEDChannelCollectionHandle;
0719     iEvent.getByToken(token, pixelFEDChannelCollectionHandle);
0720     for (const auto& disabledChannels : *pixelFEDChannelCollectionHandle) {
0721       addDetId(disabledChannels.detId());
0722     }
0723   }
0724
0725   // remove duplicates
0726   std::sort(badPixelDetsBarrel_.begin(), badPixelDetsBarrel_.end());
0727   std::sort(badPixelDetsEndcap_.begin(), badPixelDetsEndcap_.end());
0728   badPixelDetsBarrel_.erase(std::unique(badPixelDetsBarrel_.begin(), badPixelDetsBarrel_.end()),
0729                             badPixelDetsBarrel_.end());
0730   badPixelDetsEndcap_.erase(std::unique(badPixelDetsEndcap_.begin(), badPixelDetsEndcap_.end()),
0731                             badPixelDetsEndcap_.end());
0732 }
0733 // Printing functions
0734 void PixelInactiveAreaFinder::detInfo(const det_t& det, Stream& ss) {
0735   using std::fixed;
0736   using std::left;
0737   using std::noshowpos;
0738   using std::right;
0739   using std::setfill;
0740   using std::setprecision;
0741   using std::setw;
0742   using std::showpos;
0743   using std::tie;
0744   std::string deli = "; ";
0745   ss << "id:[" << det << "]" << deli;
0746   ss << "subdetid:[" << DetId(det).subdetId() << "]" << deli;
0747   if (DetId(det).subdetId() == PixelSubdetector::PixelBarrel) {
0748     unsigned int layer = trackerTopology_->pxbLayer(DetId(det));
0749     unsigned int ladder = trackerTopology_->pxbLadder(DetId(det));
0750     unsigned int module = trackerTopology_->pxbModule(DetId(det));
0751     ss << "layer:[" << layer << "]" << deli << "ladder:[" << right << setw(2) << ladder << "]" << deli << "module:["
0752        << module << "]" << deli;
0753   } else if (DetId(det).subdetId() == PixelSubdetector::PixelEndcap) {
0754     unsigned int disk = trackerTopology_->pxfDisk(DetId(det));
0755     unsigned int blade = trackerTopology_->pxfBlade(DetId(det));
0756     unsigned int panel = trackerTopology_->pxfPanel(DetId(det));
0757     ss << left << setw(6) << "disk:"
0758        << "[" << right << disk << "]" << deli << left << setw(7) << "blade:"
0759        << "[" << setw(2) << right << blade << "]" << deli << left << setw(7) << "panel:"
0760        << "[" << right << panel << "]" << deli;
0761   }
0762   float phiA, phiB, zA, zB, rA, rB;
0763   auto detSurface = trackerGeometry_->idToDet(DetId(det))->surface();
0764   tie(phiA, phiB) = detSurface.phiSpan();
0765   tie(zA, zB) = detSurface.zSpan();
0766   tie(rA, rB) = detSurface.rSpan();
0767   ss << setprecision(16) << fixed << showpos << setfill(' ') << "phi:[" << right << setw(12) << phiA << "," << left
0768      << setw(12) << phiB << "]" << deli << "z:[" << right << setw(7) << zA << "," << left << setw(7) << zB << "]"
0769      << deli << noshowpos << "r:[" << right << setw(10) << rA << "," << left << setw(10) << rB << "]" << deli;
0770 }
0771 void PixelInactiveAreaFinder::printPixelDets() {
0772   edm::LogPrint("PixelInactiveAreaFinder") << "Barrel detectors:";
0773   Stream ss;
0774   for (auto const& det : pixelDetsBarrel_) {
0775     detInfo(det, ss);
0776     edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0777     ss.str(std::string());
0778   }
0779   edm::LogPrint("PixelInactiveAreaFinder") << "Endcap detectors;";
0780   for (auto const& det : pixelDetsEndcap_) {
0781     detInfo(det, ss);
0782     edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0783     ss.str(std::string());
0784   }
0785 }
0786 void PixelInactiveAreaFinder::printBadPixelDets() {
0787   edm::LogPrint("PixelInactiveAreaFinder") << "Bad barrel detectors:";
0788   Stream ss;
0789   for (auto const& det : badPixelDetsBarrel_) {
0790     detInfo(det, ss);
0791     edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0792     ss.str(std::string());
0793   }
0794   edm::LogPrint("PixelInactiveAreaFinder") << "Endcap detectors;";
0795   for (auto const& det : badPixelDetsEndcap_) {
0796     detInfo(det, ss);
0797     edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0798     ss.str(std::string());
0799   }
0800 }
0801 void PixelInactiveAreaFinder::printBadDetGroups() {
0802   DetGroupContainer badDetGroupsBar = badDetGroupsBarrel();
0803   DetGroupContainer badDetGroupsEnd = badDetGroupsEndcap();
0804   Stream ss;
0805   for (auto const& detGroup : badDetGroupsBar) {
0806     ss << std::setfill(' ') << std::left << std::setw(16) << "DetGroup:";
0807     DetGroupSpan cspan;
0808     getPhiSpanBarrel(detGroup, cspan);
0809     getZSpan(detGroup, cspan);
0810     getRSpan(detGroup, cspan);
0811     detGroupSpanInfo(cspan, ss);
0812     ss << std::endl;
0813     for (auto const& det : detGroup) {
0814       detInfo(det, ss);
0815       ss << std::endl;
0816     }
0817     ss << std::endl;
0818   }
0819   for (auto const& detGroup : badDetGroupsEnd) {
0820     ss << std::setfill(' ') << std::left << std::setw(16) << "DetGroup:";
0821     DetGroupSpan cspan;
0822     getPhiSpanEndcap(detGroup, cspan);
0823     getZSpan(detGroup, cspan);
0824     getRSpan(detGroup, cspan);
0825     detGroupSpanInfo(cspan, ss);
0826     ss << std::endl;
0827     for (auto const& det : detGroup) {
0828       detInfo(det, ss);
0829       ss << std::endl;
0830     }
0831     ss << std::endl;
0832   }
0833   edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0834 }
0835 void PixelInactiveAreaFinder::printBadDetGroupSpans() {
0836   DetGroupSpanContainerPair cspans = detGroupSpans();
0837   Stream ss;
0838   for (auto const& cspan : cspans.first) {
0839     detGroupSpanInfo(cspan, ss);
0840     ss << std::endl;
0841   }
0842   for (auto const& cspan : cspans.second) {
0843     detGroupSpanInfo(cspan, ss);
0844     ss << std::endl;
0845   }
0846   edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0847 }
0848 void PixelInactiveAreaFinder::createPlottingFiles() {
0849   // All detectors to file DETECTORS
0850   Stream ss;
0851   std::ofstream fsDet("DETECTORS.txt");
0852   for (auto const& det : pixelDetsBarrel_) {
0853     detInfo(det, ss);
0854     ss << std::endl;
0855   }
0856   edm::LogPrint("PixelInactiveAreaFinder") << "Endcap detectors;";
0857   for (auto const& det : pixelDetsEndcap_) {
0858     detInfo(det, ss);
0859     ss << std::endl;
0860   }
0861   fsDet << ss.rdbuf();
0862   ss.str(std::string());
0863   // Bad detectors
0864   std::ofstream fsBadDet("BADDETECTORS.txt");
0865   for (auto const& det : badPixelDetsBarrel_) {
0866     detInfo(det, ss);
0867     ss << std::endl;
0868   }
0869   for (auto const& det : badPixelDetsEndcap_) {
0870     detInfo(det, ss);
0871     ss << std::endl;
0872   }
0873   fsBadDet << ss.rdbuf();
0874   ss.str(std::string());
0875   // detgroupspans
0876   std::ofstream fsSpans("DETGROUPSPANS.txt");
0877   DetGroupSpanContainerPair cspans = detGroupSpans();
0878   for (auto const& cspan : cspans.first) {
0879     detGroupSpanInfo(cspan, ss);
0880     ss << std::endl;
0881   }
0882   for (auto const& cspan : cspans.second) {
0883     detGroupSpanInfo(cspan, ss);
0884     ss << std::endl;
0885   }
0886   fsSpans << ss.rdbuf();
0887   ss.str(std::string());
0888 }
0889 // Functions for finding bad detGroups
0890 bool PixelInactiveAreaFinder::detWorks(det_t det) {
0891   return std::find(badPixelDetsBarrel_.begin(), badPixelDetsBarrel_.end(), det) == badPixelDetsBarrel_.end() &&
0892          std::find(badPixelDetsEndcap_.begin(), badPixelDetsEndcap_.end(), det) == badPixelDetsEndcap_.end();
0893 }
0894 PixelInactiveAreaFinder::DetGroup PixelInactiveAreaFinder::badAdjecentDetsBarrel(const det_t& det) {
0895   using std::remove_if;
0896
0897   DetGroup adj;
0898   auto const tTopo = trackerTopology_;
0899   auto const& detId = DetId(det);
0900   unsigned int layer = tTopo->pxbLayer(detId);
0901   unsigned int ladder = tTopo->pxbLadder(detId);
0902   unsigned int module = tTopo->pxbModule(detId);
0903   unsigned int nLads = nBPixLadders[layer];
0904   //add detectors from next and previous ladder
0905   adj.push_back(tTopo->pxbDetId(layer, ((ladder - 1) + 1) % nLads + 1, module)());
0906   adj.push_back(tTopo->pxbDetId(layer, ((ladder - 1) - 1 + nLads) % nLads + 1, module)());
0907   //add adjecent detectors from same ladder
0908   if (module == 1) {
0909     adj.push_back(tTopo->pxbDetId(layer, ladder, module + 1)());
0910   } else if (module == nModulesPerLadder) {
0911     adj.push_back(tTopo->pxbDetId(layer, ladder, module - 1)());
0912   } else {
0913     adj.push_back(tTopo->pxbDetId(layer, ladder, module + 1)());
0914     adj.push_back(tTopo->pxbDetId(layer, ladder, module - 1)());
0915   }
0916   //remove working detectors from list
0917   adj.erase(remove_if(adj.begin(), adj.end(), [&](auto c) { return this->detWorks(c); }), adj.end());
0918   return adj;
0919 }
0920 PixelInactiveAreaFinder::DetGroup PixelInactiveAreaFinder::badAdjecentDetsEndcap(const det_t& det) {
0921   // this might be faster if adjecent
0922   using std::ignore;
0923   using std::tie;
0924   DetGroup adj;
0925   Span_t phiSpan, phiSpanComp;
0926   float z, zComp;
0927   unsigned int disk, diskComp;
0928   auto const& detSurf = trackerGeometry_->idToDet(DetId(det))->surface();
0929   phiSpan = detSurf.phiSpan();
0930   tie(z, ignore) = detSurf.zSpan();
0931   disk = trackerTopology_->pxfDisk(DetId(det));
0932   // add detectors from same disk whose phi ranges overlap to the adjecent list
0933   for (auto const& detComp : badPixelDetsEndcap_) {
0934     auto const& detIdComp = DetId(detComp);
0935     auto const& detSurfComp = trackerGeometry_->idToDet(detIdComp)->surface();
0936     diskComp = trackerTopology_->pxfDisk(detIdComp);
0937     phiSpanComp = detSurfComp.phiSpan();
0938     tie(zComp, ignore) = detSurfComp.zSpan();
0939     if (det != detComp && disk == diskComp && z * zComp > 0 && phiRangesOverlap(phiSpan, phiSpanComp)) {
0940       adj.push_back(detComp);
0941     }
0942   }
0943   return adj;
0944 }
0945 PixelInactiveAreaFinder::DetGroup PixelInactiveAreaFinder::reachableDetGroup(const det_t& initDet,
0946                                                                              DetectorSet& foundDets) {
0947   DetGroup reachableDetGroup;
0948   std::queue<det_t> workQueue;
0949   det_t workDet;
0950   DetGroup badAdjDets;
0951   foundDets.insert(initDet);
0952   workQueue.push(initDet);
0953   reachableDetGroup.push_back(initDet);
0954   while (!workQueue.empty()) {
0955     workDet = workQueue.front();
0956     workQueue.pop();
0957     if (DetId(workDet).subdetId() == PixelSubdetector::PixelBarrel) {
0958       badAdjDets = this->badAdjecentDetsBarrel(workDet);
0959     } else if (DetId(workDet).subdetId() == PixelSubdetector::PixelEndcap) {
0960       badAdjDets = this->badAdjecentDetsEndcap(workDet);
0961     } else {
0962       badAdjDets = {};
0963     }
0964     for (auto const& badDet : badAdjDets) {
0965       if (foundDets.find(badDet) == foundDets.end()) {
0966         reachableDetGroup.push_back(badDet);
0967         foundDets.insert(badDet);
0968         workQueue.push(badDet);
0969       }
0970     }
0971   }
0972   return reachableDetGroup;
0973 }
0974 PixelInactiveAreaFinder::DetGroupContainer PixelInactiveAreaFinder::badDetGroupsBarrel() {
0975   DetGroupContainer detGroups;
0976   DetectorSet foundDets;
0977   for (auto const& badDet : badPixelDetsBarrel_) {
0978     if (foundDets.find(badDet) == foundDets.end()) {
0979       detGroups.push_back(this->reachableDetGroup(badDet, foundDets));
0980     }
0981   }
0982   return detGroups;
0983 }
0984 PixelInactiveAreaFinder::DetGroupContainer PixelInactiveAreaFinder::badDetGroupsEndcap() {
0985   DetGroupContainer detGroups;
0986   DetectorSet foundDets;
0987   for (auto const& badDet : badPixelDetsEndcap_) {
0988     if (foundDets.find(badDet) == foundDets.end()) {
0989       auto adjacentDets = this->reachableDetGroup(badDet, foundDets);
0990       if (ignoreSingleFPixPanelModules_ && adjacentDets.size() == 1) {
0991         // size==1 means that only a single panel of a blade was inactive
0992         // because of the large overlap with the other panel (i.e.
0993         // redundancy in the detectory) ignoring these may help to decrease fakes
0994         continue;
0995       }
0996       detGroups.push_back(adjacentDets);
0997     }
0998   }
0999   return detGroups;
1000 }
1001 // Functions for finding DetGroupSpans
1002 void PixelInactiveAreaFinder::getPhiSpanBarrel(const DetGroup& detGroup, DetGroupSpan& cspan) {
1003   // find phiSpan using ordered vector of unique ladders in detGroup
1004   if (detGroup.empty()) {
1005     cspan = DetGroupSpan();
1006     return;
1007   } else {
1008     cspan.layer = trackerTopology_->pxbLayer(DetId(detGroup[0]));
1009     cspan.disk = 0;
1010   }
1011   using uint = unsigned int;
1012   using LadderSet = std::set<uint>;
1013   using LadVec = std::vector<uint>;
1014   LadderSet lads;
1015   for (auto const& det : detGroup) {
1016     lads.insert(trackerTopology_->pxbLadder(DetId(det)));
1017   }
1018   LadVec ladv(lads.begin(), lads.end());
1019   uint nLadders = nBPixLadders[cspan.layer];
1020   // find start ladder of detGroup
1021   uint i = 0;
1022   uint currentLadder = ladv[0];
1023   uint previousLadder = ladv[(ladv.size() + i - 1) % ladv.size()];
1024   // loop until discontinuity is found from vector
1025   while ((nLadders + currentLadder - 1) % nLadders == previousLadder) {
1026     ++i;
1027     currentLadder = ladv[i % ladv.size()];
1028     previousLadder = ladv[(ladv.size() + i - 1) % ladv.size()];
1029     if (i == ladv.size()) {
1030       cspan.phiSpan.first = std::numeric_limits<float>::epsilon();
1031       cspan.phiSpan.second = -std::numeric_limits<float>::epsilon();
1032       return;
1033     }
1034   }
1035   uint startLadder = currentLadder;
1036   uint endLadder = previousLadder;
1037   auto detStart = trackerTopology_->pxbDetId(cspan.layer, startLadder, 1);
1038   auto detEnd = trackerTopology_->pxbDetId(cspan.layer, endLadder, 1);
1039   cspan.phiSpan.first = trackerGeometry_->idToDet(detStart)->surface().phiSpan().first;
1040   cspan.phiSpan.second = trackerGeometry_->idToDet(detEnd)->surface().phiSpan().second;
1041 }
1042 void PixelInactiveAreaFinder::getPhiSpanEndcap(const DetGroup& detGroup, DetGroupSpan& cspan) {
1043   // this is quite naive/bruteforce method
1044   // 1) it starts by taking one detector from detGroup and starts to compare it to others
1045   // 2) when it finds overlapping detector in clockwise direction it starts comparing
1046   //    found detector to others
1047   // 3) search stops until no overlapping detectors in clockwise detector or all detectors
1048   //    have been work detector
1049   Stream ss;
1050   bool found = false;
1051   auto const tGeom = trackerGeometry_;
1052   DetGroup::const_iterator startDetIter = detGroup.begin();
1053   Span_t phiSpan, phiSpanComp;
1054   unsigned int counter = 0;
1055   while (!found) {
1056     phiSpan = tGeom->idToDet(DetId(*startDetIter))->surface().phiSpan();
1057     for (DetGroup::const_iterator compDetIter = detGroup.begin(); compDetIter != detGroup.end(); ++compDetIter) {
1058       phiSpanComp = tGeom->idToDet(DetId(*compDetIter))->surface().phiSpan();
1059       if (phiRangesOverlap(phiSpan, phiSpanComp) && phiMoreClockwise(phiSpanComp.first, phiSpan.first) &&
1060           startDetIter != compDetIter) {
1061         ++counter;
1062         if (counter > detGroup.size()) {
1063           cspan.phiSpan.first = std::numeric_limits<float>::epsilon();
1064           cspan.phiSpan.second = -std::numeric_limits<float>::epsilon();
1065           return;
1066         }
1067         startDetIter = compDetIter;
1068         break;
1069       } else if (compDetIter == detGroup.end() - 1) {
1070         found = true;
1071       }
1072     }
1073   }
1074   cspan.phiSpan.first = phiSpan.first;
1075   // second with same method}
1076   found = false;
1077   DetGroup::const_iterator endDetIter = detGroup.begin();
1078   counter = 0;
1079   while (!found) {
1080     phiSpan = tGeom->idToDet(DetId(*endDetIter))->surface().phiSpan();
1081     for (DetGroup::const_iterator compDetIter = detGroup.begin(); compDetIter != detGroup.end(); ++compDetIter) {
1082       phiSpanComp = tGeom->idToDet(DetId(*compDetIter))->surface().phiSpan();
1083       if (phiRangesOverlap(phiSpan, phiSpanComp) && phiMoreCounterclockwise(phiSpanComp.second, phiSpan.second) &&
1084           endDetIter != compDetIter) {
1085         ++counter;
1086         if (counter > detGroup.size()) {
1087           cspan.phiSpan.first = std::numeric_limits<float>::epsilon();
1088           cspan.phiSpan.second = -std::numeric_limits<float>::epsilon();
1089           return;
1090         }
1091         endDetIter = compDetIter;
1092         break;
1093       } else if (compDetIter == detGroup.end() - 1) {
1094         found = true;
1095       }
1096     }
1097   }
1098   cspan.phiSpan.second = phiSpan.second;
1099 }
1100 void PixelInactiveAreaFinder::getZSpan(const DetGroup& detGroup, DetGroupSpan& cspan) {
1101   auto cmpFun = [this](det_t detA, det_t detB) {
1102     return trackerGeometry_->idToDet(DetId(detA))->surface().zSpan().first <
1103            trackerGeometry_->idToDet(DetId(detB))->surface().zSpan().first;
1104   };
1105
1106   auto minmaxIters = std::minmax_element(detGroup.begin(), detGroup.end(), cmpFun);
1107   cspan.zSpan.first = trackerGeometry_->idToDet(DetId(*(minmaxIters.first)))->surface().zSpan().first;
1108   cspan.zSpan.second = trackerGeometry_->idToDet(DetId(*(minmaxIters.second)))->surface().zSpan().second;
1109 }
1110 void PixelInactiveAreaFinder::getRSpan(const DetGroup& detGroup, DetGroupSpan& cspan) {
1111   auto cmpFun = [this](det_t detA, det_t detB) {
1112     return trackerGeometry_->idToDet(DetId(detA))->surface().rSpan().first <
1113            trackerGeometry_->idToDet(DetId(detB))->surface().rSpan().first;
1114   };
1115
1116   auto minmaxIters = std::minmax_element(detGroup.begin(), detGroup.end(), cmpFun);
1117   cspan.rSpan.first = trackerGeometry_->idToDet(DetId(*(minmaxIters.first)))->surface().rSpan().first;
1118   cspan.rSpan.second = trackerGeometry_->idToDet(DetId(*(minmaxIters.second)))->surface().rSpan().second;
1119 }
1120 void PixelInactiveAreaFinder::getSpan(const DetGroup& detGroup, DetGroupSpan& cspan) {
1121   auto firstDetIt = detGroup.begin();
1122   if (firstDetIt != detGroup.end()) {
1123     cspan.subdetId = DetId(*firstDetIt).subdetId();
1124     if (cspan.subdetId == 1) {
1125       cspan.layer = trackerTopology_->pxbLayer(DetId(*firstDetIt));
1126       cspan.disk = 0;
1127       getPhiSpanBarrel(detGroup, cspan);
1128     } else if (cspan.subdetId == 2) {
1129       cspan.disk = trackerTopology_->pxfDisk(DetId(*firstDetIt));
1130       cspan.layer = 0;
1131       getPhiSpanEndcap(detGroup, cspan);
1132     }
1133     getZSpan(detGroup, cspan);
1134     getRSpan(detGroup, cspan);
1135   }
1136 }
1137 PixelInactiveAreaFinder::DetGroupSpanContainerPair PixelInactiveAreaFinder::detGroupSpans() {
1138   DetGroupSpanContainer cspansBarrel;
1139   DetGroupSpanContainer cspansEndcap;
1140   DetGroupContainer badDetGroupsBar = badDetGroupsBarrel();
1141   DetGroupContainer badDetGroupsEnd = badDetGroupsEndcap();
1142   for (auto const& detGroup : badDetGroupsBar) {
1143     DetGroupSpan cspan;
1144     getSpan(detGroup, cspan);
1145     cspansBarrel.push_back(cspan);
1146   }
1147   for (auto const& detGroup : badDetGroupsEnd) {
1148     DetGroupSpan cspan;
1149     getSpan(detGroup, cspan);
1150     cspansEndcap.push_back(cspan);
1151   }
1152   return DetGroupSpanContainerPair(cspansBarrel, cspansEndcap);
1153 }