PixelInactiveAreaFinder.cc

Go to the documentation of this file.

#include "PixelInactiveAreaFinder.h"

#include "FWCore/Utilities/interface/VecArray.h"
#include "FWCore/Utilities/interface/transform.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "CondFormats/DataRecord/interface/SiPixelQualityRcd.h"
#include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "DataFormats/Math/interface/deltaPhi.h"

#include "TrackingTools/TransientTrackingRecHit/interface/SeedingLayerSetsLooper.h"

#include <fstream>
#include <queue>
#include <algorithm>

std::ostream& operator<<(std::ostream& os, SeedingLayerSetsBuilder::SeedingLayerId layer) {
  if(std::get<0>(layer) == GeomDetEnumerators::PixelBarrel) {
    os << "BPix";
  }
  else {
    os << "FPix";
  }
  os << std::get<2>(layer);
  if(std::get<1>(layer) == TrackerDetSide::PosEndcap) {
    os << "_pos";
  }
  else if(std::get<1>(layer) == TrackerDetSide::NegEndcap) {
    os << "_neg";
  }
  return os;
}

namespace {
  using LayerPair = std::pair<SeedingLayerSetsBuilder::SeedingLayerId, SeedingLayerSetsBuilder::SeedingLayerId>;
  using ActiveLayerSetToInactiveSetsMap = std::map<LayerPair, edm::VecArray<LayerPair, 5> >;
  using Stream = std::stringstream;
  using Span_t = std::pair<float,float>;

  ActiveLayerSetToInactiveSetsMap createActiveToInactiveMap() {
    ActiveLayerSetToInactiveSetsMap map;

    auto bpix = [](int layer) {
      return SeedingLayerSetsBuilder::SeedingLayerId(GeomDetEnumerators::PixelBarrel, TrackerDetSide::Barrel, layer);
    };
    auto fpix_pos = [](int disk) {
      return SeedingLayerSetsBuilder::SeedingLayerId(GeomDetEnumerators::PixelEndcap, TrackerDetSide::PosEndcap, disk);
    };
    auto fpix_neg = [](int disk) {
      return SeedingLayerSetsBuilder::SeedingLayerId(GeomDetEnumerators::PixelEndcap, TrackerDetSide::NegEndcap, disk);
    };

    auto add_permutations = [&](std::array<SeedingLayerSetsBuilder::SeedingLayerId, 4> quads) {
      do {
        // skip permutations like BPix2+BPix1 or FPix1+BPix1
        // operator> works automatically
        if(quads[0] > quads[1] || quads[2] > quads[3]) continue;

        map[std::make_pair(quads[0], quads[1])].emplace_back(quads[2], quads[3]);
      } while(std::next_permutation(quads.begin(), quads.end()));
    };

    // 4 barrel
    add_permutations({{bpix(1), bpix(2), bpix(3), bpix(4)}});

    // 3 barrel, 1 forward
    add_permutations({{bpix(1), bpix(2), bpix(3), fpix_pos(1)}});
    add_permutations({{bpix(1), bpix(2), bpix(3), fpix_neg(1)}});

    // 2 barrel, 2 forward
    add_permutations({{bpix(1), bpix(2), fpix_pos(1), fpix_pos(2)}});
    add_permutations({{bpix(1), bpix(2), fpix_neg(1), fpix_neg(2)}});

    // 1 barrel, 3 forward
    add_permutations({{bpix(1), fpix_pos(1), fpix_pos(2), fpix_pos(3)}});
    add_permutations({{bpix(1), fpix_neg(1), fpix_neg(2), fpix_neg(3)}});

#ifdef EDM_ML_DEBUG
    LogDebug("PixelInactiveAreaFinder") << "Active to inactive mapping";
    for(const auto& elem: map) {
      std::stringstream ss;
      for(const auto& layerPair: elem.second) {
        ss << layerPair.first << "+" << layerPair.second << ",";
      }
      LogTrace("PixelInactiveAreaFinder") << " " << elem.first.first << "+" << elem.first.second << " => " << ss.str();
    }
#endif

    return map;
  }



  void detGroupSpanInfo(const PixelInactiveAreaFinder::DetGroupSpan & cspan, Stream & ss){
    using std::showpos;
    using std::noshowpos;
    using std::fixed;
    using std::setprecision;
    using std::setw;
    using std::setfill;
    using std::left;
    using std::right;
    std::string deli = "; ";
    ss  << "subdetid:[" << cspan.subdetId << "]" << deli;
    if(cspan.subdetId == PixelSubdetector::PixelBarrel) {
      ss << "layer:[" << cspan.layer << "]" << deli;
    }
    else {
      ss << "disk:[" << cspan.disk << "]" << deli;
    }
    ss
      //<< setfill(' ') << setw(36) << " "
      << setprecision(16)
      << showpos
      << "phi:<" << right << setw(12) << cspan.phiSpan.first << "," << left << setw(12) << cspan.phiSpan.second << ">" << deli
      << "z:<" << right << setw(7) << cspan.zSpan.first << "," << left << setw(7) << cspan.zSpan.second << ">" << deli << noshowpos
      << "r:<" << right << setw(10) << cspan.rSpan.first << "," << left << setw(10) << cspan.rSpan.second << ">" << deli
      ;
  }
  void printOverlapSpans(const PixelInactiveAreaFinder::InactiveAreas& areasLayers) {
    auto spansLayerSets = areasLayers.spansAndLayerSets(GlobalPoint(0,0,0), std::numeric_limits<float>::infinity());

    Stream ss;
    for(auto const &spansLayers: spansLayerSets) {
      ss << "Overlapping detGroups:\n";
      for(auto const cspan : spansLayers.first) {
        detGroupSpanInfo(cspan,ss);
        ss << std::endl;
      }
    }
    edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
  }

  // Functions for finding bad detGroups
  bool phiRangesOverlap(const Span_t&phiSpanA, const Span_t&phiSpanB) {
    float x1,x2,y1,y2;
    std::tie(x1,x2) = phiSpanA;
    std::tie(y1,y2) = phiSpanB;
    // assuming phi ranges are [x1,x2] and [y1,y2] and xi,yi in [-pi,pi]
    if(x1<=x2 && y1<=y2){
      return x1<=y2 && y1 <= x2;
    }else if (( x1>x2 && y1 <= y2) || (y1 > y2 && x1 <= x2 )){
      return y1 <= x2 || x1 <= y2;
    }else if (x1 > x2 && y1 > y2){
      return true;
    }else {
      return false;
    }
  }

  // Functions for finding ranges that detGroups cover
  bool phiMoreClockwise(float phiA, float phiB) {
    // return true if a is more clockwise than b
    return reco::deltaPhi(phiA, phiB) <= 0.f;
  }
  bool phiMoreCounterclockwise(float phiA, float phiB) {
    // return true if a is more counterclockwise than b
    return reco::deltaPhi(phiA, phiB) >= 0.f;
  }

