MuonDetLayerMeasurements.cc

Go to the documentation of this file.

#include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"

#include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"

#include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

typedef MuonTransientTrackingRecHit::MuonRecHitPointer MuonRecHitPointer;
typedef MuonTransientTrackingRecHit::MuonRecHitContainer MuonRecHitContainer;

MuonDetLayerMeasurements::MuonDetLayerMeasurements(edm::InputTag dtlabel,
                                                   edm::InputTag csclabel,
                                                   edm::InputTag rpclabel,
                                                   edm::InputTag gemlabel,
                                                   edm::InputTag me0label,
                                                   edm::ConsumesCollector& iC,
                                                   bool enableDT,
                                                   bool enableCSC,
                                                   bool enableRPC,
                                                   bool enableGEM,
                                                   bool enableME0)
    : enableDTMeasurement(enableDT),
      enableCSCMeasurement(enableCSC),
      enableRPCMeasurement(enableRPC),
      enableGEMMeasurement(enableGEM),
      enableME0Measurement(enableME0),
      theDTRecHits(),
      theCSCRecHits(),
      theRPCRecHits(),
      theGEMRecHits(),
      theDTEventCacheID(0),
      theCSCEventCacheID(0),
      theRPCEventCacheID(0),
      theGEMEventCacheID(0),
      theME0EventCacheID(0),
      theEvent(nullptr) {
  dtToken_ = iC.consumes<DTRecSegment4DCollection>(dtlabel);
  cscToken_ = iC.consumes<CSCSegmentCollection>(csclabel);
  rpcToken_ = iC.consumes<RPCRecHitCollection>(rpclabel);
  gemToken_ = iC.consumes<GEMRecHitCollection>(gemlabel);
  me0Token_ = iC.consumes<ME0SegmentCollection>(me0label);

  static std::atomic<int> procInstance{0};
  std::ostringstream sDT;
  sDT << "MuonDetLayerMeasurements::checkDTRecHits::" << procInstance;
  //  theDTCheckName = sDT.str();
  std::ostringstream sRPC;
  sRPC << "MuonDetLayerMeasurements::checkRPCRecHits::" << procInstance;
  //theRPCCheckName = sRPC.str();
  std::ostringstream sCSC;
  sCSC << "MuonDetLayerMeasurements::checkCSCRecHits::" << procInstance;
  //theCSCCheckName = sCSC.str();
  std::ostringstream sGEM;
  sGEM << "MuonDetLayerMeasurements::checkGEMRecHits::" << procInstance;
  //theGEMCheckName = sGEM.str();
  std::ostringstream sME0;
  sME0 << "MuonDetLayerMeasurements::checkME0RecHits::" << procInstance;
  //theME0CheckName = sME0.str();
  procInstance++;
}

MuonDetLayerMeasurements::~MuonDetLayerMeasurements() {}

