Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:23:04 

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   getBadPixelDets(iEvent, iSetup);
0563 
0564   //write files for plotting
0565   if (createPlottingFiles_) {
0566     createPlottingFiles();
0567   }
0568 
0569   // find detGroupSpans ie phi,r,z limits for detector detGroups that are not working
0570   // returns pair where first is barrel spans and second endcap spans
0571   DetGroupSpanContainerPair cspans = detGroupSpans();
0572 
0573   // map spans to a vector with consistent indexing with inactiveLayers_ and inactiveLayerSetIndices_
0574   // TODO: try to move the inner logic towards this direction as well
0575   std::vector<DetGroupSpanContainer> spans(inactiveLayers_.size());
0576 
0577   auto doWork = [&](const DetGroupSpanContainer& container) {
0578     for (const auto& span : container) {
0579       const auto subdet = span.subdetId == PixelSubdetector::PixelBarrel ? GeomDetEnumerators::PixelBarrel
0580                                                                          : GeomDetEnumerators::PixelEndcap;
0581       const auto side = (subdet == GeomDetEnumerators::PixelBarrel
0582                              ? TrackerDetSide::Barrel
0583                              : (span.zSpan.first < 0 ? TrackerDetSide::NegEndcap : TrackerDetSide::PosEndcap));
0584       const auto layer = subdet == GeomDetEnumerators::PixelBarrel ? span.layer : span.disk;
0585       auto found = std::find(inactiveLayers_.begin(), inactiveLayers_.end(), SeedingLayerId(subdet, side, layer));
0586       if (found != inactiveLayers_.end()) {  // it is possible that this layer is ignored by the configuration
0587         spans[std::distance(inactiveLayers_.begin(), found)].push_back(span);
0588       }
0589     }
0590   };
0591   doWork(cspans.first);
0592   doWork(cspans.second);
0593 
0594   auto ret =
0595       InactiveAreas(&inactiveLayers_, std::move(spans), &inactiveLayerSetIndices_, &layerSetIndexInactiveToActive_);
0596 
0597   if (debug_) {
0598     printOverlapSpans(ret);
0599   }
0600 
0601   return ret;
0602 }
0603 
0604 // Functions for fetching date from handles
0605 void PixelInactiveAreaFinder::updatePixelDets(const edm::EventSetup& iSetup) {
0606   if (!geometryWatcher_.check(iSetup))
0607     return;
0608 
0609   // deduce the number of ladders per layer and the number of modules per ladder
0610   {
0611     // sanity checks
0612     if (trackerGeometry_->numberOfLayers(PixelSubdetector::PixelBarrel) != 4) {
0613       throw cms::Exception("NotImplemented")
0614           << "This module supports only a detector with 4 pixel barrel layers, the current geometry has "
0615           << trackerGeometry_->numberOfLayers(PixelSubdetector::PixelBarrel);
0616     }
0617     if (trackerGeometry_->numberOfLayers(PixelSubdetector::PixelEndcap) != 3) {
0618       throw cms::Exception("NotImplemented")
0619           << "This module supports only a detector with 3 pixel forward disks, the current geometry has "
0620           << trackerGeometry_->numberOfLayers(PixelSubdetector::PixelEndcap);
0621     }
0622 
0623     std::array<std::array<unsigned short, 100>, 4> counts = {};  // assume at most 100 ladders per layer
0624     for (const auto& det : trackerGeometry_->detsPXB()) {
0625       const auto layer = trackerTopology_->layer(det->geographicalId());
0626       const auto ladder = trackerTopology_->pxbLadder(det->geographicalId());
0627       if (ladder >= 100) {
0628         throw cms::Exception("LogicError")
0629             << "Got a ladder with number " << ladder
0630             << " while the expected maximum was 100; either something is wrong or the maximum has to be increased.";
0631       }
0632       counts[layer - 1][ladder - 1] += 1;  // numbering of layer and ladder starts at 1
0633     }
0634 
0635     // take number of modules per ladder from the first ladder of the first layer (such better exist)
0636     // other ladders better have the same number
0637     nModulesPerLadder = counts[0][0];
0638     if (nModulesPerLadder == 0) {
0639       throw cms::Exception("LogicError") << "Ladder 1 of layer 1 has 0 modules, something fishy is going on.";
0640     }
0641 
0642     LogDebug("PixelInactiveAreaFinder") << "Number of modules per ladder " << nModulesPerLadder
0643                                         << "; below are number of ladders per layer";
0644 
0645     // number of ladders
0646     for (unsigned layer = 0; layer < 4; ++layer) {
0647       nBPixLadders[layer] =
0648           std::count_if(counts[layer].begin(), counts[layer].end(), [](unsigned short val) { return val > 0; });
0649       LogTrace("PixelInactiveAreaFinder")
0650           << "BPix layer " << (layer + 1) << " has " << nBPixLadders[layer] << " ladders";
0651 
0652       auto fail = std::count_if(counts[layer].begin(), counts[layer].end(), [&](unsigned short val) {
0653         return val != nModulesPerLadder && val > 0;
0654       });
0655       if (fail != 0) {
0656         throw cms::Exception("LogicError")
0657             << "Layer " << (layer + 1) << " had " << fail
0658             << " ladders whose number of modules/ladder differed from the ladder 1 of layer 1 (" << nModulesPerLadder
0659             << "). Something fishy is going on.";
0660       }
0661     }
0662   }
0663 
0664   // don't bother with the rest if not needed
0665   if (!createPlottingFiles_)
0666     return;
0667 
0668   pixelDetsBarrel_.clear();
0669   pixelDetsEndcap_.clear();
0670 
0671   for (auto const& geomDetPtr : trackerGeometry_->detsPXB()) {
0672     if (geomDetPtr->geographicalId().subdetId() == PixelSubdetector::PixelBarrel) {
0673       pixelDetsBarrel_.push_back(geomDetPtr->geographicalId().rawId());
0674     }
0675   }
0676   for (auto const& geomDetPtr : trackerGeometry_->detsPXF()) {
0677     if (geomDetPtr->geographicalId().subdetId() == PixelSubdetector::PixelEndcap) {
0678       pixelDetsEndcap_.push_back(geomDetPtr->geographicalId().rawId());
0679     }
0680   }
0681   std::sort(pixelDetsBarrel_.begin(), pixelDetsBarrel_.end());
0682   std::sort(pixelDetsEndcap_.begin(), pixelDetsEndcap_.end());
0683 }
0684 void PixelInactiveAreaFinder::getBadPixelDets(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0685   auto addDetId = [&](const auto id) {
0686     const auto detid = DetId(id);
0687     const auto subdet = detid.subdetId();
0688     if (subdet == PixelSubdetector::PixelBarrel) {
0689       badPixelDetsBarrel_.push_back(detid.rawId());
0690     } else if (subdet == PixelSubdetector::PixelEndcap) {
0691       badPixelDetsEndcap_.push_back(detid.rawId());
0692     }
0693   };
0694 
0695   // SiPixelQuality
0696   auto const& pixelQuality = iSetup.getData(pixelQualityToken_);
0697 
0698   for (auto const& disabledModule : pixelQuality.getBadComponentList()) {
0699     addDetId(disabledModule.DetID);
0700   }
0701 
0702   // dynamic bad modules
0703   for (const auto& token : inactivePixelDetectorTokens_) {
0704     edm::Handle<DetIdCollection> detIds;
0705     iEvent.getByToken(token, detIds);
0706     for (const auto& id : *detIds) {
0707       addDetId(id);
0708     }
0709   }
0710 
0711   // dynamic bad ROCs ("Fed25")
0712   // TODO: consider moving to finer-grained areas for inactive ROCs
0713   for (const auto& token : badPixelFEDChannelsTokens_) {
0714     edm::Handle<PixelFEDChannelCollection> pixelFEDChannelCollectionHandle;
0715     iEvent.getByToken(token, pixelFEDChannelCollectionHandle);
0716     for (const auto& disabledChannels : *pixelFEDChannelCollectionHandle) {
0717       addDetId(disabledChannels.detId());
0718     }
0719   }
0720 
0721   // remove duplicates
0722   std::sort(badPixelDetsBarrel_.begin(), badPixelDetsBarrel_.end());
0723   std::sort(badPixelDetsEndcap_.begin(), badPixelDetsEndcap_.end());
0724   badPixelDetsBarrel_.erase(std::unique(badPixelDetsBarrel_.begin(), badPixelDetsBarrel_.end()),
0725                             badPixelDetsBarrel_.end());
0726   badPixelDetsEndcap_.erase(std::unique(badPixelDetsEndcap_.begin(), badPixelDetsEndcap_.end()),
0727                             badPixelDetsEndcap_.end());
0728 }
0729 // Printing functions
0730 void PixelInactiveAreaFinder::detInfo(const det_t& det, Stream& ss) {
0731   using std::fixed;
0732   using std::left;
0733   using std::noshowpos;
0734   using std::right;
0735   using std::setfill;
0736   using std::setprecision;
0737   using std::setw;
0738   using std::showpos;
0739   using std::tie;
0740   std::string deli = "; ";
0741   ss << "id:[" << det << "]" << deli;
0742   ss << "subdetid:[" << DetId(det).subdetId() << "]" << deli;
0743   if (DetId(det).subdetId() == PixelSubdetector::PixelBarrel) {
0744     unsigned int layer = trackerTopology_->pxbLayer(DetId(det));
0745     unsigned int ladder = trackerTopology_->pxbLadder(DetId(det));
0746     unsigned int module = trackerTopology_->pxbModule(DetId(det));
0747     ss << "layer:[" << layer << "]" << deli << "ladder:[" << right << setw(2) << ladder << "]" << deli << "module:["
0748        << module << "]" << deli;
0749   } else if (DetId(det).subdetId() == PixelSubdetector::PixelEndcap) {
0750     unsigned int disk = trackerTopology_->pxfDisk(DetId(det));
0751     unsigned int blade = trackerTopology_->pxfBlade(DetId(det));
0752     unsigned int panel = trackerTopology_->pxfPanel(DetId(det));
0753     ss << left << setw(6) << "disk:"
0754        << "[" << right << disk << "]" << deli << left << setw(7) << "blade:"
0755        << "[" << setw(2) << right << blade << "]" << deli << left << setw(7) << "panel:"
0756        << "[" << right << panel << "]" << deli;
0757   }
0758   float phiA, phiB, zA, zB, rA, rB;
0759   auto detSurface = trackerGeometry_->idToDet(DetId(det))->surface();
0760   tie(phiA, phiB) = detSurface.phiSpan();
0761   tie(zA, zB) = detSurface.zSpan();
0762   tie(rA, rB) = detSurface.rSpan();
0763   ss << setprecision(16) << fixed << showpos << setfill(' ') << "phi:[" << right << setw(12) << phiA << "," << left
0764      << setw(12) << phiB << "]" << deli << "z:[" << right << setw(7) << zA << "," << left << setw(7) << zB << "]"
0765      << deli << noshowpos << "r:[" << right << setw(10) << rA << "," << left << setw(10) << rB << "]" << deli;
0766 }
0767 void PixelInactiveAreaFinder::printPixelDets() {
0768   edm::LogPrint("PixelInactiveAreaFinder") << "Barrel detectors:";
0769   Stream ss;
0770   for (auto const& det : pixelDetsBarrel_) {
0771     detInfo(det, ss);
0772     edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0773     ss.str(std::string());
0774   }
0775   edm::LogPrint("PixelInactiveAreaFinder") << "Endcap detectors;";
0776   for (auto const& det : pixelDetsEndcap_) {
0777     detInfo(det, ss);
0778     edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0779     ss.str(std::string());
0780   }
0781 }
0782 void PixelInactiveAreaFinder::printBadPixelDets() {
0783   edm::LogPrint("PixelInactiveAreaFinder") << "Bad barrel detectors:";
0784   Stream ss;
0785   for (auto const& det : badPixelDetsBarrel_) {
0786     detInfo(det, ss);
0787     edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0788     ss.str(std::string());
0789   }
0790   edm::LogPrint("PixelInactiveAreaFinder") << "Endcap detectors;";
0791   for (auto const& det : badPixelDetsEndcap_) {
0792     detInfo(det, ss);
0793     edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0794     ss.str(std::string());
0795   }
0796 }
0797 void PixelInactiveAreaFinder::printBadDetGroups() {
0798   DetGroupContainer badDetGroupsBar = badDetGroupsBarrel();
0799   DetGroupContainer badDetGroupsEnd = badDetGroupsEndcap();
0800   Stream ss;
0801   for (auto const& detGroup : badDetGroupsBar) {
0802     ss << std::setfill(' ') << std::left << std::setw(16) << "DetGroup:";
0803     DetGroupSpan cspan;
0804     getPhiSpanBarrel(detGroup, cspan);
0805     getZSpan(detGroup, cspan);
0806     getRSpan(detGroup, cspan);
0807     detGroupSpanInfo(cspan, ss);
0808     ss << std::endl;
0809     for (auto const& det : detGroup) {
0810       detInfo(det, ss);
0811       ss << std::endl;
0812     }
0813     ss << std::endl;
0814   }
0815   for (auto const& detGroup : badDetGroupsEnd) {
0816     ss << std::setfill(' ') << std::left << std::setw(16) << "DetGroup:";
0817     DetGroupSpan cspan;
0818     getPhiSpanEndcap(detGroup, cspan);
0819     getZSpan(detGroup, cspan);
0820     getRSpan(detGroup, cspan);
0821     detGroupSpanInfo(cspan, ss);
0822     ss << std::endl;
0823     for (auto const& det : detGroup) {
0824       detInfo(det, ss);
0825       ss << std::endl;
0826     }
0827     ss << std::endl;
0828   }
0829   edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0830 }
0831 void PixelInactiveAreaFinder::printBadDetGroupSpans() {
0832   DetGroupSpanContainerPair cspans = detGroupSpans();
0833   Stream ss;
0834   for (auto const& cspan : cspans.first) {
0835     detGroupSpanInfo(cspan, ss);
0836     ss << std::endl;
0837   }
0838   for (auto const& cspan : cspans.second) {
0839     detGroupSpanInfo(cspan, ss);
0840     ss << std::endl;
0841   }
0842   edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
0843 }
0844 void PixelInactiveAreaFinder::createPlottingFiles() {
0845   // All detectors to file DETECTORS
0846   Stream ss;
0847   std::ofstream fsDet("DETECTORS.txt");
0848   for (auto const& det : pixelDetsBarrel_) {
0849     detInfo(det, ss);
0850     ss << std::endl;
0851   }
0852   edm::LogPrint("PixelInactiveAreaFinder") << "Endcap detectors;";
0853   for (auto const& det : pixelDetsEndcap_) {
0854     detInfo(det, ss);
0855     ss << std::endl;
0856   }
0857   fsDet << ss.rdbuf();
0858   ss.str(std::string());
0859   // Bad detectors
0860   std::ofstream fsBadDet("BADDETECTORS.txt");
0861   for (auto const& det : badPixelDetsBarrel_) {
0862     detInfo(det, ss);
0863     ss << std::endl;
0864   }
0865   for (auto const& det : badPixelDetsEndcap_) {
0866     detInfo(det, ss);
0867     ss << std::endl;
0868   }
0869   fsBadDet << ss.rdbuf();
0870   ss.str(std::string());
0871   // detgroupspans
0872   std::ofstream fsSpans("DETGROUPSPANS.txt");
0873   DetGroupSpanContainerPair cspans = detGroupSpans();
0874   for (auto const& cspan : cspans.first) {
0875     detGroupSpanInfo(cspan, ss);
0876     ss << std::endl;
0877   }
0878   for (auto const& cspan : cspans.second) {
0879     detGroupSpanInfo(cspan, ss);
0880     ss << std::endl;
0881   }
0882   fsSpans << ss.rdbuf();
0883   ss.str(std::string());
0884 }
0885 // Functions for finding bad detGroups
0886 bool PixelInactiveAreaFinder::detWorks(det_t det) {
0887   return std::find(badPixelDetsBarrel_.begin(), badPixelDetsBarrel_.end(), det) == badPixelDetsBarrel_.end() &&
0888          std::find(badPixelDetsEndcap_.begin(), badPixelDetsEndcap_.end(), det) == badPixelDetsEndcap_.end();
0889 }
0890 PixelInactiveAreaFinder::DetGroup PixelInactiveAreaFinder::badAdjecentDetsBarrel(const det_t& det) {
0891   using std::remove_if;
0892 
0893   DetGroup adj;
0894   auto const tTopo = trackerTopology_;
0895   auto const& detId = DetId(det);
0896   unsigned int layer = tTopo->pxbLayer(detId);
0897   unsigned int ladder = tTopo->pxbLadder(detId);
0898   unsigned int module = tTopo->pxbModule(detId);
0899   unsigned int nLads = nBPixLadders[layer];
0900   //add detectors from next and previous ladder
0901   adj.push_back(tTopo->pxbDetId(layer, ((ladder - 1) + 1) % nLads + 1, module)());
0902   adj.push_back(tTopo->pxbDetId(layer, ((ladder - 1) - 1 + nLads) % nLads + 1, module)());
0903   //add adjecent detectors from same ladder
0904   if (module == 1) {
0905     adj.push_back(tTopo->pxbDetId(layer, ladder, module + 1)());
0906   } else if (module == nModulesPerLadder) {
0907     adj.push_back(tTopo->pxbDetId(layer, ladder, module - 1)());
0908   } else {
0909     adj.push_back(tTopo->pxbDetId(layer, ladder, module + 1)());
0910     adj.push_back(tTopo->pxbDetId(layer, ladder, module - 1)());
0911   }
0912   //remove working detectors from list
0913   adj.erase(remove_if(adj.begin(), adj.end(), [&](auto c) { return this->detWorks(c); }), adj.end());
0914   return adj;
0915 }
0916 PixelInactiveAreaFinder::DetGroup PixelInactiveAreaFinder::badAdjecentDetsEndcap(const det_t& det) {
0917   // this might be faster if adjecent
0918   using std::ignore;
0919   using std::tie;
0920   DetGroup adj;
0921   Span_t phiSpan, phiSpanComp;
0922   float z, zComp;
0923   unsigned int disk, diskComp;
0924   auto const& detSurf = trackerGeometry_->idToDet(DetId(det))->surface();
0925   phiSpan = detSurf.phiSpan();
0926   tie(z, ignore) = detSurf.zSpan();
0927   disk = trackerTopology_->pxfDisk(DetId(det));
0928   // add detectors from same disk whose phi ranges overlap to the adjecent list
0929   for (auto const& detComp : badPixelDetsEndcap_) {
0930     auto const& detIdComp = DetId(detComp);
0931     auto const& detSurfComp = trackerGeometry_->idToDet(detIdComp)->surface();
0932     diskComp = trackerTopology_->pxfDisk(detIdComp);
0933     phiSpanComp = detSurfComp.phiSpan();
0934     tie(zComp, ignore) = detSurfComp.zSpan();
0935     if (det != detComp && disk == diskComp && z * zComp > 0 && phiRangesOverlap(phiSpan, phiSpanComp)) {
0936       adj.push_back(detComp);
0937     }
0938   }
0939   return adj;
0940 }
0941 PixelInactiveAreaFinder::DetGroup PixelInactiveAreaFinder::reachableDetGroup(const det_t& initDet,
0942                                                                              DetectorSet& foundDets) {
0943   DetGroup reachableDetGroup;
0944   std::queue<det_t> workQueue;
0945   det_t workDet;
0946   DetGroup badAdjDets;
0947   foundDets.insert(initDet);
0948   workQueue.push(initDet);
0949   reachableDetGroup.push_back(initDet);
0950   while (!workQueue.empty()) {
0951     workDet = workQueue.front();
0952     workQueue.pop();
0953     if (DetId(workDet).subdetId() == PixelSubdetector::PixelBarrel) {
0954       badAdjDets = this->badAdjecentDetsBarrel(workDet);
0955     } else if (DetId(workDet).subdetId() == PixelSubdetector::PixelEndcap) {
0956       badAdjDets = this->badAdjecentDetsEndcap(workDet);
0957     } else {
0958       badAdjDets = {};
0959     }
0960     for (auto const& badDet : badAdjDets) {
0961       if (foundDets.find(badDet) == foundDets.end()) {
0962         reachableDetGroup.push_back(badDet);
0963         foundDets.insert(badDet);
0964         workQueue.push(badDet);
0965       }
0966     }
0967   }
0968   return reachableDetGroup;
0969 }
0970 PixelInactiveAreaFinder::DetGroupContainer PixelInactiveAreaFinder::badDetGroupsBarrel() {
0971   DetGroupContainer detGroups;
0972   DetectorSet foundDets;
0973   for (auto const& badDet : badPixelDetsBarrel_) {
0974     if (foundDets.find(badDet) == foundDets.end()) {
0975       detGroups.push_back(this->reachableDetGroup(badDet, foundDets));
0976     }
0977   }
0978   return detGroups;
0979 }
0980 PixelInactiveAreaFinder::DetGroupContainer PixelInactiveAreaFinder::badDetGroupsEndcap() {
0981   DetGroupContainer detGroups;
0982   DetectorSet foundDets;
0983   for (auto const& badDet : badPixelDetsEndcap_) {
0984     if (foundDets.find(badDet) == foundDets.end()) {
0985       auto adjacentDets = this->reachableDetGroup(badDet, foundDets);
0986       if (ignoreSingleFPixPanelModules_ && adjacentDets.size() == 1) {
0987         // size==1 means that only a single panel of a blade was inactive
0988         // because of the large overlap with the other panel (i.e.
0989         // redundancy in the detectory) ignoring these may help to decrease fakes
0990         continue;
0991       }
0992       detGroups.push_back(adjacentDets);
0993     }
0994   }
0995   return detGroups;
0996 }
0997 // Functions for finding DetGroupSpans
0998 void PixelInactiveAreaFinder::getPhiSpanBarrel(const DetGroup& detGroup, DetGroupSpan& cspan) {
0999   // find phiSpan using ordered vector of unique ladders in detGroup
1000   if (detGroup.empty()) {
1001     cspan = DetGroupSpan();
1002     return;
1003   } else {
1004     cspan.layer = trackerTopology_->pxbLayer(DetId(detGroup[0]));
1005     cspan.disk = 0;
1006   }
1007   using uint = unsigned int;
1008   using LadderSet = std::set<uint>;
1009   using LadVec = std::vector<uint>;
1010   LadderSet lads;
1011   for (auto const& det : detGroup) {
1012     lads.insert(trackerTopology_->pxbLadder(DetId(det)));
1013   }
1014   LadVec ladv(lads.begin(), lads.end());
1015   uint nLadders = nBPixLadders[cspan.layer];
1016   // find start ladder of detGroup
1017   uint i = 0;
1018   uint currentLadder = ladv[0];
1019   uint previousLadder = ladv[(ladv.size() + i - 1) % ladv.size()];
1020   // loop until discontinuity is found from vector
1021   while ((nLadders + currentLadder - 1) % nLadders == previousLadder) {
1022     ++i;
1023     currentLadder = ladv[i % ladv.size()];
1024     previousLadder = ladv[(ladv.size() + i - 1) % ladv.size()];
1025     if (i == ladv.size()) {
1026       cspan.phiSpan.first = std::numeric_limits<float>::epsilon();
1027       cspan.phiSpan.second = -std::numeric_limits<float>::epsilon();
1028       return;
1029     }
1030   }
1031   uint startLadder = currentLadder;
1032   uint endLadder = previousLadder;
1033   auto detStart = trackerTopology_->pxbDetId(cspan.layer, startLadder, 1);
1034   auto detEnd = trackerTopology_->pxbDetId(cspan.layer, endLadder, 1);
1035   cspan.phiSpan.first = trackerGeometry_->idToDet(detStart)->surface().phiSpan().first;
1036   cspan.phiSpan.second = trackerGeometry_->idToDet(detEnd)->surface().phiSpan().second;
1037 }
1038 void PixelInactiveAreaFinder::getPhiSpanEndcap(const DetGroup& detGroup, DetGroupSpan& cspan) {
1039   // this is quite naive/bruteforce method
1040   // 1) it starts by taking one detector from detGroup and starts to compare it to others
1041   // 2) when it finds overlapping detector in clockwise direction it starts comparing
1042   //    found detector to others
1043   // 3) search stops until no overlapping detectors in clockwise detector or all detectors
1044   //    have been work detector
1045   Stream ss;
1046   bool found = false;
1047   auto const tGeom = trackerGeometry_;
1048   DetGroup::const_iterator startDetIter = detGroup.begin();
1049   Span_t phiSpan, phiSpanComp;
1050   unsigned int counter = 0;
1051   while (!found) {
1052     phiSpan = tGeom->idToDet(DetId(*startDetIter))->surface().phiSpan();
1053     for (DetGroup::const_iterator compDetIter = detGroup.begin(); compDetIter != detGroup.end(); ++compDetIter) {
1054       phiSpanComp = tGeom->idToDet(DetId(*compDetIter))->surface().phiSpan();
1055       if (phiRangesOverlap(phiSpan, phiSpanComp) && phiMoreClockwise(phiSpanComp.first, phiSpan.first) &&
1056           startDetIter != compDetIter) {
1057         ++counter;
1058         if (counter > detGroup.size()) {
1059           cspan.phiSpan.first = std::numeric_limits<float>::epsilon();
1060           cspan.phiSpan.second = -std::numeric_limits<float>::epsilon();
1061           return;
1062         }
1063         startDetIter = compDetIter;
1064         break;
1065       } else if (compDetIter == detGroup.end() - 1) {
1066         found = true;
1067       }
1068     }
1069   }
1070   cspan.phiSpan.first = phiSpan.first;
1071   // second with same method}
1072   found = false;
1073   DetGroup::const_iterator endDetIter = detGroup.begin();
1074   counter = 0;
1075   while (!found) {
1076     phiSpan = tGeom->idToDet(DetId(*endDetIter))->surface().phiSpan();
1077     for (DetGroup::const_iterator compDetIter = detGroup.begin(); compDetIter != detGroup.end(); ++compDetIter) {
1078       phiSpanComp = tGeom->idToDet(DetId(*compDetIter))->surface().phiSpan();
1079       if (phiRangesOverlap(phiSpan, phiSpanComp) && phiMoreCounterclockwise(phiSpanComp.second, phiSpan.second) &&
1080           endDetIter != compDetIter) {
1081         ++counter;
1082         if (counter > detGroup.size()) {
1083           cspan.phiSpan.first = std::numeric_limits<float>::epsilon();
1084           cspan.phiSpan.second = -std::numeric_limits<float>::epsilon();
1085           return;
1086         }
1087         endDetIter = compDetIter;
1088         break;
1089       } else if (compDetIter == detGroup.end() - 1) {
1090         found = true;
1091       }
1092     }
1093   }
1094   cspan.phiSpan.second = phiSpan.second;
1095 }
1096 void PixelInactiveAreaFinder::getZSpan(const DetGroup& detGroup, DetGroupSpan& cspan) {
1097   auto cmpFun = [this](det_t detA, det_t detB) {
1098     return trackerGeometry_->idToDet(DetId(detA))->surface().zSpan().first <
1099            trackerGeometry_->idToDet(DetId(detB))->surface().zSpan().first;
1100   };
1101 
1102   auto minmaxIters = std::minmax_element(detGroup.begin(), detGroup.end(), cmpFun);
1103   cspan.zSpan.first = trackerGeometry_->idToDet(DetId(*(minmaxIters.first)))->surface().zSpan().first;
1104   cspan.zSpan.second = trackerGeometry_->idToDet(DetId(*(minmaxIters.second)))->surface().zSpan().second;
1105 }
1106 void PixelInactiveAreaFinder::getRSpan(const DetGroup& detGroup, DetGroupSpan& cspan) {
1107   auto cmpFun = [this](det_t detA, det_t detB) {
1108     return trackerGeometry_->idToDet(DetId(detA))->surface().rSpan().first <
1109            trackerGeometry_->idToDet(DetId(detB))->surface().rSpan().first;
1110   };
1111 
1112   auto minmaxIters = std::minmax_element(detGroup.begin(), detGroup.end(), cmpFun);
1113   cspan.rSpan.first = trackerGeometry_->idToDet(DetId(*(minmaxIters.first)))->surface().rSpan().first;
1114   cspan.rSpan.second = trackerGeometry_->idToDet(DetId(*(minmaxIters.second)))->surface().rSpan().second;
1115 }
1116 void PixelInactiveAreaFinder::getSpan(const DetGroup& detGroup, DetGroupSpan& cspan) {
1117   auto firstDetIt = detGroup.begin();
1118   if (firstDetIt != detGroup.end()) {
1119     cspan.subdetId = DetId(*firstDetIt).subdetId();
1120     if (cspan.subdetId == 1) {
1121       cspan.layer = trackerTopology_->pxbLayer(DetId(*firstDetIt));
1122       cspan.disk = 0;
1123       getPhiSpanBarrel(detGroup, cspan);
1124     } else if (cspan.subdetId == 2) {
1125       cspan.disk = trackerTopology_->pxfDisk(DetId(*firstDetIt));
1126       cspan.layer = 0;
1127       getPhiSpanEndcap(detGroup, cspan);
1128     }
1129     getZSpan(detGroup, cspan);
1130     getRSpan(detGroup, cspan);
1131   }
1132 }
1133 PixelInactiveAreaFinder::DetGroupSpanContainerPair PixelInactiveAreaFinder::detGroupSpans() {
1134   DetGroupSpanContainer cspansBarrel;
1135   DetGroupSpanContainer cspansEndcap;
1136   DetGroupContainer badDetGroupsBar = badDetGroupsBarrel();
1137   DetGroupContainer badDetGroupsEnd = badDetGroupsEndcap();
1138   for (auto const& detGroup : badDetGroupsBar) {
1139     DetGroupSpan cspan;
1140     getSpan(detGroup, cspan);
1141     cspansBarrel.push_back(cspan);
1142   }
1143   for (auto const& detGroup : badDetGroupsEnd) {
1144     DetGroupSpan cspan;
1145     getSpan(detGroup, cspan);
1146     cspansEndcap.push_back(cspan);
1147   }
1148   return DetGroupSpanContainerPair(cspansBarrel, cspansEndcap);
1149 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!