MuonShowerInformationFiller.cc

Go to the documentation of this file.

#include "RecoMuon/MuonIdentification/interface/MuonShowerInformationFiller.h"

// system include files
#include <memory>
#include <algorithm>
#include <iostream>
#include <numeric>

// user include files
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/PluginManager/interface/PluginManager.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
#include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
#include "DataFormats/GeometrySurface/interface/BoundDisk.h"
#include "DataFormats/GeometrySurface/interface/Bounds.h"
#include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"
#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
#include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"
#include "RecoMuon/TrackingTools/interface/MuonTrajectoryBuilder.h"
#include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
#include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"

#include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"

#include "DataFormats/MuonReco/interface/MuonShower.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"

using namespace std;
using namespace edm;

namespace {
  //
  // Sort function optimized for the case where the desired comparison
  // is extravagantly expensive.  To use, replace
  //
  // auto LessPhi = [](const auto& lhs, const auto& rhs) {
  //   return (lhs->globalPosition().barePhi() < rhs->globalPosition().barePhi())
  // }
  // stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessPhi);
  //
  // with
  //
  // auto calcPhi = [](const auto& hit) { return hit->globalPosition().barePhi(); };
  // muonRecHits[stat] = sort_all_indexed(muonRecHits[stat], std::less(), calcPhi);
  //
  // This calculates the values to be compared once in O(n) time, instead of
  // O(n*log(n)) times in the comparison function.
  //
  template <typename RandomAccessSequence, typename Predicate, typename Transform>
  RandomAccessSequence sort_all_indexed(const RandomAccessSequence& s, Predicate p, Transform t) {
    std::vector<size_t> idx(s.size());
    std::iota(idx.begin(), idx.end(), 0);

    // fill the cache of the values to be sorted
    using valueCacheType = std::invoke_result_t<decltype(t), typename RandomAccessSequence::value_type>;
    std::vector<valueCacheType> valcache(s.size());
    std::transform(s.begin(), s.end(), valcache.begin(), t);

    // sort the indices of the value cache
    auto idxComp = [&valcache, p](auto i1, auto i2) { return p(valcache[i1], valcache[i2]); };
    std::stable_sort(idx.begin(), idx.end(), idxComp);

    // fill the sorted output vector
    RandomAccessSequence r(s.size());
    for (size_t i = 0; i < s.size(); ++i) {
      r[i] = s[idx[i]];
    }
    return r;
  }
}  // namespace

//
// Constructor
//
MuonShowerInformationFiller::MuonShowerInformationFiller(const edm::ParameterSet& par, edm::ConsumesCollector& iC)
    : theService(nullptr),
      theDTRecHitLabel(par.getParameter<InputTag>("DTRecSegmentLabel")),
      theCSCRecHitLabel(par.getParameter<InputTag>("CSCRecSegmentLabel")),
      theCSCSegmentsLabel(par.getParameter<InputTag>("CSCSegmentLabel")),
      theDT4DRecSegmentLabel(par.getParameter<InputTag>("DT4DRecSegmentLabel")) {
  theDTRecHitToken = iC.consumes<DTRecHitCollection>(theDTRecHitLabel);
  theCSCRecHitToken = iC.consumes<CSCRecHit2DCollection>(theCSCRecHitLabel);
  theCSCSegmentsToken = iC.consumes<CSCSegmentCollection>(theCSCSegmentsLabel);
  theDT4DRecSegmentToken = iC.consumes<DTRecSegment4DCollection>(theDT4DRecSegmentLabel);

  theTrackerToken = iC.esConsumes();
  theTrackingGeometryToken = iC.esConsumes();
  theFieldToken = iC.esConsumes();
  theCSCGeometryToken = iC.esConsumes();
  theDTGeometryToken = iC.esConsumes();

  edm::ParameterSet serviceParameters = par.getParameter<edm::ParameterSet>("ServiceParameters");
  theService = new MuonServiceProxy(serviceParameters, edm::ConsumesCollector(iC));

  auto trackerRecHitBuilderName = par.getParameter<string>("TrackerRecHitBuilder");
  theTrackerRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", trackerRecHitBuilderName));
  auto muonRecHitBuilderName = par.getParameter<string>("MuonRecHitBuilder");
  theMuonRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", muonRecHitBuilderName));

  theCacheId_TRH = 0;
  theCacheId_MT = 0;

  category_ = "MuonShowerInformationFiller";

  for (int istat = 0; istat < 4; istat++) {
    theStationShowerDeltaR.push_back(0.);
    theStationShowerTSize.push_back(0.);
    theAllStationHits.push_back(0);
    theCorrelatedStationHits.push_back(0);
  }
}