  // Functions for findind overlapping functions
  float zAxisIntersection(const float zrPointA[2], const float zrPointB[2]){
    return (zrPointB[0]-zrPointA[0])/(zrPointB[1]-zrPointA[1])*(-zrPointA[1])+zrPointA[0];
  }
  bool getZAxisOverlapRangeBarrel(const PixelInactiveAreaFinder::DetGroupSpan & cspanA, const PixelInactiveAreaFinder::DetGroupSpan & cspanB, std::pair<float,float> & range){
    PixelInactiveAreaFinder::DetGroupSpan cspanUpper;
    PixelInactiveAreaFinder::DetGroupSpan cspanLower;
    if(cspanA.rSpan.second < cspanB.rSpan.first){
      cspanLower = cspanA;
      cspanUpper = cspanB;
    }else if(cspanA.rSpan.first > cspanB.rSpan.second){
      cspanUpper = cspanA;
      cspanLower = cspanB;
    }else{
      return false;
    }
    float lower = 0;
    float upper = 0;
    if(cspanUpper.zSpan.second < cspanLower.zSpan.first){
      // lower intersectionpoint, point = {z,r} in cylindrical coordinates
      const float pointUpperDetGroupL[2] = {cspanUpper.zSpan.second, cspanUpper.rSpan.second};
      const float pointLowerDetGroupL[2] = {cspanLower.zSpan.first,  cspanLower.rSpan.first};
      lower = zAxisIntersection(pointUpperDetGroupL,pointLowerDetGroupL);
      // upper intersectionpoint
      const float pointUpperDetGroupU[2] = {cspanUpper.zSpan.first, cspanUpper.rSpan.first};
      const float pointLowerDetGroupU[2] = {cspanLower.zSpan.second,  cspanLower.rSpan.second};
      upper = zAxisIntersection(pointUpperDetGroupU,pointLowerDetGroupU);
    }else if (cspanUpper.zSpan.first <= cspanLower.zSpan.second && cspanLower.zSpan.first <= cspanUpper.zSpan.second){
      // lower intersectionpoint, point = {z,r} in cylindrical coordinates
      const float pointUpperDetGroupL[2] = {cspanUpper.zSpan.second, cspanUpper.rSpan.first};
      const float pointLowerDetGroupL[2] = {cspanLower.zSpan.first,  cspanLower.rSpan.second};
      lower = zAxisIntersection(pointUpperDetGroupL,pointLowerDetGroupL);
      // upper intersectionpoint
      const float pointUpperDetGroupU[2] = {cspanUpper.zSpan.first, cspanUpper.rSpan.first};
      const float pointLowerDetGroupU[2] = {cspanLower.zSpan.second,  cspanLower.rSpan.second};
      upper = zAxisIntersection(pointUpperDetGroupU,pointLowerDetGroupU);
    }else if (cspanUpper.zSpan.first > cspanLower.zSpan.second){
      // lower intersectionpoint, point = {z,r} in cylindrical coordinates
      const float pointUpperDetGroupL[2] = {cspanUpper.zSpan.second, cspanUpper.rSpan.first};
      const float pointLowerDetGroupL[2] = {cspanLower.zSpan.first,  cspanLower.rSpan.second};
      lower = zAxisIntersection(pointUpperDetGroupL,pointLowerDetGroupL);
      // upper intersectionpoint
      const float pointUpperDetGroupU[2] = {cspanUpper.zSpan.first, cspanUpper.rSpan.second};
      const float pointLowerDetGroupU[2] = {cspanLower.zSpan.second,  cspanLower.rSpan.first};
      upper = zAxisIntersection(pointUpperDetGroupU,pointLowerDetGroupU);
    }else{
      //something wrong
      return false;
    }
    range = std::pair<float,float>(lower,upper);
    return true;
  }
  bool getZAxisOverlapRangeEndcap(const PixelInactiveAreaFinder::DetGroupSpan & cspanA, const PixelInactiveAreaFinder::DetGroupSpan & cspanB, std::pair<float,float> & range){
    // While on left hand side of pixel detector
    PixelInactiveAreaFinder::DetGroupSpan cspanNearer;
    PixelInactiveAreaFinder::DetGroupSpan cspanFurther;
    float lower = 0;
    float upper = 0;
    if(cspanA.zSpan.first < 0 && cspanB.zSpan.first < 0){
      if(cspanA.zSpan.second < cspanB.zSpan.first){
        cspanFurther = cspanA;
        cspanNearer = cspanB;
      }else if (cspanB.zSpan.second < cspanA.zSpan.first){
        cspanFurther = cspanB;
        cspanNearer = cspanA;
      }else {
#ifdef EDM_ML_DEBUG
        LogTrace("PixelInactiveAreaFinder") << "No overlap, same disk propably. Spans:";
        Stream ss;
        detGroupSpanInfo(cspanA,ss);ss<<std::endl;detGroupSpanInfo(cspanB,ss);ss<<std::endl;
        LogTrace("PixelInactiveAreaFinder") << ss.str();ss.str(std::string());
        LogTrace("PixelInactiveAreaFinder") << "**";
#endif
        return false;
      }
      if(cspanFurther.rSpan.second > cspanNearer.rSpan.first){
        const float pointA[2] = {cspanFurther.zSpan.second, cspanFurther.rSpan.second};
        const float pointB[2] = {cspanNearer.zSpan.first, cspanNearer.rSpan.first};
        lower = zAxisIntersection(pointA,pointB);
        if(cspanFurther.rSpan.first > cspanNearer.rSpan.second){
          const float pointC[2] = {cspanFurther.zSpan.first, cspanFurther.rSpan.first};
          const float pointD[2] = {cspanNearer.zSpan.second, cspanFurther.rSpan.second};
          upper = zAxisIntersection(pointC,pointD);
        }else{
          upper = std::numeric_limits<float>::infinity();
        }
      }else{
#ifdef EDM_ML_DEBUG
        LogTrace("PixelInactiveAreaFinder") << "No overlap, further detGroup is lower. Spans:";
        Stream ss;
        detGroupSpanInfo(cspanA,ss);ss<<std::endl;detGroupSpanInfo(cspanB,ss);ss<<std::endl;
        LogTrace("PixelInactiveAreaFinder") << ss.str();ss.str(std::string());
        LogTrace("PixelInactiveAreaFinder") << "**";
#endif
        return false;
      }
    }else if(cspanA.zSpan.first > 0 && cspanB.zSpan.first > 0){
      if(cspanA.zSpan.first > cspanB.zSpan.second ){
        cspanFurther = cspanA;
        cspanNearer = cspanB;
      }else if(cspanB.zSpan.first > cspanA.zSpan.second){
        cspanFurther = cspanB;
        cspanNearer = cspanA;
      }else{
#ifdef EDM_ML_DEBUG
        LogTrace("PixelInactiveAreaFinder") << "No overlap, same disk propably. Spans:";
        Stream ss;
        detGroupSpanInfo(cspanA,ss);ss<<std::endl;detGroupSpanInfo(cspanB,ss);ss<<std::endl;
        LogTrace("PixelInactiveAreaFinder") << ss.str();ss.str(std::string());
        LogTrace("PixelInactiveAreaFinder") << "**";
#endif
        return false;
      }
      if(cspanFurther.rSpan.second > cspanNearer.rSpan.first){
        const float pointA[2] = {cspanFurther.zSpan.first, cspanFurther.rSpan.second};
        const float pointB[2] = {cspanNearer.zSpan.second, cspanNearer.rSpan.first};
        upper = zAxisIntersection(pointA,pointB);
        if(cspanFurther.rSpan.first > cspanNearer.rSpan.second){
          const float pointC[2] = {cspanFurther.zSpan.second, cspanFurther.rSpan.first};
          const float pointD[2] = {cspanNearer.zSpan.first, cspanFurther.rSpan.second};
          lower = zAxisIntersection(pointC,pointD);
        }else{
          lower = -std::numeric_limits<float>::infinity();
        }
      }else{
#ifdef EDM_ML_DEBUG
        LogTrace("PixelInactiveAreaFinder") << "No overlap, further detGroup lower. Spans:";
        Stream ss;
        detGroupSpanInfo(cspanA,ss);ss<<std::endl;detGroupSpanInfo(cspanB,ss);ss<<std::endl;
        LogTrace("PixelInactiveAreaFinder") << ss.str();ss.str(std::string());
        LogTrace("PixelInactiveAreaFinder") << "**";
#endif
        return false;
      }
    }else{
#ifdef EDM_ML_DEBUG
      LogTrace("PixelInactiveAreaFinder") << "No overlap, different sides of z axis. Spans:";
      Stream ss;
      detGroupSpanInfo(cspanA,ss);ss<<std::endl;detGroupSpanInfo(cspanB,ss);ss<<std::endl;
      LogTrace("PixelInactiveAreaFinder") << ss.str();ss.str(std::string());
      LogTrace("PixelInactiveAreaFinder") << "**";
#endif
      return false;
    }
    range = std::pair<float,float>(lower,upper);
    return true;
  }