MuonRecHitContainer MuonDetLayerMeasurements::recHits(const GeomDet* geomDet, const edm::Event& iEvent) {
  DetId geoId = geomDet->geographicalId();
  theEvent = &iEvent;
  MuonRecHitContainer result;

  if (geoId.subdetId() == MuonSubdetId::DT) {
    if (enableDTMeasurement) {
      checkDTRecHits();

      // Create the ChamberId
      DTChamberId chamberId(geoId.rawId());
      LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(DT): " << chamberId << std::endl;

      // Get the DT-Segment which relies on this chamber
      DTRecSegment4DCollection::range range = theDTRecHits->get(chamberId);

      // Create the MuonTransientTrackingRechit
      for (DTRecSegment4DCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit)
        result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
    }
  }

  else if (geoId.subdetId() == MuonSubdetId::CSC) {
    if (enableCSCMeasurement) {
      checkCSCRecHits();

      // Create the chamber Id
      CSCDetId chamberId(geoId.rawId());
      LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(CSC): " << chamberId << std::endl;

      // Get the CSC-Segment which relies on this chamber
      CSCSegmentCollection::range range = theCSCRecHits->get(chamberId);

      // Create the MuonTransientTrackingRecHit
      for (CSCSegmentCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit)
        result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
    }
  }

  else if (geoId.subdetId() == MuonSubdetId::RPC) {
    if (enableRPCMeasurement) {
      checkRPCRecHits();

      // Create the chamber Id
      RPCDetId chamberId(geoId.rawId());
      LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(RPC): " << chamberId << std::endl;

      // Get the RPC-Segment which relies on this chamber
      RPCRecHitCollection::range range = theRPCRecHits->get(chamberId);

      // Create the MuonTransientTrackingRecHit
      for (RPCRecHitCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit)
        result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
    }
  } else if (geoId.subdetId() == MuonSubdetId::GEM) {
    if (enableGEMMeasurement) {
      checkGEMRecHits();

      // Create the chamber Id
      GEMDetId chamberId(geoId.rawId());

      LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(GEM): " << chamberId << std::endl;

      // Get the GEM-Segment which relies on this chamber
      GEMRecHitCollection::range range = theGEMRecHits->get(chamberId);

      LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements")
          << "Number of GEM rechits available =  " << theGEMRecHits->size() << ", from chamber: " << chamberId
          << std::endl;

      // Create the MuonTransientTrackingRecHit
      for (GEMRecHitCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit)
        result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
      LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "Number of GEM rechits = " << result.size() << std::endl;
    }
  }

  else if (geoId.subdetId() == MuonSubdetId::ME0) {
    LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(ME0): identified" << std::endl;
    if (enableME0Measurement) {
      checkME0RecHits();

      // Create the chamber Id
      ME0DetId chamberId(geoId.rawId());

      // Get the ME0-Segment which relies on this chamber
      // Getting rechits right now, not segments - maybe it should be segments?
      ME0SegmentCollection::range range = theME0RecHits->get(chamberId);

      LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements")
          << "Number of ME0 rechits available =  " << theME0RecHits->size() << ", from chamber: " << chamberId
          << std::endl;

      // Create the MuonTransientTrackingRecHit
      for (ME0SegmentCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit) {
        LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "On ME0 iteration " << std::endl;
        result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
      }
      LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "Number of ME0 rechits = " << result.size() << std::endl;
    }
  } else {
    // wrong type
    throw cms::Exception("MuonDetLayerMeasurements") << "The DetLayer with det " << geoId.det() << " subdet "
                                                     << geoId.subdetId() << " is not a valid Muon DetLayer. ";
    LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "The DetLayer with det " << geoId.det() << " subdet "
                                                       << geoId.subdetId() << " is not a valid Muon DetLayer. ";
  }
  if (enableME0Measurement) {
    LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(ME0): enabled" << std::endl;
  }

  if (enableGEMMeasurement) {
    LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(GEM): enabled" << std::endl;
  }
  return result;
}

void MuonDetLayerMeasurements::checkDTRecHits() {
  checkEvent();
  auto const cacheID = theEvent->cacheIdentifier();
  if (cacheID == theDTEventCacheID)
    return;

  { theEvent->getByToken(dtToken_, theDTRecHits); }
  if (!theDTRecHits.isValid()) {
    throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get DT RecHits";
  }
}

void MuonDetLayerMeasurements::checkCSCRecHits() {
  checkEvent();
  auto cacheID = theEvent->cacheIdentifier();
  if (cacheID == theCSCEventCacheID)
    return;

  {
    theEvent->getByToken(cscToken_, theCSCRecHits);
    theCSCEventCacheID = cacheID;
  }
  if (!theCSCRecHits.isValid()) {
    throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get CSC RecHits";
  }
}

void MuonDetLayerMeasurements::checkRPCRecHits() {
  checkEvent();
  auto cacheID = theEvent->cacheIdentifier();
  if (cacheID == theRPCEventCacheID)
    return;

  {
    theEvent->getByToken(rpcToken_, theRPCRecHits);
    theRPCEventCacheID = cacheID;
  }
  if (!theRPCRecHits.isValid()) {
    throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get RPC RecHits";
  }
}

