SubCollectionProducers/src/TrackClusterRemover.cc

0001 #include "FWCore/Framework/interface/global/EDProducer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0005 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0006
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008
0009 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0010 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0011 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0012
0013 #include "DataFormats/Common/interface/ValueMap.h"
0014 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0015 #include "DataFormats/Provenance/interface/ProductID.h"
0016 #include "DataFormats/Common/interface/ContainerMask.h"
0017
0018 #include "DataFormats/DetId/interface/DetId.h"
0019
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021 #include "DataFormats/TrackerRecHit2D/interface/ClusterRemovalInfo.h"
0022
0023 #include "TrackingTools/PatternTools/interface/TrackCollectionTokens.h"
0024
0025 #include "RecoTracker/TransientTrackingRecHit/interface/Traj2TrackHits.h"
0026
0027 #include <limits>
0028
0029 namespace {
0030
0031   class TrackClusterRemover final : public edm::global::EDProducer<> {
0032   public:
0033     TrackClusterRemover(const edm::ParameterSet& iConfig);
0034     ~TrackClusterRemover() override {}
0035     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0036
0037   private:
0038     void produce(edm::StreamID, edm::Event& evt, const edm::EventSetup&) const override;
0039
0040     using PixelMaskContainer = edm::ContainerMask<edmNew::DetSetVector<SiPixelCluster>>;
0041     using StripMaskContainer = edm::ContainerMask<edmNew::DetSetVector<SiStripCluster>>;
0042
0043     using QualityMaskCollection = std::vector<unsigned char>;
0044
0045     const unsigned char maxChi2x5_;
0046     const int minNumberOfLayersWithMeasBeforeFiltering_;
0047     const reco::TrackBase::TrackQuality trackQuality_;
0048
0049     const TrackCollectionTokens trajectories_;
0050     edm::EDGetTokenT<QualityMaskCollection> srcQuals;
0051
0052     edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster>> pixelClusters_;
0053     edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster>> stripClusters_;
0054
0055     edm::EDGetTokenT<PixelMaskContainer> oldPxlMaskToken_;
0056     edm::EDGetTokenT<StripMaskContainer> oldStrMaskToken_;
0057
0058     // backward compatibility during transition period
0059     edm::EDGetTokenT<edm::ValueMap<int>> overrideTrkQuals_;
0060   };
0061
0062   void TrackClusterRemover::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0063     edm::ParameterSetDescription desc;
0064     desc.add<edm::InputTag>("trajectories", edm::InputTag());
0065     desc.add<edm::InputTag>("trackClassifier", edm::InputTag("", "QualityMasks"));
0066     desc.add<edm::InputTag>("pixelClusters", edm::InputTag("siPixelClusters"));
0067     desc.add<edm::InputTag>("stripClusters", edm::InputTag("siStripClusters"));
0068     desc.add<edm::InputTag>("oldClusterRemovalInfo", edm::InputTag());
0069
0070     desc.add<std::string>("TrackQuality", "highPurity");
0071     desc.add<double>("maxChi2", 30.);
0072     desc.add<int>("minNumberOfLayersWithMeasBeforeFiltering", 0);
0073     // old mode
0074     desc.add<edm::InputTag>("overrideTrkQuals", edm::InputTag());
0075
0076     descriptions.add("trackClusterRemover", desc);
0077   }
0078
0079   TrackClusterRemover::TrackClusterRemover(const edm::ParameterSet& iConfig)
0080       : maxChi2x5_(Traj2TrackHits::toChi2x5(iConfig.getParameter<double>("maxChi2"))),
0081         minNumberOfLayersWithMeasBeforeFiltering_(
0082             iConfig.getParameter<int>("minNumberOfLayersWithMeasBeforeFiltering")),
0083         trackQuality_(reco::TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"))),
0084
0085         trajectories_(iConfig.getParameter<edm::InputTag>("trajectories"), consumesCollector())
0086
0087   {
0088     auto const& pixelClusters = iConfig.getParameter<edm::InputTag>("pixelClusters");
0089     auto const& stripClusters = iConfig.getParameter<edm::InputTag>("stripClusters");
0090
0091     if (pixelClusters.label().empty() && stripClusters.label().empty())
0092       throw edm::Exception(edm::errors::Configuration)
0093           << "Configuration Error: TrackClusterRemover used without input cluster collections";
0094     if (!pixelClusters.label().empty()) {
0095       pixelClusters_ = consumes<edmNew::DetSetVector<SiPixelCluster>>(pixelClusters);
0096       produces<edm::ContainerMask<edmNew::DetSetVector<SiPixelCluster>>>();
0097     }
0098     if (!stripClusters.label().empty()) {
0099       stripClusters_ = consumes<edmNew::DetSetVector<SiStripCluster>>(stripClusters);
0100       produces<edm::ContainerMask<edmNew::DetSetVector<SiStripCluster>>>();
0101     }
0102     // old mode
0103     auto const& overrideTrkQuals = iConfig.getParameter<edm::InputTag>("overrideTrkQuals");
0104     if (!overrideTrkQuals.label().empty())
0105       overrideTrkQuals_ = consumes<edm::ValueMap<int>>(overrideTrkQuals);
0106
0107     auto const& classifier = iConfig.getParameter<edm::InputTag>("trackClassifier");
0108     if (!classifier.label().empty())
0109       srcQuals = consumes<QualityMaskCollection>(classifier);
0110
0111     auto const& oldClusterRemovalInfo = iConfig.getParameter<edm::InputTag>("oldClusterRemovalInfo");
0112     if (!oldClusterRemovalInfo.label().empty()) {
0113       oldPxlMaskToken_ = consumes<PixelMaskContainer>(oldClusterRemovalInfo);
0114       oldStrMaskToken_ = consumes<StripMaskContainer>(oldClusterRemovalInfo);
0115     }
0116   }
0117
0118   void TrackClusterRemover::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0119     edm::Handle<edmNew::DetSetVector<SiPixelCluster>> pixelClusters;
0120     if (!pixelClusters_.isUninitialized())
0121       iEvent.getByToken(pixelClusters_, pixelClusters);
0122     edm::Handle<edmNew::DetSetVector<SiStripCluster>> stripClusters;
0123     if (!stripClusters_.isUninitialized())
0124       iEvent.getByToken(stripClusters_, stripClusters);
0125
0126     std::vector<bool> collectedStrips;
0127     std::vector<bool> collectedPixels;
0128
0129     if (!oldPxlMaskToken_.isUninitialized()) {
0130       edm::Handle<PixelMaskContainer> oldPxlMask;
0131       edm::Handle<StripMaskContainer> oldStrMask;
0132       iEvent.getByToken(oldPxlMaskToken_, oldPxlMask);
0133       iEvent.getByToken(oldStrMaskToken_, oldStrMask);
0134       LogDebug("TrackClusterRemover") << "to merge in, " << oldStrMask->size() << " strp, " << oldPxlMask->size()
0135                                       << " pxl";
0136
0137       // std::cout <<"TrackClusterRemover "<<"to merge in, "<<oldStrMask->size()<<" strp and "<<oldPxlMask->size()<<" pxl" << std::endl;
0138       oldStrMask->copyMaskTo(collectedStrips);
0139       oldPxlMask->copyMaskTo(collectedPixels);
0140       if (!stripClusters_.isUninitialized()) {
0141         assert(stripClusters->dataSize() >= collectedStrips.size());
0142         collectedStrips.resize(stripClusters->dataSize(), false);
0143         //std::cout << "TrackClusterRemover " <<"total strip already to skip: "
0144         //  <<std::count(collectedStrips.begin(),collectedStrips.end(),true) <<std::endl;
0145       }
0146     } else {
0147       if (!stripClusters_.isUninitialized())
0148         collectedStrips.resize(stripClusters->dataSize(), false);
0149       if (!pixelClusters_.isUninitialized())
0150         collectedPixels.resize(pixelClusters->dataSize(), false);
0151     }
0152
0153     // loop over trajectories, filter, mask clusters../
0154
0155     unsigned char qualMask = ~0;
0156     if (trackQuality_ != reco::TrackBase::undefQuality)
0157       qualMask = 1 << trackQuality_;
0158
0159     auto const& tracks = trajectories_.tracks(iEvent);
0160     auto s = tracks.size();
0161
0162     // assert(s==trajs.size());
0163
0164     QualityMaskCollection oldStyle;
0165     QualityMaskCollection const* pquals = nullptr;
0166
0167     if (!overrideTrkQuals_.isUninitialized()) {
0168       edm::Handle<edm::ValueMap<int>> quals;
0169       iEvent.getByToken(overrideTrkQuals_, quals);
0170       assert(s == (*quals).size());
0171
0172       oldStyle.resize(s, 0);
0173       for (auto i = 0U; i < s; ++i)
0174         if ((*quals).get(i) > 0)
0175           oldStyle[i] = (255) & (*quals).get(i);
0176       pquals = &oldStyle;
0177     }
0178
0179     if (!srcQuals.isUninitialized()) {
0180       edm::Handle<QualityMaskCollection> hqual;
0181       iEvent.getByToken(srcQuals, hqual);
0182       pquals = hqual.product();
0183     }
0184     // if (!pquals) std::cout << "no qual collection" << std::endl;
0185     for (auto i = 0U; i < s; ++i) {
0186       const reco::Track& track = tracks[i];
0187       bool goodTk = (pquals) ? (*pquals)[i] & qualMask : track.quality(trackQuality_);
0188       if (!goodTk)
0189         continue;
0190       if (track.hitPattern().trackerLayersWithMeasurement() < minNumberOfLayersWithMeasBeforeFiltering_)
0191         continue;
0192
0193       auto const& chi2sX5 = track.extra()->chi2sX5();
0194       assert(chi2sX5.size() == track.recHitsSize());
0195       auto hb = track.recHitsBegin();
0196       for (unsigned int h = 0; h < track.recHitsSize(); h++) {
0197         auto recHit = *(hb + h);
0198         auto const& hit = *recHit;
0199         if (!trackerHitRTTI::isFromDet(hit))
0200           continue;
0201         if (chi2sX5[h] > maxChi2x5_)
0202           continue;  // skip outliers
0203         auto const& thit = reinterpret_cast<BaseTrackerRecHit const&>(hit);
0204         auto const& cluster = thit.firstClusterRef();
0205         if (!stripClusters_.isUninitialized() && cluster.isStrip())
0206           collectedStrips[cluster.key()] = true;
0207         if (!pixelClusters_.isUninitialized() && cluster.isPixel())
0208           collectedPixels[cluster.key()] = true;
0209         if (trackerHitRTTI::isMatched(thit))
0210           collectedStrips[reinterpret_cast<SiStripMatchedRecHit2D const&>(hit).stereoClusterRef().key()] = true;
0211       }
0212     }
0213
0214     // std::cout << " => collectedStrips: " << collectedStrips.size() << std::endl;
0215     if (!stripClusters_.isUninitialized()) {
0216       auto removedStripClusterMask = std::make_unique<StripMaskContainer>(
0217           edm::RefProd<edmNew::DetSetVector<SiStripCluster>>(stripClusters), collectedStrips);
0218       LogDebug("TrackClusterRemover") << "total strip to skip: "
0219                                       << std::count(collectedStrips.begin(), collectedStrips.end(), true);
0220       // std::cout << "TrackClusterRemover " <<"total strip to skip: "<<std::count(collectedStrips.begin(),collectedStrips.end(),true) <<std::endl;
0221       iEvent.put(std::move(removedStripClusterMask));
0222     }
0223     if (!pixelClusters_.isUninitialized()) {
0224       auto removedPixelClusterMask = std::make_unique<PixelMaskContainer>(
0225           edm::RefProd<edmNew::DetSetVector<SiPixelCluster>>(pixelClusters), collectedPixels);
0226       LogDebug("TrackClusterRemover") << "total pxl to skip: "
0227                                       << std::count(collectedPixels.begin(), collectedPixels.end(), true);
0228       iEvent.put(std::move(removedPixelClusterMask));
0229     }
0230   }
0231
0232 }  // namespace
0233
0234 #include "FWCore/PluginManager/interface/ModuleDef.h"
0235 #include "FWCore/Framework/interface/MakerMacros.h"
0236 DEFINE_FWK_MODULE(TrackClusterRemover);