  bool getZAxisOverlapRangeBarrelEndcap(const PixelInactiveAreaFinder::DetGroupSpan & cspanBar, const PixelInactiveAreaFinder::DetGroupSpan & cspanEnd, std::pair<float,float> & range){
    float lower = 0;
    float upper = 0;
    if(cspanEnd.rSpan.second > cspanBar.rSpan.first){
      if(cspanEnd.zSpan.second < cspanBar.zSpan.first){
        // if we are on the left hand side of pixel detector
        const float pointA[2] = {cspanEnd.zSpan.second, cspanEnd.rSpan.second};
        const float pointB[2] = {cspanBar.zSpan.first, cspanBar.rSpan.first};
        lower = zAxisIntersection(pointA,pointB);
        if(cspanEnd.rSpan.first > cspanBar.rSpan.second){
          // if does not overlap, then there is also upper limit
          const float pointC[2] = {cspanEnd.zSpan.first, cspanEnd.rSpan.first};
          const float pointD[2] = {cspanBar.zSpan.second, cspanBar.rSpan.second};
          upper = zAxisIntersection(pointC,pointD);
        }else{
          upper = std::numeric_limits<float>::infinity();
        }
      }else if (cspanEnd.zSpan.first > cspanBar.zSpan.second){
        // if we are on the right hand side of pixel detector
        const float pointA[2] = {cspanEnd.zSpan.first, cspanEnd.rSpan.second};
        const float pointB[2] = {cspanBar.zSpan.second, cspanBar.rSpan.first};
        upper = zAxisIntersection(pointA,pointB);
        if(cspanEnd.rSpan.first > cspanBar.rSpan.second){
          const float pointC[2] = {cspanEnd.zSpan.second,cspanEnd.rSpan.first};
          const float pointD[2] = {cspanBar.zSpan.first, cspanBar.rSpan.second};
          lower = zAxisIntersection(pointC,pointD);
        }else{
          lower = - std::numeric_limits<float>::infinity();
        }
      }else {
        return false;
      }
    }else{
      return false;
    } 
    range =  std::pair<float,float>(lower,upper);
    return true;
  }
}

std::vector<std::pair<edm::VecArray<PixelInactiveAreaFinder::Area, 2>, std::vector<PixelInactiveAreaFinder::LayerSetIndex> > >
PixelInactiveAreaFinder::InactiveAreas::areasAndLayerSets(const GlobalPoint& point, float zwidth) const {
  auto spansLayerSets = spansAndLayerSets(point, zwidth);

  // TODO: try to remove this conversion...
  std::vector<std::pair<VecArray2<Area>, std::vector<LayerSetIndex> > > ret;
  for(auto& item: spansLayerSets) {
    auto& innerSpan = item.first[0];
    auto& outerSpan = item.first[1];
    VecArray2<Area> areas;
    areas.emplace_back(innerSpan.rSpan.first, innerSpan.rSpan.second,
                       innerSpan.phiSpan.first, innerSpan.phiSpan.second,
                       innerSpan.zSpan.first, innerSpan.zSpan.second);
    areas.emplace_back(outerSpan.rSpan.first, outerSpan.rSpan.second,
                       outerSpan.phiSpan.first, outerSpan.phiSpan.second,
                       outerSpan.zSpan.first, outerSpan.zSpan.second);
    ret.emplace_back(std::move(areas), std::move(item.second));
  }

  return ret;
}

std::vector<std::pair<edm::VecArray<PixelInactiveAreaFinder::DetGroupSpan, 2>, std::vector<PixelInactiveAreaFinder::LayerSetIndex> > >
PixelInactiveAreaFinder::InactiveAreas::spansAndLayerSets(const GlobalPoint& point, float zwidth) const {
  // TODO: in the future use 2D-r for the origin for the phi overlap check
  const float zmin = point.z()-zwidth;
  const float zmax = point.z()+zwidth;

  std::vector<std::pair<VecArray2<DetGroupSpan>, std::vector<LayerSetIndex> > > ret;

  LogDebug("PixelInactiveAreaFinder") << "Origin at " << point.x() << "," << point.y() << "," << point.z() << " z half width " << zwidth;

  for(LayerSetIndex i=0, end=inactiveLayerPairIndices_->size(); i<end; ++i) {
    const auto& layerIdxPair = (*inactiveLayerPairIndices_)[i];
    const auto& innerSpans = inactiveSpans_[layerIdxPair.first];
    const auto& outerSpans = inactiveSpans_[layerIdxPair.second];

    for(const auto& innerSpan: innerSpans) {
      for(const auto& outerSpan: outerSpans) {

        if(phiRangesOverlap(innerSpan.phiSpan, outerSpan.phiSpan)) {
          std::pair<float,float> range(0,0);

          bool zOverlap = false;
          const auto innerDet = std::get<0>((*inactiveLayers_)[layerIdxPair.first]);
          const auto outerDet = std::get<0>((*inactiveLayers_)[layerIdxPair.second]);
          if(innerDet == GeomDetEnumerators::PixelBarrel) {
            if(outerDet == GeomDetEnumerators::PixelBarrel)
              zOverlap = getZAxisOverlapRangeBarrel(innerSpan, outerSpan, range);
            else
              zOverlap = getZAxisOverlapRangeBarrelEndcap(innerSpan, outerSpan, range);
          }
          else {
            if(outerDet == GeomDetEnumerators::PixelEndcap)
              zOverlap = getZAxisOverlapRangeEndcap(innerSpan, outerSpan, range);
            else
              throw cms::Exception("LogicError") << "Forward->barrel transition is not supported";
          }

          if(zOverlap && zmin <= range.second && range.first <= zmax) {
#ifdef EDM_ML_DEBUG
            Stream ss;
            for(auto ind: (*layerSetIndexInactiveToActive_)[i]) {
              ss << ind << ",";
            }
            ss << "\n  ";
            detGroupSpanInfo(innerSpan, ss);
            ss << "\n  ";
            detGroupSpanInfo(outerSpan, ss);
            LogTrace("PixelInactiveAreaFinder") << " adding areas for active layer sets " << ss.str();
#endif
            VecArray2<DetGroupSpan> vec;
            vec.emplace_back(innerSpan);
            vec.emplace_back(outerSpan);
            ret.emplace_back(std::move(vec), (*layerSetIndexInactiveToActive_)[i]);
          }
        }
      }
    }
  }

  return ret;
}