void MuonDetLayerMeasurements::checkGEMRecHits() {
  checkEvent();
  auto cacheID = theEvent->cacheIdentifier();
  if (cacheID == theGEMEventCacheID)
    return;

  {
    theEvent->getByToken(gemToken_, theGEMRecHits);
    theGEMEventCacheID = cacheID;
  }
  if (!theGEMRecHits.isValid()) {
    throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get GEM RecHits";
  }
}

void MuonDetLayerMeasurements::checkME0RecHits() {
  LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "Checking ME0 RecHits";
  checkEvent();
  auto cacheID = theEvent->cacheIdentifier();
  if (cacheID == theME0EventCacheID)
    return;

  {
    theEvent->getByToken(me0Token_, theME0RecHits);
    theME0EventCacheID = cacheID;
  }
  if (!theME0RecHits.isValid()) {
    throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get ME0 RecHits";
    LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "Cannot get ME0 RecHits";
  }
}

MeasurementContainer MuonDetLayerMeasurements::measurements(const DetLayer* layer,
                                                            const TrajectoryStateOnSurface& startingState,
                                                            const Propagator& prop,
                                                            const MeasurementEstimator& est) {
  checkEvent();
  return measurements(layer, startingState, prop, est, *theEvent);
}

MeasurementContainer MuonDetLayerMeasurements::measurements(const DetLayer* layer,
                                                            const TrajectoryStateOnSurface& startingState,
                                                            const Propagator& prop,
                                                            const MeasurementEstimator& est,
                                                            const edm::Event& iEvent) {
  MeasurementContainer result;

  std::vector<DetWithState> dss = layer->compatibleDets(startingState, prop, est);
  LogDebug("RecoMuon") << "compatibleDets: " << dss.size() << std::endl;

  for (std::vector<DetWithState>::const_iterator detWithStateItr = dss.begin(); detWithStateItr != dss.end();
       ++detWithStateItr) {
    MeasurementContainer detMeasurements =
        measurements(layer, detWithStateItr->first, detWithStateItr->second, est, iEvent);
    result.insert(result.end(), detMeasurements.begin(), detMeasurements.end());
  }

  if (!result.empty())
    sort(result.begin(), result.end(), TrajMeasLessEstim());

  return result;
}

MeasurementContainer MuonDetLayerMeasurements::measurements(const DetLayer* layer,
                                                            const GeomDet* det,
                                                            const TrajectoryStateOnSurface& stateOnDet,
                                                            const MeasurementEstimator& est,
                                                            const edm::Event& iEvent) {
  MeasurementContainer result;

  DetId geoId = det->geographicalId();
  //no chamber here that is actually an me0....
  if (geoId.subdetId() == MuonSubdetId::ME0) {
    if (enableME0Measurement) {
      ME0DetId chamberId(geoId.rawId());
      LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements")
          << "ME0 Chamber ID in measurements: " << chamberId << std::endl;
    }
  }

  // Get the Segments which relies on the GeomDet given by compatibleDets
  MuonRecHitContainer muonRecHits = recHits(det, iEvent);

  // Create the Trajectory Measurement
  for (MuonRecHitContainer::const_iterator rechit = muonRecHits.begin(); rechit != muonRecHits.end(); ++rechit) {
    MeasurementEstimator::HitReturnType estimate = est.estimate(stateOnDet, **rechit);
    LogDebug("RecoMuon") << "Dimension: " << (*rechit)->dimension() << " Chi2: " << estimate.second << std::endl;
    if (estimate.first) {
      result.push_back(TrajectoryMeasurement(stateOnDet, *rechit, estimate.second, layer));
    }
  }

  if (!result.empty())
    sort(result.begin(), result.end(), TrajMeasLessEstim());

  return result;
}

