HitExtractorSTRP.cc

Go to the documentation of this file.

#include "HitExtractorSTRP.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"

#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/Event.h"

#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"

#include "DataFormats/Common/interface/ContainerMask.h"

#include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h"
#include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"

#include <tuple>

#include <iostream>

using namespace ctfseeding;
using namespace std;
using namespace edm;

HitExtractorSTRP::HitExtractorSTRP(GeomDetEnumerators::SubDetector subdet,
                                   TrackerDetSide side,
                                   int idLayer,
                                   float iminGoodCharge)
    : theLayerSubDet(subdet),
      theSide(side),
      theIdLayer(idLayer),
      minAbsZ(0),
      theMinRing(1),
      theMaxRing(0),
      hasMatchedHits(false),
      hasRPhiHits(false),
      hasStereoHits(false),
      hasRingSelector(false),
      hasSimpleRphiHitsCleaner(true) {
  minGoodCharge = iminGoodCharge;
  if (minGoodCharge > 0)
    skipClusters = true;
}

void HitExtractorSTRP::useSkipClusters_(const edm::InputTag& m, edm::ConsumesCollector& iC) {
  theSkipClusters = iC.consumes<SkipClustersCollection>(m);
}

void HitExtractorSTRP::useRingSelector(int minRing, int maxRing) {
  hasRingSelector = true;
  theMinRing = minRing;
  theMaxRing = maxRing;
}

bool HitExtractorSTRP::ringRange(int ring) const {
  if (!hasRingSelector)
    return true;
  return (ring >= theMinRing) & (ring <= theMaxRing);
}

bool HitExtractorSTRP::skipThis(
    DetId id,
    OmniClusterRef const& clus,
    edm::Handle<edm::ContainerMask<edmNew::DetSetVector<SiStripCluster> > >& stripClusterMask) const {
  if (maskCluster && (stripClusterMask->mask(clus.key())))
    return true;

  if
    UNLIKELY(minGoodCharge <= 0) return false;
  return siStripClusterTools::chargePerCM(id, *clus.cluster_strip()) <= minGoodCharge;
}

std::pair<bool, ProjectedSiStripRecHit2D*> HitExtractorSTRP::skipThis(
    const TkTransientTrackingRecHitBuilder& ttrhBuilder,
    TkHitRef matched,
    edm::Handle<edm::ContainerMask<edmNew::DetSetVector<SiStripCluster> > >& stripClusterMask) const {
  const SiStripMatchedRecHit2D& hit = (SiStripMatchedRecHit2D const&)(matched);

  assert(dynamic_cast<SiStripMatchedRecHit2D const*>(&matched));

  auto id = hit.geographicalId();
  ProjectedSiStripRecHit2D* replaceMe = nullptr;
  bool rejectSt = skipThis(id, hit.stereoClusterRef(), stripClusterMask);
  bool rejectMono = skipThis(id, hit.monoClusterRef(), stripClusterMask);

  if ((!rejectSt) & (!rejectMono)) {
    // keepit
    return std::make_pair(false, replaceMe);
  }

  if (failProjection || (rejectSt & rejectMono)) {
    //only skip if both hits are done
    return std::make_pair(true, replaceMe);
  }

  // replace with one

  auto cloner = ttrhBuilder.cloner();
  replaceMe = cloner.project(hit, rejectSt, TrajectoryStateOnSurface()).release();
  if (rejectSt)
    LogDebug("HitExtractorSTRP") << "a matched hit is partially masked, and the mono hit got projected onto: "
                                 << replaceMe->geographicalId().rawId() << " key: " << hit.monoClusterRef().key();
  else
    LogDebug("HitExtractorSTRP") << "a matched hit is partially masked, and the stereo hit got projected onto: "
                                 << replaceMe->geographicalId().rawId() << " key: " << hit.stereoClusterRef().key();

  return std::make_pair(true, replaceMe);
}