PixelInactiveAreaFinder::PixelInactiveAreaFinder(const edm::ParameterSet& iConfig, const std::vector<SeedingLayerSetsBuilder::SeedingLayerId>& seedingLayers,
                                                 const SeedingLayerSetsLooper& seedingLayerSetsLooper,
                                                 edm::ConsumesCollector&& iC):
  debug_(iConfig.getUntrackedParameter<bool>("debug")),
  createPlottingFiles_(iConfig.getUntrackedParameter<bool>("createPlottingFiles")),
  ignoreSingleFPixPanelModules_(iConfig.getParameter<bool>("ignoreSingleFPixPanelModules")),
  inactivePixelDetectorTokens_(edm::vector_transform(iConfig.getParameter<std::vector<edm::InputTag> >("inactivePixelDetectorLabels"), [&](const auto& tag) { return iC.consumes<DetIdCollection>(tag); })),
  badPixelFEDChannelsTokens_(edm::vector_transform(iConfig.getParameter<std::vector<edm::InputTag> >("badPixelFEDChannelCollectionLabels"), [&](const auto& tag) { return iC.consumes<PixelFEDChannelCollection>(tag); }))
{
#ifdef EDM_ML_DEBUG
  for(const auto& layer: seedingLayers) {
    LogTrace("PixelInactiveAreaFinder") << "Input layer subdet " << std::get<0>(layer) << " side " << static_cast<unsigned int>(std::get<1>(layer)) << " layer " << std::get<2>(layer);
  }
#endif

  auto findOrAdd = [&](SeedingLayerId layer) -> unsigned short {
    auto found = std::find(inactiveLayers_.cbegin(), inactiveLayers_.cend(), layer);
    if(found == inactiveLayers_.cend()) {
      auto ret = inactiveLayers_.size();
      inactiveLayers_.push_back(layer);
      return ret;
    }
    return std::distance(inactiveLayers_.cbegin(), found);
  };

  // mapping from active layer pairs to inactive layer pairs
  const auto activeToInactiveMap = createActiveToInactiveMap();

  // convert input layer pairs (that are for active layers) to layer
  // pairs to look for inactive areas
  LayerSetIndex i=0;
  for(const auto& layerSet: seedingLayerSetsLooper.makeRange(seedingLayers)) {
    assert(layerSet.size() == 2);
    auto found = activeToInactiveMap.find(std::make_pair(layerSet[0], layerSet[1]));
    if(found == activeToInactiveMap.end()) {
      throw cms::Exception("Configuration") << "Encountered layer pair " << layerSet[0] << "+" << layerSet[1] << " not found from the internal 'active layer pairs' to 'inactive layer pairs' mapping; either fix the input or the mapping (in PixelInactiveAreaFinder.cc)";
    }

    LogTrace("PixelInactiveAreaFinder") << "Input layer set " << layerSet[0] << "+" << layerSet[1];
    for(const auto& inactiveLayerSet: found->second) {
      auto innerInd = findOrAdd(inactiveLayerSet.first);
      auto outerInd = findOrAdd(inactiveLayerSet.second);

      auto found = std::find(inactiveLayerSetIndices_.cbegin(), inactiveLayerSetIndices_.cend(), std::make_pair(innerInd, outerInd));
      if(found == inactiveLayerSetIndices_.end()) {
        inactiveLayerSetIndices_.emplace_back(innerInd, outerInd);
        layerSetIndexInactiveToActive_.push_back(std::vector<LayerSetIndex>{i});
      }
      else {
        layerSetIndexInactiveToActive_.at(std::distance(inactiveLayerSetIndices_.cbegin(), found)).push_back(i); // TODO: move to operator[] once finished
      }

      LogTrace("PixelInactiveAreaFinder") << " inactive layer set " << inactiveLayerSet.first << "+" << inactiveLayerSet.second;
    }

    ++i;
  }
  

#ifdef EDM_ML_DEBUG
  LogDebug("PixelInactiveAreaFinder") << "All inactive layer sets";
  for(const auto& idxPair: inactiveLayerSetIndices_) {
    LogTrace("PixelInactiveAreaFinder") << " " << inactiveLayers_[idxPair.first] << "+" << inactiveLayers_[idxPair.second];
  }
#endif
}

void PixelInactiveAreaFinder::fillDescriptions(edm::ParameterSetDescription& desc) {
  desc.add<std::vector<edm::InputTag>>("inactivePixelDetectorLabels", std::vector<edm::InputTag>{{edm::InputTag("siPixelDigis")}})->setComment("One or more DetIdCollections of modules to mask on the fly for a given event");
  desc.add<std::vector<edm::InputTag>>("badPixelFEDChannelCollectionLabels", std::vector<edm::InputTag>{{edm::InputTag("siPixelDigis")}})->setComment("One or more PixelFEDChannelCollections of modules+ROCs to mask on the fly for a given event");
  desc.add<bool>("ignoreSingleFPixPanelModules", false);

  desc.addUntracked<bool>("debug", false);
  desc.addUntracked<bool>("createPlottingFiles", false);
}