//
// Destructor
//
MuonShowerInformationFiller::~MuonShowerInformationFiller() {
  if (theService)
    delete theService;
}

//
//Fill the MuonShower struct
//
reco::MuonShower MuonShowerInformationFiller::fillShowerInformation(const reco::Muon& muon,
                                                                    const edm::Event& iEvent,
                                                                    const edm::EventSetup& iSetup) {
  reco::MuonShower returnShower;

  // Update the services
  theService->update(iSetup);
  setEvent(iEvent);
  setServices(theService->eventSetup());

  fillHitsByStation(muon);
  std::vector<int> nStationHits = theAllStationHits;
  std::vector<int> nStationCorrelatedHits = theCorrelatedStationHits;
  std::vector<float> stationShowerSizeT = theStationShowerTSize;
  std::vector<float> stationShowerDeltaR = theStationShowerDeltaR;

  returnShower.nStationHits = nStationHits;
  returnShower.nStationCorrelatedHits = nStationCorrelatedHits;
  returnShower.stationShowerSizeT = stationShowerSizeT;
  returnShower.stationShowerDeltaR = stationShowerDeltaR;

  return returnShower;
}

//
// Set Event
//
void MuonShowerInformationFiller::setEvent(const edm::Event& event) {
  // get all the necesary products
  event.getByToken(theDTRecHitToken, theDTRecHits);
  event.getByToken(theCSCRecHitToken, theCSCRecHits);
  event.getByToken(theCSCSegmentsToken, theCSCSegments);
  event.getByToken(theDT4DRecSegmentToken, theDT4DRecSegments);

  for (int istat = 0; istat < 4; istat++) {
    theStationShowerDeltaR.at(istat) = 0.;
    theStationShowerTSize.at(istat) = 0.;
    theAllStationHits.at(istat) = 0;
    theCorrelatedStationHits.at(istat) = 0;
  }
}

//
// Set services
//
void MuonShowerInformationFiller::setServices(const EventSetup& setup) {
  // DetLayer Geometry
  theTrackingGeometry = setup.getHandle(theTrackingGeometryToken);
  theField = setup.getHandle(theFieldToken);
  theTracker = setup.getHandle(theTrackerToken);
  theCSCGeometry = setup.getHandle(theCSCGeometryToken);
  theDTGeometry = setup.getHandle(theDTGeometryToken);

  // Transient Rechit Builders
  unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
  if (newCacheId_TRH != theCacheId_TRH) {
    theTrackerRecHitBuilder = setup.getHandle(theTrackerRecHitBuilderToken);
    theMuonRecHitBuilder = setup.getHandle(theMuonRecHitBuilderToken);
  }
}

//
// Get hits owned by segments
//
TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::hitsFromSegments(
    const GeomDet* geomDet,
    edm::Handle<DTRecSegment4DCollection> dtSegments,
    edm::Handle<CSCSegmentCollection> cscSegments) const {
  MuonTransientTrackingRecHit::MuonRecHitContainer segments;

  DetId geoId = geomDet->geographicalId();

  if (geoId.subdetId() == MuonSubdetId::DT) {
    DTChamberId chamberId(geoId.rawId());

    // loop on segments 4D
    DTRecSegment4DCollection::id_iterator chamberIdIt;
    for (chamberIdIt = dtSegments->id_begin(); chamberIdIt != dtSegments->id_end(); ++chamberIdIt) {
      if (*chamberIdIt != chamberId)
        continue;

      // Get the range for the corresponding ChamberId
      DTRecSegment4DCollection::range range = dtSegments->get((*chamberIdIt));

      for (DTRecSegment4DCollection::const_iterator iseg = range.first; iseg != range.second; ++iseg) {
        if (iseg->dimension() != 4)
          continue;
        segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*iseg));
      }
    }
  } else if (geoId.subdetId() == MuonSubdetId::CSC) {
    CSCDetId did(geoId.rawId());

    for (CSCSegmentCollection::id_iterator chamberId = cscSegments->id_begin(); chamberId != cscSegments->id_end();
         ++chamberId) {
      if ((*chamberId).chamber() != did.chamber())
        continue;

      // Get the range for the corresponding ChamberId
      CSCSegmentCollection::range range = cscSegments->get((*chamberId));

      for (CSCSegmentCollection::const_iterator iseg = range.first; iseg != range.second; ++iseg) {
        if (iseg->dimension() != 3)
          continue;
        segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*iseg));
      }
    }
  } else {
    LogTrace(category_) << "Segments are not built in RPCs" << endl;
  }

  TransientTrackingRecHit::ConstRecHitContainer allhitscorrelated;

  if (segments.empty())
    return allhitscorrelated;

  TransientTrackingRecHit::ConstRecHitPointer muonRecHit(segments.front());
  allhitscorrelated = MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit, 2);

  if (segments.size() == 1)
    return allhitscorrelated;

  for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseg = segments.begin() + 1;
       iseg != segments.end();
       ++iseg) {
    TransientTrackingRecHit::ConstRecHitPointer muonRecHit((*iseg));
    TransientTrackingRecHit::ConstRecHitContainer hits1 =
        MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit, 2);

    for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit1 = hits1.begin(); ihit1 != hits1.end();
         ++ihit1) {
      bool usedbefore = false;
      //unused      DetId thisID = (*ihit1)->geographicalId();
      //LocalPoint lp1dinsegHit = (*ihit1)->localPosition();
      GlobalPoint gp1dinsegHit = (*ihit1)->globalPosition();

      for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit2 = allhitscorrelated.begin();
           ihit2 != allhitscorrelated.end();
           ++ihit2) {
        //unused        DetId thisID2 = (*ihit2)->geographicalId();
        //LocalPoint lp1dinsegHit2 = (*ihit2)->localPosition();
        GlobalPoint gp1dinsegHit2 = (*ihit2)->globalPosition();

        if ((gp1dinsegHit2 - gp1dinsegHit).mag() < 1.0)
          usedbefore = true;
      }
      if (!usedbefore)
        allhitscorrelated.push_back(*ihit1);
    }
  }

  return allhitscorrelated;
}