MeasurementContainer MuonDetLayerMeasurements::fastMeasurements(const DetLayer* layer,
                                                                const TrajectoryStateOnSurface& theStateOnDet,
                                                                const TrajectoryStateOnSurface& startingState,
                                                                const Propagator& prop,
                                                                const MeasurementEstimator& est,
                                                                const edm::Event& iEvent) {
  MeasurementContainer result;
  MuonRecHitContainer rhs = recHits(layer, iEvent);
  for (MuonRecHitContainer::const_iterator irh = rhs.begin(); irh != rhs.end(); irh++) {
    MeasurementEstimator::HitReturnType estimate = est.estimate(theStateOnDet, (**irh));
    if (estimate.first) {
      result.push_back(TrajectoryMeasurement(theStateOnDet, (*irh), estimate.second, layer));
    }
  }

  if (!result.empty()) {
    sort(result.begin(), result.end(), TrajMeasLessEstim());
  }

  return result;
}

MeasurementContainer MuonDetLayerMeasurements::fastMeasurements(const DetLayer* layer,
                                                                const TrajectoryStateOnSurface& theStateOnDet,
                                                                const TrajectoryStateOnSurface& startingState,
                                                                const Propagator& prop,
                                                                const MeasurementEstimator& est) {
  checkEvent();
  return fastMeasurements(layer, theStateOnDet, startingState, prop, est, *theEvent);
}

std::vector<TrajectoryMeasurementGroup> MuonDetLayerMeasurements::groupedMeasurements(
    const DetLayer* layer,
    const TrajectoryStateOnSurface& startingState,
    const Propagator& prop,
    const MeasurementEstimator& est) {
  checkEvent();
  return groupedMeasurements(layer, startingState, prop, est, *theEvent);
}

std::vector<TrajectoryMeasurementGroup> MuonDetLayerMeasurements::groupedMeasurements(
    const DetLayer* layer,
    const TrajectoryStateOnSurface& startingState,
    const Propagator& prop,
    const MeasurementEstimator& est,
    const edm::Event& iEvent) {
  std::vector<TrajectoryMeasurementGroup> result;
  // if we want to use the concept of InvalidRecHits,
  // we can reuse LayerMeasurements from TrackingTools/MeasurementDet
  std::vector<DetGroup> groups(layer->groupedCompatibleDets(startingState, prop, est));

  // this should be fixed either in RecoMuon/MeasurementDet/MuonDetLayerMeasurements or
  // RecoMuon/DetLayers/MuRingForwardDoubleLayer
  // and removed the reverse operation in StandAloneMuonFilter::findBestMeasurements

  for (std::vector<DetGroup>::const_iterator grp = groups.begin(); grp != groups.end(); ++grp) {
    std::vector<TrajectoryMeasurement> groupMeasurements;
    for (DetGroup::const_iterator detAndStateItr = grp->begin(); detAndStateItr != grp->end(); ++detAndStateItr) {
      std::vector<TrajectoryMeasurement> detMeasurements =
          measurements(layer, detAndStateItr->det(), detAndStateItr->trajectoryState(), est, iEvent);
      groupMeasurements.insert(groupMeasurements.end(), detMeasurements.begin(), detMeasurements.end());
    }

    if (!groupMeasurements.empty())
      std::sort(groupMeasurements.begin(), groupMeasurements.end(), TrajMeasLessEstim());

    result.push_back(TrajectoryMeasurementGroup(groupMeasurements, *grp));
  }

  return result;
}

void MuonDetLayerMeasurements::setEvent(const edm::Event& event) { theEvent = &event; }

void MuonDetLayerMeasurements::checkEvent() const {
  if (!theEvent)
    throw cms::Exception("MuonDetLayerMeasurements") << "The event has not been set";
}

MuonRecHitContainer MuonDetLayerMeasurements::recHits(const DetLayer* layer, const edm::Event& iEvent) {
  MuonRecHitContainer rhs;

  std::vector<const GeomDet*> gds = layer->basicComponents();

  for (std::vector<const GeomDet*>::const_iterator igd = gds.begin(); igd != gds.end(); igd++) {
    MuonRecHitContainer detHits = recHits(*igd, iEvent);
    rhs.insert(rhs.end(), detHits.begin(), detHits.end());
  }
  return rhs;
}

MuonRecHitContainer MuonDetLayerMeasurements::recHits(const DetLayer* layer) {
  checkEvent();
  return recHits(layer, *theEvent);
}