PixelInactiveAreaFinder::InactiveAreas
PixelInactiveAreaFinder::inactiveAreas(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  // Set data to handles
  {
    edm::ESHandle<TrackerGeometry> trackerGeometry;
    iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
    trackerGeometry_ = trackerGeometry.product();

    edm::ESHandle<TrackerTopology> trackerTopology;
    iSetup.get<TrackerTopologyRcd>().get(trackerTopology);
    trackerTopology_ = trackerTopology.product();
  }

  // assign data to instance variables
  updatePixelDets(iSetup);
  getBadPixelDets(iEvent, iSetup);

  //write files for plotting
  if(createPlottingFiles_) {
    createPlottingFiles();
  }

  // find detGroupSpans ie phi,r,z limits for detector detGroups that are not working
  // returns pair where first is barrel spans and second endcap spans
  DetGroupSpanContainerPair cspans = detGroupSpans();

  // map spans to a vector with consistent indexing with inactiveLayers_ and inactiveLayerSetIndices_
  // TODO: try to move the inner logic towards this direction as well
  std::vector<DetGroupSpanContainer> spans(inactiveLayers_.size());

  auto doWork = [&](const DetGroupSpanContainer& container) {
    for(const auto& span: container) {
      const auto subdet = span.subdetId == PixelSubdetector::PixelBarrel ? GeomDetEnumerators::PixelBarrel : GeomDetEnumerators::PixelEndcap;
      const auto side = (subdet == GeomDetEnumerators::PixelBarrel ? TrackerDetSide::Barrel :
                         (span.zSpan.first < 0 ? TrackerDetSide::NegEndcap : TrackerDetSide::PosEndcap));
      const auto layer = subdet == GeomDetEnumerators::PixelBarrel ? span.layer : span.disk;
      auto found = std::find(inactiveLayers_.begin(), inactiveLayers_.end(), SeedingLayerId(subdet, side, layer));
      if(found != inactiveLayers_.end()) { // it is possible that this layer is ignored by the configuration
        spans[std::distance(inactiveLayers_.begin(), found)].push_back(span);
      }
    }
  };
  doWork(cspans.first);
  doWork(cspans.second);

  auto ret = InactiveAreas(&inactiveLayers_,
                           std::move(spans),
                           &inactiveLayerSetIndices_,
                           &layerSetIndexInactiveToActive_);

  if(debug_) {
    printOverlapSpans(ret);
  }

  return ret;
}