//
// Find cluster
//

TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::findThetaCluster(
    const TransientTrackingRecHit::ConstRecHitContainer& muonRecHitsIn, const GlobalPoint& refpoint) const {
  if (muonRecHitsIn.empty())
    return muonRecHitsIn;

  //clustering step by theta
  float step = 0.05;
  TransientTrackingRecHit::ConstRecHitContainer result;

  auto muonRecHitsTmp = sort_all_indexed(muonRecHitsIn, std::less(), AbsDThetaTransform(refpoint));

  for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = muonRecHitsTmp.begin();
       ihit != muonRecHitsTmp.end() - 1;
       ++ihit) {
    if (fabs((*(ihit + 1))->globalPosition().theta() - (*ihit)->globalPosition().theta()) < step) {
      result.push_back(*ihit);
    } else {
      break;
    }
  }

  return result;
}

//
//Used to treat overlap region
//
MuonTransientTrackingRecHit::MuonRecHitContainer MuonShowerInformationFiller::findPerpCluster(
    const MuonTransientTrackingRecHit::MuonRecHitContainer& muonRecHitsIn) const {
  if (muonRecHitsIn.empty())
    return muonRecHitsIn;

  auto muonRecHitsTmp = sort_all_indexed(muonRecHitsIn, std::less(), PerpTransform());

  MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator seedhit =
      min_element(muonRecHitsTmp.begin(), muonRecHitsTmp.end(), LessPerp());

  MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihigh = seedhit;
  MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ilow = seedhit;

  float step = 0.1;
  while (ihigh != muonRecHitsTmp.end() - 1 &&
         (fabs((*(ihigh + 1))->globalPosition().perp() - (*ihigh)->globalPosition().perp()) < step)) {
    ihigh++;
  }
  while (ilow != muonRecHitsTmp.begin() &&
         (fabs((*ilow)->globalPosition().perp() - (*(ilow - 1))->globalPosition().perp()) < step)) {
    ilow--;
  }

  MuonTransientTrackingRecHit::MuonRecHitContainer result(ilow, ihigh);

  return result;
}