void HitExtractorSTRP::cleanedOfClusters(const TkTransientTrackingRecHitBuilder& ttrhBuilder,
                                         const edm::Event& ev,
                                         HitExtractor::Hits& hits,
                                         bool matched,
                                         unsigned int cleanFrom) const {
  LogDebug("HitExtractorPIX") << "getting: " << hits.size() << " in input.";
  edm::Handle<SkipClustersCollection> stripClusterMask;
  if (maskCluster)
    ev.getByToken(theSkipClusters, stripClusterMask);
  unsigned int skipped = 0;
  unsigned int projected = 0;
  for (unsigned int iH = cleanFrom; iH < hits.size(); ++iH) {
    assert(hits[iH]->isValid());
    auto id = hits[iH]->geographicalId();
    if (matched) {
      bool replace;
      ProjectedSiStripRecHit2D* replaceMe;
      std::tie(replace, replaceMe) = skipThis(ttrhBuilder, *hits[iH], stripClusterMask);
      if (replace) {
        if (!replaceMe) {
          LogDebug("HitExtractorSTRP") << "skipping a matched hit on :" << hits[iH]->geographicalId().rawId();
          skipped++;
        } else
          projected++;
        hits[iH].reset(replaceMe);
        if (replaceMe == nullptr)
          assert(hits[iH].empty());
        else
          assert(hits[iH].isOwn());
      }
    } else if (skipThis(id, hits[iH]->firstClusterRef(), stripClusterMask)) {
      LogDebug("HitExtractorSTRP") << "skipping a hit on :" << hits[iH]->geographicalId().rawId() << " key: ";
      skipped++;
      hits[iH].reset();
    }
  }
  //  remove empty elements...
  auto last = std::remove_if(hits.begin() + cleanFrom, hits.end(), [](HitPointer const& p) { return p.empty(); });
  hits.resize(last - hits.begin());

  // std::cout << "HitExtractorSTRP " <<"skipped :"<<skipped<<" strip rechits because of clusters and projected: "<<projected << std::endl;
  LogDebug("HitExtractorSTRP") << "skipped :" << skipped
                               << " strip rechits because of clusters and projected: " << projected;
}