// Functions for fetching date from handles
void PixelInactiveAreaFinder::updatePixelDets(const edm::EventSetup& iSetup) {
  if(!geometryWatcher_.check(iSetup))
    return;

  // deduce the number of ladders per layer and the number of modules per ladder
  {
    // sanity checks
    if(trackerGeometry_->numberOfLayers(PixelSubdetector::PixelBarrel) != 4) {
      throw cms::Exception("NotImplemented") << "This module supports only a detector with 4 pixel barrel layers, the current geometry has " << trackerGeometry_->numberOfLayers(PixelSubdetector::PixelBarrel);
    }
    if(trackerGeometry_->numberOfLayers(PixelSubdetector::PixelEndcap) != 3) {
      throw cms::Exception("NotImplemented") << "This module supports only a detector with 3 pixel forward disks, the current geometry has " << trackerGeometry_->numberOfLayers(PixelSubdetector::PixelEndcap);
    }

    std::array<std::array<unsigned short, 100>, 4> counts = {}; // assume at most 100 ladders per layer
    for(const auto& det: trackerGeometry_->detsPXB()) {
      const auto layer = trackerTopology_->layer(det->geographicalId());
      const auto ladder = trackerTopology_->pxbLadder(det->geographicalId());
      if(ladder >= 100) {
        throw cms::Exception("LogicError") << "Got a ladder with number " << ladder << " while the expected maximum was 100; either something is wrong or the maximum has to be increased.";
      }
      counts[layer-1][ladder-1] += 1; // numbering of layer and ladder starts at 1
    }

    // take number of modules per ladder from the first ladder of the first layer (such better exist)
    // other ladders better have the same number
    nModulesPerLadder = counts[0][0];
    if(nModulesPerLadder == 0) {
      throw cms::Exception("LogicError") << "Ladder 1 of layer 1 has 0 modules, something fishy is going on.";
    }

    LogDebug("PixelInactiveAreaFinder") << "Number of modules per ladder " << nModulesPerLadder << "; below are number of ladders per layer";

    // number of ladders
    for(unsigned layer=0; layer<4; ++layer) {
      nBPixLadders[layer] = std::count_if(counts[layer].begin(), counts[layer].end(), [](unsigned short val) { return val > 0; });
      LogTrace("PixelInactiveAreaFinder") << "BPix layer " << (layer+1) << " has " << nBPixLadders[layer] << " ladders";

      auto fail = std::count_if(counts[layer].begin(), counts[layer].end(), [&](unsigned short val) { return val != nModulesPerLadder && val > 0; });
      if(fail != 0) {
        throw cms::Exception("LogicError") << "Layer " << (layer+1) << " had " << fail << " ladders whose number of modules/ladder differed from the ladder 1 of layer 1 (" << nModulesPerLadder << "). Something fishy is going on.";
      }
    }
  }


  // don't bother with the rest if not needed
  if(!createPlottingFiles_)
    return;

  pixelDetsBarrel_.clear();
  pixelDetsEndcap_.clear();

  for(auto const & geomDetPtr : trackerGeometry_->detsPXB() ) {
    if(geomDetPtr->geographicalId().subdetId() == PixelSubdetector::PixelBarrel){
      pixelDetsBarrel_.push_back(geomDetPtr->geographicalId().rawId());
    }
  }
  for(auto const & geomDetPtr : trackerGeometry_->detsPXF() ) {
    if(geomDetPtr->geographicalId().subdetId() == PixelSubdetector::PixelEndcap){
      pixelDetsEndcap_.push_back(geomDetPtr->geographicalId().rawId());
    }
  }
  std::sort(pixelDetsBarrel_.begin(),pixelDetsBarrel_.end());
  std::sort(pixelDetsEndcap_.begin(),pixelDetsEndcap_.end());
}
void PixelInactiveAreaFinder::getBadPixelDets(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  auto addDetId = [&](const auto id) {
    const auto detid = DetId(id);
    const auto subdet = detid.subdetId();
    if(subdet == PixelSubdetector::PixelBarrel) {
      badPixelDetsBarrel_.push_back(detid.rawId());
    }
    else if(subdet == PixelSubdetector::PixelEndcap) {
      badPixelDetsEndcap_.push_back(detid.rawId());
    }
  };

  // SiPixelQuality
  edm::ESHandle<SiPixelQuality> pixelQuality;
  iSetup.get<SiPixelQualityRcd>().get(pixelQuality);

  for(auto const & disabledModule : pixelQuality->getBadComponentList() ){
    addDetId(disabledModule.DetID);
  }

  // dynamic bad modules
  for(const auto& token: inactivePixelDetectorTokens_) {
    edm::Handle<DetIdCollection> detIds;
    iEvent.getByToken(token, detIds);
    for(const auto& id: *detIds) {
      addDetId(id);
    }
  }

  // dynamic bad ROCs ("Fed25")
  // TODO: consider moving to finer-grained areas for inactive ROCs
  for(const auto& token: badPixelFEDChannelsTokens_) {
    edm::Handle<PixelFEDChannelCollection> pixelFEDChannelCollectionHandle;
    iEvent.getByToken(token, pixelFEDChannelCollectionHandle);
    for(const auto& disabledChannels: *pixelFEDChannelCollectionHandle) {
      addDetId(disabledChannels.detId());
    }
  }

  // remove duplicates
  std::sort(badPixelDetsBarrel_.begin(),badPixelDetsBarrel_.end());
  std::sort(badPixelDetsEndcap_.begin(),badPixelDetsEndcap_.end());
  badPixelDetsBarrel_.erase(std::unique(badPixelDetsBarrel_.begin(),badPixelDetsBarrel_.end()),badPixelDetsBarrel_.end());
  badPixelDetsEndcap_.erase(std::unique(badPixelDetsEndcap_.begin(),badPixelDetsEndcap_.end()),badPixelDetsEndcap_.end());
}
// Printing functions
void PixelInactiveAreaFinder::detInfo(const det_t & det, Stream & ss){
  using std::tie;
  using std::setw;
  using std::showpos;
  using std::noshowpos;
  using std::fixed;
  using std::right;
  using std::left;
  using std::setprecision;
  using std::setfill;
  std::string deli = "; ";
  ss << "id:[" << det << "]" <<deli;
  ss << "subdetid:[" << DetId(det).subdetId() << "]" << deli;
  if(DetId(det).subdetId()==PixelSubdetector::PixelBarrel){
    unsigned int layer  = trackerTopology_->pxbLayer (DetId(det));
    unsigned int ladder = trackerTopology_->pxbLadder(DetId(det));
    unsigned int module = trackerTopology_->pxbModule(DetId(det));
    ss  << "layer:["                      << layer  << "]" << deli 
        << "ladder:[" << right << setw(2) << ladder << "]" << deli 
        << "module:["                     << module << "]" << deli;
  }else if(DetId(det).subdetId()==PixelSubdetector::PixelEndcap){
    unsigned int disk  = trackerTopology_->pxfDisk (DetId(det));
    unsigned int blade = trackerTopology_->pxfBlade(DetId(det));
    unsigned int panel = trackerTopology_->pxfPanel(DetId(det));
    ss  << left << setw(6) << "disk:"  << "["            << right << disk  << "]" << deli 
        << left << setw(7) << "blade:" << "[" << setw(2) << right << blade << "]" << deli 
        << left << setw(7) << "panel:" << "["            << right << panel << "]" << deli;
  }
  float phiA,phiB,zA,zB,rA,rB;
  auto detSurface = trackerGeometry_->idToDet(DetId(det))->surface();
  tie(phiA,phiB) = detSurface.phiSpan();
  tie(zA,zB) = detSurface.zSpan();
  tie(rA,rB) = detSurface.rSpan();
  ss 
    << setprecision(16) 
    << fixed 
    << showpos
    << setfill(' ')
    << "phi:[" << right << setw(12) << phiA << "," << left << setw(12) << phiB << "]" << deli
    << "z:["   << right << setw(7)  << zA   << "," << left << setw(7)  << zB   << "]" << deli << noshowpos
    << "r:["   << right << setw(10) << rA   << "," << left << setw(10) << rB   << "]" << deli;

}
void PixelInactiveAreaFinder::printPixelDets(){
  edm::LogPrint("PixelInactiveAreaFinder") << "Barrel detectors:";
  Stream ss;
  for(auto const & det : pixelDetsBarrel_){
    detInfo(det,ss);
    edm::LogPrint("PixelInactiveAreaFinder") << ss.str();ss.str(std::string());
  }
  edm::LogPrint("PixelInactiveAreaFinder") << "Endcap detectors;";
  for(auto const & det : pixelDetsEndcap_){
    detInfo(det,ss);
    edm::LogPrint("PixelInactiveAreaFinder") << ss.str();ss.str(std::string());
  }
}
void PixelInactiveAreaFinder::printBadPixelDets(){
  edm::LogPrint("PixelInactiveAreaFinder") << "Bad barrel detectors:";
  Stream ss;
  for(auto const & det : badPixelDetsBarrel_){
    detInfo(det,ss);
    edm::LogPrint("PixelInactiveAreaFinder") << ss.str();ss.str(std::string());
  }
  edm::LogPrint("PixelInactiveAreaFinder") << "Endcap detectors;";
  for(auto const & det : badPixelDetsEndcap_){
    detInfo(det,ss);
    edm::LogPrint("PixelInactiveAreaFinder") << ss.str();ss.str(std::string());
  }
}
void PixelInactiveAreaFinder::printBadDetGroups(){
  DetGroupContainer badDetGroupsBar = badDetGroupsBarrel();
  DetGroupContainer badDetGroupsEnd = badDetGroupsEndcap();
  Stream ss;
  for(auto const & detGroup : badDetGroupsBar){
    ss << std::setfill(' ') << std::left << std::setw(16) << "DetGroup:";
    DetGroupSpan cspan;
    getPhiSpanBarrel(detGroup,cspan);
    getZSpan(detGroup,cspan);
    getRSpan(detGroup,cspan);
    detGroupSpanInfo(cspan,ss);
    ss<<std::endl;
    for(auto const & det : detGroup){
      detInfo(det,ss);ss<<std::endl;
    }
    ss<<std::endl;
  }
  for(auto const & detGroup : badDetGroupsEnd){
    ss << std::setfill(' ') << std::left << std::setw(16) << "DetGroup:";
    DetGroupSpan cspan;
    getPhiSpanEndcap(detGroup,cspan);
    getZSpan(detGroup,cspan);
    getRSpan(detGroup,cspan);
    detGroupSpanInfo(cspan,ss);
    ss << std::endl;
    for(auto const & det : detGroup){
      detInfo(det,ss);ss<<std::endl;
    }
    ss << std::endl;
  }
  edm::LogPrint("PixelInactiveAreaFinder")<<ss.str();
}
void PixelInactiveAreaFinder::printBadDetGroupSpans(){
  DetGroupSpanContainerPair cspans = detGroupSpans();
  Stream ss;
  for(auto const & cspan : cspans.first){
    detGroupSpanInfo(cspan,ss);ss<<std::endl;
  }
  for(auto const & cspan : cspans.second){
    detGroupSpanInfo(cspan,ss);ss<<std::endl;
  }
  edm::LogPrint("PixelInactiveAreaFinder") << ss.str();
}
void PixelInactiveAreaFinder::createPlottingFiles(){
  // All detectors to file DETECTORS
  Stream ss;
  std::ofstream fsDet("DETECTORS.txt");
  for(auto const & det : pixelDetsBarrel_){
    detInfo(det,ss);
    ss << std::endl;
  }
  edm::LogPrint("PixelInactiveAreaFinder") << "Endcap detectors;";
  for(auto const & det : pixelDetsEndcap_){
    detInfo(det,ss);
    ss << std::endl;
  }
  fsDet << ss.rdbuf();
  ss.str(std::string());
  // Bad detectors
  std::ofstream fsBadDet("BADDETECTORS.txt");
  for(auto const & det : badPixelDetsBarrel_){
    detInfo(det,ss);
    ss<<std::endl;
  }
  for(auto const & det : badPixelDetsEndcap_){
    detInfo(det,ss);
    ss<<std::endl;
  }
  fsBadDet << ss.rdbuf();
  ss.str(std::string());
  // detgroupspans
  std::ofstream fsSpans("DETGROUPSPANS.txt");
  DetGroupSpanContainerPair cspans = detGroupSpans();
  for(auto const & cspan : cspans.first){
    detGroupSpanInfo(cspan,ss);ss<<std::endl;
  }
  for(auto const & cspan : cspans.second){
    detGroupSpanInfo(cspan,ss);ss<<std::endl;
  }
  fsSpans << ss.rdbuf();
  ss.str(std::string());

}
// Functions for finding bad detGroups
bool PixelInactiveAreaFinder::detWorks(det_t det){
  return 
    std::find(badPixelDetsBarrel_.begin(),badPixelDetsBarrel_.end(),det)
    == badPixelDetsBarrel_.end()
    &&
    std::find(badPixelDetsEndcap_.begin(),badPixelDetsEndcap_.end(),det)
    == badPixelDetsEndcap_.end()
    ;
}
PixelInactiveAreaFinder::DetGroup PixelInactiveAreaFinder::badAdjecentDetsBarrel(const det_t & det){
  using std::remove_if;

  DetGroup adj;
  auto const tTopo = trackerTopology_;
  auto const & detId = DetId(det);
  unsigned int layer  = tTopo->pxbLayer (detId);
  unsigned int ladder = tTopo->pxbLadder(detId);
  unsigned int module = tTopo->pxbModule(detId);
  unsigned int nLads = nBPixLadders[layer];
  //add detectors from next and previous ladder
  adj.push_back( tTopo->pxbDetId( layer, ((ladder-1)+1)%nLads+1, module )() );
  adj.push_back( tTopo->pxbDetId( layer, ((ladder-1)-1+nLads)%nLads+1, module )() );
  //add adjecent detectors from same ladder
  if(module == 1) {
    adj.push_back( tTopo->pxbDetId( layer, ladder, module+1 )() );
  }
  else if(module == nModulesPerLadder) {
    adj.push_back( tTopo->pxbDetId( layer, ladder, module-1 )() );
  }
  else {
    adj.push_back( tTopo->pxbDetId( layer, ladder, module+1 )() );
    adj.push_back( tTopo->pxbDetId( layer, ladder, module-1 )() );
  }
  //remove working detectors from list
  adj.erase(remove_if(adj.begin(),adj.end(), [&](auto c) {return this->detWorks(c);}), adj.end());
  return adj;
}
PixelInactiveAreaFinder::DetGroup PixelInactiveAreaFinder::badAdjecentDetsEndcap(const det_t & det){
  // this might be faster if adjecent 
  using std::tie;
  using std::ignore;
  DetGroup adj;
  Span_t  phiSpan, phiSpanComp;
  float z, zComp;
  unsigned int disk, diskComp;
  auto const & detSurf = trackerGeometry_->idToDet(DetId(det))->surface();
  phiSpan = detSurf.phiSpan();
  tie(z,ignore) = detSurf.zSpan();
  disk = trackerTopology_->pxfDisk(DetId(det));
  // add detectors from same disk whose phi ranges overlap to the adjecent list
  for(auto const & detComp : badPixelDetsEndcap_){
    auto const & detIdComp = DetId(detComp);
    auto const & detSurfComp = trackerGeometry_->idToDet(detIdComp)->surface();
    diskComp = trackerTopology_->pxfDisk(detIdComp);
    phiSpanComp = detSurfComp.phiSpan();
    tie(zComp,ignore) = detSurfComp.zSpan();
    if(det != detComp && disk == diskComp && z*zComp > 0
       && phiRangesOverlap(phiSpan,phiSpanComp)){
      adj.push_back(detComp);
    }
  }
  return adj;
}
PixelInactiveAreaFinder::DetGroup PixelInactiveAreaFinder::reachableDetGroup(const det_t & initDet, DetectorSet & foundDets){
  DetGroup reachableDetGroup;
  std::queue<det_t> workQueue;
  det_t workDet;
  DetGroup badAdjDets;
  foundDets.insert(initDet);
  workQueue.push(initDet);
  reachableDetGroup.push_back(initDet);
  while(!workQueue.empty()){
    workDet = workQueue.front();workQueue.pop();
    if(DetId(workDet).subdetId() == PixelSubdetector::PixelBarrel){
      badAdjDets = this->badAdjecentDetsBarrel(workDet);
    }else if(DetId(workDet).subdetId() == PixelSubdetector::PixelEndcap){
      badAdjDets = this->badAdjecentDetsEndcap(workDet);
    }else {
      badAdjDets = {};
    }
    for(auto const & badDet : badAdjDets){
      if(foundDets.find(badDet) == foundDets.end()){
        reachableDetGroup.push_back(badDet);
        foundDets.insert(badDet);
        workQueue.push(badDet);
      }
    }
  }
  return reachableDetGroup;
}
PixelInactiveAreaFinder::DetGroupContainer PixelInactiveAreaFinder::badDetGroupsBarrel(){
  DetGroupContainer detGroups;
  DetectorSet foundDets;
  for(auto const & badDet : badPixelDetsBarrel_){
    if(foundDets.find(badDet) == foundDets.end()){
      detGroups.push_back(this->reachableDetGroup(badDet,foundDets));
    } 
  }
  return detGroups;
}
PixelInactiveAreaFinder::DetGroupContainer PixelInactiveAreaFinder::badDetGroupsEndcap(){
  DetGroupContainer detGroups;
  DetectorSet foundDets;
  for(auto const & badDet : badPixelDetsEndcap_){
    if(foundDets.find(badDet) == foundDets.end()){
      auto adjacentDets = this->reachableDetGroup(badDet,foundDets);
      if(ignoreSingleFPixPanelModules_ && adjacentDets.size() == 1) {
        // size==1 means that only a single panel of a blade was inactive
        // because of the large overlap with the other panel (i.e.
        // redundancy in the detectory) ignoring these may help to decrease fakes
        continue;
      }
      detGroups.push_back(adjacentDets);
    } 
  }
  return detGroups;
}
// Functions for finding DetGroupSpans
void PixelInactiveAreaFinder::getPhiSpanBarrel(const DetGroup & detGroup, DetGroupSpan & cspan){
  // find phiSpan using ordered vector of unique ladders in detGroup
  if(detGroup.empty()){
    cspan = DetGroupSpan();
    return;
  } else{
    cspan.layer = trackerTopology_->pxbLayer(DetId(detGroup[0]));
    cspan.disk = 0;
  }
  using uint = unsigned int;
  using LadderSet = std::set<uint>;
  using LadVec = std::vector<uint>;
  LadderSet lads;
  for(auto const & det : detGroup){
    lads.insert(trackerTopology_->pxbLadder(DetId(det)));
  }
  LadVec ladv(lads.begin(),lads.end());
  uint nLadders = nBPixLadders[cspan.layer];
  // find start ladder of detGroup
  uint i = 0;
  uint currentLadder = ladv[0];
  uint previousLadder = ladv[ (ladv.size()+i-1) % ladv.size() ];
  // loop until discontinuity is found from vector
  while ( (nLadders+currentLadder-1)%nLadders  == previousLadder  ){
    ++i;
    currentLadder = ladv[i%ladv.size()];
    previousLadder = ladv[ (ladv.size()+i-1)%ladv.size() ];
    if(i == ladv.size()){
      cspan.phiSpan.first  =  std::numeric_limits<float>::epsilon();
      cspan.phiSpan.second = -std::numeric_limits<float>::epsilon();
      return;
    }
  }
  uint startLadder = currentLadder;
  uint endLadder = previousLadder;
  auto detStart = trackerTopology_->pxbDetId(cspan.layer,startLadder,1);
  auto detEnd = trackerTopology_->pxbDetId(cspan.layer,endLadder,1);
  cspan.phiSpan.first  = trackerGeometry_->idToDet(detStart)->surface().phiSpan().first;
  cspan.phiSpan.second = trackerGeometry_->idToDet(detEnd)->surface().phiSpan().second;
}
void PixelInactiveAreaFinder::getPhiSpanEndcap(const DetGroup & detGroup, DetGroupSpan & cspan){
  // this is quite naive/bruteforce method
  // 1) it starts by taking one detector from detGroup and starts to compare it to others
  // 2) when it finds overlapping detector in clockwise direction it starts comparing 
  //    found detector to others
  // 3) search stops until no overlapping detectors in clockwise detector or all detectors
  //    have been work detector
  Stream ss;
  bool found = false;
  auto const tGeom = trackerGeometry_;
  DetGroup::const_iterator startDetIter = detGroup.begin();
  Span_t phiSpan,phiSpanComp;
  unsigned int counter = 0;
  while(!found){
    phiSpan = tGeom->idToDet(DetId(*startDetIter))->surface().phiSpan();
    for(DetGroup::const_iterator compDetIter=detGroup.begin();compDetIter!=detGroup.end();++compDetIter){
      phiSpanComp = tGeom->idToDet(DetId(*compDetIter))->surface().phiSpan();
      if(phiRangesOverlap(phiSpan,phiSpanComp)
         && phiMoreClockwise(phiSpanComp.first,phiSpan.first)
         && startDetIter != compDetIter)
        {
          ++counter;
          if(counter > detGroup.size()){
            cspan.phiSpan.first  =  std::numeric_limits<float>::epsilon();
            cspan.phiSpan.second = -std::numeric_limits<float>::epsilon();
            return;
          }
          startDetIter = compDetIter;break;
        } else if (compDetIter == detGroup.end()-1){
        found = true;
      }
    }
  }
  cspan.phiSpan.first = phiSpan.first;
  // second with same method}
  found = false;
  DetGroup::const_iterator endDetIter = detGroup.begin();
  counter = 0;
  while(!found){
    phiSpan = tGeom->idToDet(DetId(*endDetIter))->surface().phiSpan();
    for(DetGroup::const_iterator compDetIter=detGroup.begin();compDetIter!=detGroup.end();++compDetIter){
      phiSpanComp = tGeom->idToDet(DetId(*compDetIter))->surface().phiSpan();
      if(phiRangesOverlap(phiSpan,phiSpanComp)
         && phiMoreCounterclockwise(phiSpanComp.second,phiSpan.second)
         && endDetIter != compDetIter)
        {
          ++counter;
          if(counter > detGroup.size()){
            cspan.phiSpan.first  =  std::numeric_limits<float>::epsilon();
            cspan.phiSpan.second = -std::numeric_limits<float>::epsilon();
            return;
          }
          endDetIter = compDetIter;break;
        } else if (compDetIter == detGroup.end()-1){
        found = true;
      }
    }
  }
  cspan.phiSpan.second = phiSpan.second;
}
void PixelInactiveAreaFinder::getZSpan(const DetGroup & detGroup, DetGroupSpan & cspan){
  auto cmpFun = [this] (det_t detA, det_t detB){
    return
    trackerGeometry_->idToDet(DetId(detA))->surface().zSpan().first
    <
    trackerGeometry_->idToDet(DetId(detB))->surface().zSpan().first
    ;
  };
    
  auto minmaxIters = std::minmax_element(detGroup.begin(),detGroup.end(),cmpFun);
  cspan.zSpan.first = trackerGeometry_->idToDet(DetId(*(minmaxIters.first)))->surface().zSpan().first;
  cspan.zSpan.second = trackerGeometry_->idToDet(DetId(*(minmaxIters.second)))->surface().zSpan().second;
}
void PixelInactiveAreaFinder::getRSpan(const DetGroup & detGroup, DetGroupSpan & cspan){
  auto cmpFun = [this] (det_t detA, det_t detB){
    return
    trackerGeometry_->idToDet(DetId(detA))->surface().rSpan().first
    <
    trackerGeometry_->idToDet(DetId(detB))->surface().rSpan().first
    ;
  };
    
  auto minmaxIters = std::minmax_element(detGroup.begin(),detGroup.end(),cmpFun);
  cspan.rSpan.first = trackerGeometry_->idToDet(DetId(*(minmaxIters.first)))->surface().rSpan().first;
  cspan.rSpan.second = trackerGeometry_->idToDet(DetId(*(minmaxIters.second)))->surface().rSpan().second;
}
void PixelInactiveAreaFinder::getSpan(const DetGroup & detGroup, DetGroupSpan & cspan){
  auto firstDetIt = detGroup.begin();
  if(firstDetIt != detGroup.end()){
    cspan.subdetId = DetId(*firstDetIt).subdetId();
    if(cspan.subdetId == 1){
      cspan.layer = trackerTopology_->pxbLayer(DetId(*firstDetIt));
      cspan.disk = 0;
      getPhiSpanBarrel(detGroup,cspan);    
    }else if(cspan.subdetId == 2){
      cspan.disk = trackerTopology_->pxfDisk(DetId(*firstDetIt));
      cspan.layer = 0;
      getPhiSpanEndcap(detGroup,cspan);
    }
    getZSpan(detGroup,cspan);
    getRSpan(detGroup,cspan);
  }
}
PixelInactiveAreaFinder::DetGroupSpanContainerPair PixelInactiveAreaFinder::detGroupSpans(){
  DetGroupSpanContainer cspansBarrel;
  DetGroupSpanContainer cspansEndcap;
  DetGroupContainer badDetGroupsBar = badDetGroupsBarrel();
  DetGroupContainer badDetGroupsEnd = badDetGroupsEndcap();
  for(auto const & detGroup : badDetGroupsBar){
    DetGroupSpan cspan;
    getSpan(detGroup,cspan);
    cspansBarrel.push_back(cspan);
  }
  for(auto const & detGroup : badDetGroupsEnd){
    DetGroupSpan cspan;
    getSpan(detGroup,cspan);
    cspansEndcap.push_back(cspan);
  }
  return DetGroupSpanContainerPair(cspansBarrel,cspansEndcap);
}