//
// Get compatible dets
//
vector<const GeomDet*> MuonShowerInformationFiller::getCompatibleDets(const reco::Track& track) const {
  vector<const GeomDet*> total;
  total.reserve(1000);

  LogTrace(category_) << "Consider a track " << track.p() << " eta: " << track.eta() << " phi " << track.phi() << endl;

  TrajectoryStateOnSurface innerTsos = trajectoryStateTransform::innerStateOnSurface(
      track, *theService->trackingGeometry(), &*theService->magneticField());
  TrajectoryStateOnSurface outerTsos = trajectoryStateTransform::outerStateOnSurface(
      track, *theService->trackingGeometry(), &*theService->magneticField());

  GlobalPoint innerPos = innerTsos.globalPosition();
  GlobalPoint outerPos = outerTsos.globalPosition();

  vector<GlobalPoint> allCrossingPoints;

  const vector<const DetLayer*>& dtlayers = theService->detLayerGeometry()->allDTLayers();

  for (auto iLayer = dtlayers.begin(); iLayer != dtlayers.end(); ++iLayer) {
    // crossing points of track with cylinder
    GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const BarrelDetLayer*>(*iLayer));

    // check if point is inside the detector
    if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500) &&
        (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0)))
      allCrossingPoints.push_back(xPoint);
  }

  allCrossingPoints = sort_all_indexed(allCrossingPoints, std::less(), MagTransform(innerPos));

  vector<const GeomDet*> tempDT;

  for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
    tempDT = dtPositionToDets(*ipos);
    vector<const GeomDet*>::const_iterator begin = tempDT.begin();
    vector<const GeomDet*>::const_iterator end = tempDT.end();

    for (; begin != end; ++begin) {
      total.push_back(*begin);
    }
  }
  allCrossingPoints.clear();

  const vector<const DetLayer*>& csclayers = theService->detLayerGeometry()->allCSCLayers();
  for (auto iLayer = csclayers.begin(); iLayer != csclayers.end(); ++iLayer) {
    GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const ForwardDetLayer*>(*iLayer));

    // check if point is inside the detector
    if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500.0) &&
        (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0)))
      allCrossingPoints.push_back(xPoint);
  }
  allCrossingPoints = sort_all_indexed(allCrossingPoints, std::less(), MagTransform(innerPos));

  vector<const GeomDet*> tempCSC;
  for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
    tempCSC = cscPositionToDets(*ipos);
    vector<const GeomDet*>::const_iterator begin = tempCSC.begin();
    vector<const GeomDet*>::const_iterator end = tempCSC.end();

    for (; begin != end; ++begin) {
      total.push_back(*begin);
    }
  }

  return total;
}

//
// Intersection point of track with barrel layer
//
GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
                                                       const GlobalPoint& p2,
                                                       const BarrelDetLayer* dl) const {
  const BoundCylinder& bc = dl->specificSurface();
  return crossingPoint(p1, p2, bc);
}

GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
                                                       const GlobalPoint& p2,
                                                       const Cylinder& cyl) const {
  float radius = cyl.radius();

  GlobalVector dp = p1 - p2;
  float slope = dp.x() / dp.y();
  float a = p1.x() - slope * p1.y();

  float n2 = (1 + slope * slope);
  float n1 = 2 * a * slope;
  float n0 = a * a - radius * radius;

  float y1 = 9999;
  float y2 = 9999;
  if (n1 * n1 - 4 * n2 * n0 > 0) {
    y1 = (-n1 + sqrt(n1 * n1 - 4 * n2 * n0)) / (2 * n2);
    y2 = (-n1 - sqrt(n1 * n1 - 4 * n2 * n0)) / (2 * n2);
  }

  float x1 = p1.x() + slope * (y1 - p1.y());
  float x2 = p1.x() + slope * (y2 - p1.y());

  float slopeZ = dp.z() / dp.y();

  float z1 = p1.z() + slopeZ * (y1 - p1.y());
  float z2 = p1.z() + slopeZ * (y2 - p1.y());

  // there are two crossing points, return the one that is in the same quadrant as point of extrapolation
  if ((p2.x() * x1 > 0) && (y1 * p2.y() > 0) && (z1 * p2.z() > 0)) {
    return GlobalPoint(x1, y1, z1);
  } else {
    return GlobalPoint(x2, y2, z2);
  }
}

//
// Intersection point of track with a forward layer
//
GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
                                                       const GlobalPoint& p2,
                                                       const ForwardDetLayer* dl) const {
  const BoundDisk& bc = dl->specificSurface();
  return crossingPoint(p1, p2, bc);
}

GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
                                                       const GlobalPoint& p2,
                                                       const BoundDisk& disk) const {
  float diskZ = disk.position().z();
  int endcap = diskZ > 0 ? 1 : (diskZ < 0 ? -1 : 0);
  diskZ = diskZ + endcap * dynamic_cast<const SimpleDiskBounds&>(disk.bounds()).thickness() / 2.;

  GlobalVector dp = p1 - p2;

  float slopeZ = dp.z() / dp.y();
  float y1 = diskZ / slopeZ;

  float slopeX = dp.z() / dp.x();
  float x1 = diskZ / slopeX;

  float z1 = diskZ;

  if (p2.z() * z1 > 0) {
    return GlobalPoint(x1, y1, z1);
  } else {
    return GlobalPoint(0, 0, 0);
  }
}