HitExtractor::Hits HitExtractorSTRP::hits(const TkTransientTrackingRecHitBuilder& ttrhBuilder,
                                          const edm::Event& ev,
                                          const edm::EventSetup& es) const {
  HitExtractor::Hits result;
  unsigned int cleanFrom = 0;

  //Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopoHandle;
  es.get<TrackerTopologyRcd>().get(tTopoHandle);
  const TrackerTopology* const tTopo = tTopoHandle.product();

  //
  // TIB
  //
  if (theLayerSubDet == GeomDetEnumerators::TIB) {
    if (hasMatchedHits) {
      edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
      ev.getByToken(theMatchedHits, matchedHits);
      if (skipClusters)
        cleanFrom = result.size();
      range2SeedingHits(*matchedHits, result, tTopo->tibDetIdLayerComparator(theIdLayer));
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, true, cleanFrom);
    }
    if (hasRPhiHits) {
      edm::Handle<SiStripRecHit2DCollection> rphiHits;
      ev.getByToken(theRPhiHits, rphiHits);
      if (hasMatchedHits) {
        if (!hasSimpleRphiHitsCleaner) {  // this is a brutal "cleaning". Add something smarter in the future
          if (skipClusters)
            cleanFrom = result.size();
          range2SeedingHits(*rphiHits, result, tTopo->tibDetIdLayerComparator(theIdLayer));
          if (skipClusters)
            cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
        }
      } else {
        if (skipClusters)
          cleanFrom = result.size();
        range2SeedingHits(*rphiHits, result, tTopo->tibDetIdLayerComparator(theIdLayer));
        if (skipClusters)
          cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
      }
    }
    if (hasStereoHits) {
      edm::Handle<SiStripRecHit2DCollection> stereoHits;
      ev.getByToken(theStereoHits, stereoHits);
      if (skipClusters)
        cleanFrom = result.size();
      range2SeedingHits(*stereoHits, result, tTopo->tibDetIdLayerComparator(theIdLayer));
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
    }
  }

  //
  // TID
  //
  else if (theLayerSubDet == GeomDetEnumerators::TID) {
    if (hasMatchedHits) {
      edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
      ev.getByToken(theMatchedHits, matchedHits);
      if (skipClusters)
        cleanFrom = result.size();
      auto getter = tTopo->tidDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
      SiStripMatchedRecHit2DCollection::Range range = matchedHits->equal_range(getter.first, getter.second);
      for (SiStripMatchedRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
        int ring = tTopo->tidRing(it->detId());
        if (!ringRange(ring))
          continue;
        for (SiStripMatchedRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end;
             ++hit) {
          result.emplace_back(*hit);
        }
      }
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, true, cleanFrom);
    }
    if (hasRPhiHits) {
      edm::Handle<SiStripRecHit2DCollection> rphiHits;
      ev.getByToken(theRPhiHits, rphiHits);
      if (skipClusters)
        cleanFrom = result.size();
      auto getter = tTopo->tidDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
      SiStripRecHit2DCollection::Range range = rphiHits->equal_range(getter.first, getter.second);
      for (SiStripRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
        int ring = tTopo->tidRing(it->detId());
        if (!ringRange(ring))
          continue;
        if ((SiStripDetId(it->detId()).partnerDetId() != 0) && hasSimpleRphiHitsCleaner)
          continue;  // this is a brutal "cleaning". Add something smarter in the future
        for (SiStripRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
          result.emplace_back(*hit);
        }
      }
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
    }
    if (hasStereoHits) {
      edm::Handle<SiStripRecHit2DCollection> stereoHits;
      ev.getByToken(theStereoHits, stereoHits);
      if (skipClusters)
        cleanFrom = result.size();
      auto getter = tTopo->tidDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
      SiStripRecHit2DCollection::Range range = stereoHits->equal_range(getter.first, getter.second);
      for (SiStripRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
        int ring = tTopo->tidRing(it->detId());
        if (!ringRange(ring))
          continue;
        for (SiStripRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
          result.emplace_back(*hit);
        }
      }
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
    }
  }
  //
  // TOB
  //
  else if (theLayerSubDet == GeomDetEnumerators::TOB) {
    if (hasMatchedHits) {
      edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
      ev.getByToken(theMatchedHits, matchedHits);
      if (skipClusters)
        cleanFrom = result.size();
      if (minAbsZ > 0.) {
        auto getter = tTopo->tobDetIdLayerComparator(theIdLayer);
        SiStripMatchedRecHit2DCollection::Range range = matchedHits->equal_range(getter.first, getter.second);
        for (SiStripMatchedRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
          for (SiStripMatchedRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end;
               ++hit) {
            if (fabs(hit->globalPosition().z()) >= minAbsZ)
              result.emplace_back(*hit);
          }
        }
      } else {
        range2SeedingHits(*matchedHits, result, tTopo->tobDetIdLayerComparator(theIdLayer));
      }
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, true, cleanFrom);
    }
    if (hasRPhiHits) {
      edm::Handle<SiStripRecHit2DCollection> rphiHits;
      ev.getByToken(theRPhiHits, rphiHits);
      if (hasMatchedHits) {
        if (!hasSimpleRphiHitsCleaner) {  // this is a brutal "cleaning". Add something smarter in the future
          if (skipClusters)
            cleanFrom = result.size();
          range2SeedingHits(*rphiHits, result, tTopo->tobDetIdLayerComparator(theIdLayer));
          if (skipClusters)
            cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
        }
      } else {
        if (skipClusters)
          cleanFrom = result.size();
        range2SeedingHits(*rphiHits, result, tTopo->tobDetIdLayerComparator(theIdLayer));
        if (skipClusters)
          cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
      }
    }
    if (hasStereoHits) {
      edm::Handle<SiStripRecHit2DCollection> stereoHits;
      ev.getByToken(theStereoHits, stereoHits);
      if (skipClusters)
        cleanFrom = result.size();
      range2SeedingHits(*stereoHits, result, tTopo->tobDetIdLayerComparator(theIdLayer));
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
    }
  }

  //
  // TEC
  //
  else if (theLayerSubDet == GeomDetEnumerators::TEC) {
    if (hasMatchedHits) {
      edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
      ev.getByToken(theMatchedHits, matchedHits);
      if (skipClusters)
        cleanFrom = result.size();
      auto getter = tTopo->tecDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
      SiStripMatchedRecHit2DCollection::Range range = matchedHits->equal_range(getter.first, getter.second);
      for (SiStripMatchedRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
        int ring = tTopo->tecRing(it->detId());
        if (!ringRange(ring))
          continue;
        for (SiStripMatchedRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end;
             ++hit) {
          result.emplace_back(*hit);
        }
      }
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, true, cleanFrom);
    }
    if (hasRPhiHits) {
      edm::Handle<SiStripRecHit2DCollection> rphiHits;
      ev.getByToken(theRPhiHits, rphiHits);
      if (skipClusters)
        cleanFrom = result.size();
      auto getter = tTopo->tecDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
      SiStripRecHit2DCollection::Range range = rphiHits->equal_range(getter.first, getter.second);
      for (SiStripRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
        int ring = tTopo->tecRing(it->detId());
        if (!ringRange(ring))
          continue;
        if ((SiStripDetId(it->detId()).partnerDetId() != 0) && hasSimpleRphiHitsCleaner)
          continue;  // this is a brutal "cleaning". Add something smarter in the future
        for (SiStripRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
          result.emplace_back(*hit);
        }
      }
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
    }
    if (hasStereoHits) {
      edm::Handle<SiStripRecHit2DCollection> stereoHits;
      ev.getByToken(theStereoHits, stereoHits);
      if (skipClusters)
        cleanFrom = result.size();
      auto getter = tTopo->tecDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
      SiStripRecHit2DCollection::Range range = stereoHits->equal_range(getter.first, getter.second);
      for (SiStripRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
        int ring = tTopo->tecRing(it->detId());
        if (!ringRange(ring))
          continue;
        for (SiStripRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
          result.emplace_back(*hit);
        }
      }
      if (skipClusters)
        cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
    }
  }

  LogDebug("HitExtractorSTRP") << " giving: " << result.size() << " out";
  // std::cout << "HitExtractorSTRP "<<" giving: "<<result.size() << " for charge cut " << minGoodCharge << std::endl;
  return result;
}