//
// GeomDets along the track in DT
//
vector<const GeomDet*> MuonShowerInformationFiller::dtPositionToDets(const GlobalPoint& gp) const {
  int minwheel = -3;
  int maxwheel = -3;
  if (gp.z() < -680.0) {
    minwheel = -3;
    maxwheel = -3;
  } else if (gp.z() < -396.0) {
    minwheel = -2;
    maxwheel = -1;
  } else if (gp.z() < -126.8) {
    minwheel = -2;
    maxwheel = 0;
  } else if (gp.z() < 126.8) {
    minwheel = -1;
    maxwheel = 1;
  } else if (gp.z() < 396.0) {
    minwheel = 0;
    maxwheel = 2;
  } else if (gp.z() < 680.0) {
    minwheel = 1;
    maxwheel = 2;
  } else {
    minwheel = 3;
    maxwheel = 3;
  }

  int station = 5;
  if (gp.perp() > 680.0 && gp.perp() < 755.0)
    station = 4;
  else if (gp.perp() > 580.0)
    station = 3;
  else if (gp.perp() > 480.0)
    station = 2;
  else if (gp.perp() > 380.0)
    station = 1;
  else
    station = 0;

  vector<int> sectors;

  float phistep = M_PI / 6;

  float phigp = (float)gp.barePhi();

  if (fabs(deltaPhi(phigp, 0 * phistep)) < phistep)
    sectors.push_back(1);
  if (fabs(deltaPhi(phigp, phistep)) < phistep)
    sectors.push_back(2);
  if (fabs(deltaPhi(phigp, 2 * phistep)) < phistep)
    sectors.push_back(3);
  if (fabs(deltaPhi(phigp, 3 * phistep)) < phistep) {
    sectors.push_back(4);
    if (station == 4)
      sectors.push_back(13);
  }
  if (fabs(deltaPhi(phigp, 4 * phistep)) < phistep)
    sectors.push_back(5);
  if (fabs(deltaPhi(phigp, 5 * phistep)) < phistep)
    sectors.push_back(6);
  if (fabs(deltaPhi(phigp, 6 * phistep)) < phistep)
    sectors.push_back(7);
  if (fabs(deltaPhi(phigp, 7 * phistep)) < phistep)
    sectors.push_back(8);
  if (fabs(deltaPhi(phigp, 8 * phistep)) < phistep)
    sectors.push_back(9);
  if (fabs(deltaPhi(phigp, 9 * phistep)) < phistep) {
    sectors.push_back(10);
    if (station == 4)
      sectors.push_back(14);
  }
  if (fabs(deltaPhi(phigp, 10 * phistep)) < phistep)
    sectors.push_back(11);
  if (fabs(deltaPhi(phigp, 11 * phistep)) < phistep)
    sectors.push_back(12);

  LogTrace(category_) << "DT position to dets" << endl;
  LogTrace(category_) << "number of sectors to consider: " << sectors.size() << endl;
  LogTrace(category_) << "station: " << station << " wheels: " << minwheel << " " << maxwheel << endl;

  vector<const GeomDet*> result;
  if (station > 4 || station < 1)
    return result;
  if (minwheel > 2 || maxwheel < -2)
    return result;

  for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
    for (int iwheel = minwheel; iwheel != maxwheel + 1; ++iwheel) {
      DTChamberId chamberid(iwheel, station, (*isector));
      result.push_back(theService->trackingGeometry()->idToDet(chamberid));
    }
  }

  LogTrace(category_) << "number of GeomDets for this track: " << result.size() << endl;

  return result;
}

//
// GeomDets along the track in CSC
//
vector<const GeomDet*> MuonShowerInformationFiller::cscPositionToDets(const GlobalPoint& gp) const {
  // determine the endcap side
  int endcap = 0;
  if (gp.z() > 0) {
    endcap = 1;
  } else {
    endcap = 2;
  }

  // determine the csc station and range of rings
  int station = 5;

  // check all rings in a station
  if (fabs(gp.z()) > 1000. && fabs(gp.z()) < 1055.0) {
    station = 4;
  } else if (fabs(gp.z()) > 910.0 && fabs(gp.z()) < 965.0) {
    station = 3;
  } else if (fabs(gp.z()) > 800.0 && fabs(gp.z()) < 860.0) {
    station = 2;
  } else if (fabs(gp.z()) > 570.0 && fabs(gp.z()) < 730.0) {
    station = 1;
  }

  vector<int> sectors;

  float phistep1 = M_PI / 18.;  //for all the rings except first rings for stations > 1
  float phistep2 = M_PI / 9.;
  float phigp = (float)gp.barePhi();

  int ring = -1;

  // determine the ring
  if (station == 1) {
    //FIX ME!!!
    if (gp.perp() > 100 && gp.perp() < 270)
      ring = 1;
    else if (gp.perp() > 270 && gp.perp() < 450)
      ring = 2;
    else if (gp.perp() > 450 && gp.perp() < 695)
      ring = 3;
    else if (gp.perp() > 100 && gp.perp() < 270)
      ring = 4;

  } else if (station == 2) {
    if (gp.perp() > 140 && gp.perp() < 350)
      ring = 1;
    else if (gp.perp() > 350 && gp.perp() < 700)
      ring = 2;

  } else if (station == 3) {
    if (gp.perp() > 160 && gp.perp() < 350)
      ring = 1;
    else if (gp.perp() > 350 && gp.perp() < 700)
      ring = 2;

  } else if (station == 4) {
    if (gp.perp() > 175 && gp.perp() < 350)
      ring = 1;
    else if (gp.perp() > 350 && gp.perp() < 700)
      ring = 2;
  }

  if (station > 1 && ring == 1) {
    // we have 18 sectors in that case
    for (int i = 0; i < 18; i++) {
      if (fabs(deltaPhi(phigp, i * phistep2)) < phistep2)
        sectors.push_back(i + 1);
    }

  } else {
    // we have 36 sectors in that case
    for (int i = 0; i < 36; i++) {
      if (fabs(deltaPhi(phigp, i * phistep1)) < phistep1)
        sectors.push_back(i + 1);
    }
  }

  LogTrace(category_) << "CSC position to dets" << endl;
  LogTrace(category_) << "ring: " << ring << endl;
  LogTrace(category_) << "endcap: " << endcap << endl;
  LogTrace(category_) << "station: " << station << endl;
  LogTrace(category_) << "CSC number of sectors to consider: " << sectors.size() << endl;

  // check exceptional cases
  vector<const GeomDet*> result;
  if (station > 4 || station < 1)
    return result;
  if (endcap == 0)
    return result;
  if (ring == -1)
    return result;

  int minlayer = 1;
  int maxlayer = 6;

  for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
    for (int ilayer = minlayer; ilayer != maxlayer + 1; ++ilayer) {
      CSCDetId cscid(endcap, station, ring, (*isector), ilayer);
      result.push_back(theService->trackingGeometry()->idToDet(cscid));
    }
  }

  return result;
}

//
//Set class members
//
void MuonShowerInformationFiller::fillHitsByStation(const reco::Muon& muon) {
  reco::TrackRef track;
  if (muon.isGlobalMuon())
    track = muon.globalTrack();
  else if (muon.isStandAloneMuon())
    track = muon.outerTrack();
  else
    return;

  // split 1D rechits by station
  vector<MuonRecHitContainer> muonRecHits(4);

  // split rechits from segs by station
  vector<TransientTrackingRecHit::ConstRecHitContainer> muonCorrelatedHits(4);

  // get vector of GeomDets compatible with a track
  vector<const GeomDet*> compatibleLayers = getCompatibleDets(*track);

  // for special cases: CSC station 1
  MuonRecHitContainer tmpCSC1;
  bool dtOverlapToCheck = false;
  bool cscOverlapToCheck = false;

  for (vector<const GeomDet*>::const_iterator igd = compatibleLayers.begin(); igd != compatibleLayers.end(); igd++) {
    // get det id
    DetId geoId = (*igd)->geographicalId();

    // skip tracker hits
    if (geoId.det() != DetId::Muon)
      continue;

    // DT
    if (geoId.subdetId() == MuonSubdetId::DT) {
      // get station
      DTChamberId detid(geoId.rawId());
      int station = detid.station();
      int wheel = detid.wheel();

      // get rechits from segments per station
      TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp =
          hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);
      TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
      TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();

      for (; hits_begin != hits_end; ++hits_begin) {
        muonCorrelatedHits.at(station - 1).push_back(*hits_begin);
      }

      //check overlap certain wheels and stations
      if (abs(wheel) == 2 && station != 4 && station != 1)
        dtOverlapToCheck = true;

      // loop over all superlayers of a DT chamber
      for (int isuperlayer = DTChamberId::minSuperLayerId; isuperlayer != DTChamberId::maxSuperLayerId + 1;
           ++isuperlayer) {
        // loop over all layers inside the superlayer
        for (int ilayer = DTChamberId::minLayerId; ilayer != DTChamberId::maxLayerId + 1; ++ilayer) {
          DTLayerId lid(detid, isuperlayer, ilayer);
          DTRecHitCollection::range dRecHits = theDTRecHits->get(lid);
          for (DTRecHitCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second; ++rechit) {
            vector<const TrackingRecHit*> subrechits = (*rechit).recHits();
            for (vector<const TrackingRecHit*>::iterator irechit = subrechits.begin(); irechit != subrechits.end();
                 ++irechit) {
              muonRecHits.at(station - 1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd), &**irechit));
            }
          }
        }
      }
    } else if (geoId.subdetId() == MuonSubdetId::CSC) {
      // get station
      CSCDetId did(geoId.rawId());
      int station = did.station();
      int ring = did.ring();

      //get rechits from segments by station
      TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp =
          hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);
      TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
      TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();

      for (; hits_begin != hits_end; ++hits_begin) {
        muonCorrelatedHits.at(station - 1).push_back(*hits_begin);
      }

      if ((station == 1 && ring == 3) && dtOverlapToCheck)
        cscOverlapToCheck = true;

      // split 1D rechits by station
      CSCRecHit2DCollection::range dRecHits = theCSCRecHits->get(did);
      for (CSCRecHit2DCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second; ++rechit) {
        if (!cscOverlapToCheck) {
          muonRecHits.at(station - 1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd), &*rechit));
        } else {
          tmpCSC1.push_back(MuonTransientTrackingRecHit::specificBuild((&**igd), &*rechit));

          //sort by perp, then insert to appropriate container
          MuonRecHitContainer temp = findPerpCluster(tmpCSC1);
          if (temp.empty())
            continue;

          float center;
          if (temp.size() > 1) {
            center = (temp.front()->globalPosition().perp() + temp.back()->globalPosition().perp()) / 2.;
          } else {
            center = temp.front()->globalPosition().perp();
          }
          temp.clear();

          if (center > 550.) {
            muonRecHits.at(2).insert(muonRecHits.at(2).end(), tmpCSC1.begin(), tmpCSC1.end());
          } else {
            muonRecHits.at(1).insert(muonRecHits.at(1).end(), tmpCSC1.begin(), tmpCSC1.end());
          }
          tmpCSC1.clear();
        }
      }
    } else if (geoId.subdetId() == MuonSubdetId::RPC) {
      LogTrace(category_) << "Wrong subdet id" << endl;
    }
  }  //loop over GeomDets compatible with a track

  // calculate number of all and correlated hits
  for (int stat = 0; stat < 4; stat++) {
    theCorrelatedStationHits[stat] = muonCorrelatedHits.at(stat).size();
    theAllStationHits[stat] = muonRecHits[stat].size();
  }
  LogTrace(category_) << "Hits used by the segments, by station " << theCorrelatedStationHits.at(0) << " "
                      << theCorrelatedStationHits.at(1) << " " << theCorrelatedStationHits.at(2) << " "
                      << theCorrelatedStationHits.at(3) << endl;

  LogTrace(category_) << "All DT 1D/CSC 2D  hits, by station " << theAllStationHits.at(0) << " "
                      << theAllStationHits.at(1) << " " << theAllStationHits.at(2) << " " << theAllStationHits.at(3)
                      << endl;

  //station shower sizes
  MuonTransientTrackingRecHit::MuonRecHitContainer muonRecHitsPhiBest;
  TransientTrackingRecHit::ConstRecHitContainer muonRecHitsThetaTemp, muonRecHitsThetaBest;
  // send station hits to the clustering algorithm
  for (int stat = 0; stat != 4; stat++) {
    const size_t nhit = muonRecHits[stat].size();
    if (nhit < 2)
      continue;  // Require at least 2 hits
    muonRecHitsPhiBest.clear();
    muonRecHitsPhiBest.reserve(nhit);

    // Cluster seeds by global position phi. Best cluster is chosen to give greatest dphi
    // Sort by phi (complexity = NLogN with enough memory, or = NLog^2N for insufficient mem)
    muonRecHits[stat] = sort_all_indexed(muonRecHits[stat], std::less(), PhiTransform());

    // Search for gaps (complexity = N)
    std::vector<size_t> clUppers;
    for (size_t ihit = 0; ihit < nhit; ++ihit) {
      const size_t jhit = (ihit + 1) % nhit;
      const double phi1 = muonRecHits[stat].at(ihit)->globalPosition().barePhi();
      const double phi2 = muonRecHits[stat].at(jhit)->globalPosition().barePhi();

      const double dphi = std::abs(deltaPhi(phi1, phi2));
      if (dphi >= 0.05)
        clUppers.push_back(ihit);
    }

    //station shower sizes
    double dphimax = 0;
    if (clUppers.empty()) {
      // No gaps - there is only one cluster. Take all of them
      const double refPhi = muonRecHits[stat].at(0)->globalPosition().barePhi();
      double dphilo = 0, dphihi = 0;
      for (auto& hit : muonRecHits[stat]) {
        muonRecHitsPhiBest.push_back(hit);
        const double phi = hit->globalPosition().barePhi();
        dphilo = std::min(dphilo, deltaPhi(refPhi, phi));
        dphihi = std::max(dphihi, deltaPhi(refPhi, phi));
      }
      dphimax = std::abs(dphihi + dphilo);
    } else {
      // Loop over gaps to find the one with maximum dphi(begin, end)
      // By construction, number of gap must be greater than 1.
      size_t bestUpper = 0, bestLower = 0;
      for (auto icl = clUppers.begin(); icl != clUppers.end(); ++icl) {
        // upper bound of this cluster
        const size_t upper = *icl;
        // lower bound is +1 of preceeding upper bound
        const auto prevCl = (icl == clUppers.begin()) ? clUppers.end() : icl;
        const size_t lower = (*(prevCl - 1) + 1) % nhit;

        // At least two hit for a cluster
        if (upper == lower)
          continue;

        const double phi1 = muonRecHits[stat].at(upper)->globalPosition().barePhi();
        const double phi2 = muonRecHits[stat].at(lower)->globalPosition().barePhi();

        const double dphi = std::abs(deltaPhi(phi1, phi2));
        if (dphimax < dphi) {
          dphimax = dphi;
          bestUpper = upper;
          bestLower = lower;
        }
      }

      if (bestUpper > bestLower) {
        muonRecHitsPhiBest.reserve(bestUpper - bestLower + 1);
        std::copy(muonRecHits[stat].begin() + bestLower,
                  muonRecHits[stat].begin() + bestUpper + 1,
                  std::back_inserter(muonRecHitsPhiBest));
      } else if (bestUpper < bestLower) {
        muonRecHitsPhiBest.reserve(bestUpper + (nhit - bestLower) + 1);
        std::copy(muonRecHits[stat].begin(),
                  muonRecHits[stat].begin() + bestUpper + 1,
                  std::back_inserter(muonRecHitsPhiBest));
        std::copy(
            muonRecHits[stat].begin() + bestLower, muonRecHits[stat].end(), std::back_inserter(muonRecHitsPhiBest));
      }
    }

    //fill showerTs
    if (!muonRecHitsPhiBest.empty()) {
      muonRecHits[stat] = sort_all_indexed(muonRecHitsPhiBest, std::less(), AbsMagTransform());
      GlobalPoint refpoint = muonRecHits[stat].front()->globalPosition();
      theStationShowerTSize.at(stat) = refpoint.mag() * dphimax;
    }

    //for theta
    if (!muonCorrelatedHits.at(stat).empty()) {
      float dthetamax = 0;
      for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iseed = muonCorrelatedHits.at(stat).begin();
           iseed != muonCorrelatedHits.at(stat).end();
           ++iseed) {
        if (!(*iseed)->isValid())
          continue;
        GlobalPoint refpoint = (*iseed)->globalPosition();  //starting from the one with smallest value of phi
        muonRecHitsThetaTemp.clear();
        muonRecHitsThetaTemp = findThetaCluster(muonCorrelatedHits.at(stat), refpoint);
        if (muonRecHitsThetaTemp.size() > 1) {
          float dtheta = fabs((float)muonRecHitsThetaTemp.back()->globalPosition().theta() -
                              (float)muonRecHitsThetaTemp.front()->globalPosition().theta());
          if (dtheta > dthetamax) {
            dthetamax = dtheta;
            muonRecHitsThetaBest = muonRecHitsThetaTemp;
          }
        }  //at least two hits
      }    //loop over seeds
    }      //not empty container2

    //fill deltaRs
    if (muonRecHitsThetaBest.size() > 1 && muonRecHitsPhiBest.size() > 1)
      theStationShowerDeltaR.at(stat) = sqrt(pow(deltaPhi(muonRecHitsPhiBest.front()->globalPosition().barePhi(),
                                                          muonRecHitsPhiBest.back()->globalPosition().barePhi()),
                                                 2) +
                                             pow(muonRecHitsThetaBest.front()->globalPosition().theta() -
                                                     muonRecHitsThetaBest.back()->globalPosition().theta(),
                                                 2));

  }  //loop over station

  LogTrace(category_) << "deltaR around a track containing all the station hits, by station "
                      << theStationShowerDeltaR.at(0) << " " << theStationShowerDeltaR.at(1) << " "
                      << theStationShowerDeltaR.at(2) << " " << theStationShowerDeltaR.at(3) << endl;

  LogTrace(category_) << "Transverse cluster size, by station " << theStationShowerTSize.at(0) << " "
                      << theStationShowerTSize.at(1) << " " << theStationShowerTSize.at(2) << " "
                      << theStationShowerTSize.at(3) << endl;

  return;
}