NearbyPixelClustersProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    Calibration/TkAlCaRecoProducers
// Class:      NearbyPixelClustersProducer
//
//
// Original Author:  Marco Musich
//         Created:  Mon, 29 Mar 2021 12:29:30 GMT
//
//

// system include files
#include <memory>
#include <map>

// user include files
#include "DataFormats/Common/interface/DetSetVector.h"
#include "DataFormats/Common/interface/DetSetVectorNew.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
#include "DataFormats/TrackerCommon/interface/PixelEndcapName.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/StreamID.h"
#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
#include "Geometry/CommonTopologies/interface/PixelTopology.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
#include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
#include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
#include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
#include "CondFormats/DataRecord/interface/SiPixelQualityFromDbRcd.h"

using trajCrossings_t = std::map<uint32_t, std::vector<LocalPoint>>;

//
// class declaration
//

class NearbyPixelClustersProducer : public edm::stream::EDProducer<> {
public:
  explicit NearbyPixelClustersProducer(const edm::ParameterSet&);
  ~NearbyPixelClustersProducer() override = default;

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
  void produce(edm::Event&, const edm::EventSetup&) override;
  const trajCrossings_t findAllTrajectoriesCrossings(
      const edm::Handle<TrajTrackAssociationCollection>& trajTrackCollectionHandle);

  const std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> findAllNearbyClusters(
      const SiPixelClusterCollectionNew::const_iterator& clusterSet,
      const uint32_t rawId,
      const std::vector<LocalPoint>& vLocalPos);

  const std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> findAllNearbyClusters(
      const SiPixelClusterCollectionNew& clusterSet, const uint32_t rawId, const std::vector<LocalPoint>& vLocalPos);

  TrajectoryStateOnSurface getTrajectoryStateOnSurface(const TrajectoryMeasurement& measurement);
  bool detidIsOnPixel(const DetId& detid);

  // ----------member data ---------------------------

  // switches
  const bool throwBadComponents_;
  const bool dumpWholeDetId_;

  // esTokens
  const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
  const edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> pixelCPEEsToken_;
  const edm::ESGetToken<SiPixelQuality, SiPixelQualityFromDbRcd> badModuleToken_;

  // edToken
  edm::EDGetTokenT<SiPixelClusterCollectionNew> clustersToken_;
  edm::EDGetTokenT<TrajTrackAssociationCollection> trajTrackCollectionToken_;

  // putToken
  edm::EDPutTokenT<SiPixelClusterCollectionNew> clusterPutToken_;

  // event setup
  const TrackerGeometry* trackerGeometry_;
  const PixelClusterParameterEstimator* pixelCPE_;
};

//
// constructors and destructor
//
NearbyPixelClustersProducer::NearbyPixelClustersProducer(const edm::ParameterSet& iConfig)
    : throwBadComponents_(iConfig.getParameter<bool>("throwBadComponents")),
      dumpWholeDetId_(iConfig.getParameter<bool>("dumpWholeDetIds")),
      geomEsToken_(esConsumes()),
      pixelCPEEsToken_(esConsumes(edm::ESInputTag("", "PixelCPEGeneric"))),
      badModuleToken_(esConsumes()),
      clustersToken_(consumes<SiPixelClusterCollectionNew>(iConfig.getParameter<edm::InputTag>("clusterCollection"))),
      trajTrackCollectionToken_(
          consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("trajectoryInput"))),
      clusterPutToken_(produces<SiPixelClusterCollectionNew>()) {
  if (dumpWholeDetId_) {
    edm::LogPrint("NearbyPixelClustersProducer") << "WARNING: selected to dump the whole DetId's worth of clusters.\n "
                                                    "This will have consequences on the event size!"
                                                 << std::endl;
  }
}

//
// member functions
//

// ------------ method called to produce the data  ------------
void NearbyPixelClustersProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  using namespace edm;

  auto outputClusters = std::make_unique<SiPixelClusterCollectionNew>();

  // get the Tracker geometry from event setup
  trackerGeometry_ = &iSetup.getData(geomEsToken_);

  // get the Pixel CPE from event setup
  pixelCPE_ = &iSetup.getData(pixelCPEEsToken_);

  const auto& SiPixelBadModule_ = &iSetup.getData(badModuleToken_);

  // get cluster collection
  const auto& clusterCollectionHandle = iEvent.getHandle(clustersToken_);
  const SiPixelClusterCollectionNew& clusterCollection = *clusterCollectionHandle;

  // get Traj-Track Collection
  const auto& trajTrackCollectionHandle = iEvent.getHandle(trajTrackCollectionToken_);
  if (!trajTrackCollectionHandle.isValid())
    return;

  // find all trajectory crossings in the event
  const auto& allCrossings = this->findAllTrajectoriesCrossings(trajTrackCollectionHandle);

  LogDebug("NearbyPixelClustersProducer") << allCrossings.size() << std::endl;

  // now find all nearby clusters
  for (const auto& [id, vLocalPos] : allCrossings) {
    // prepare the filler
    SiPixelClusterCollectionNew::FastFiller spc(*outputClusters, id);

    // retrieve the clusters of the right detId
    const auto& clustersOnDet = clusterCollection.find(DetId(id));

    if (clustersOnDet == clusterCollection.end())
      continue;

    // if the cluster DetSet is not valid, move on
    if (!(*clustersOnDet).isValid())
      continue;

    // if the module is bad continue
    if (throwBadComponents_ && SiPixelBadModule_->IsModuleBad(id))
      continue;

    // find all the clusters to put into the event
    const auto& clustersToPut = this->findAllNearbyClusters(clustersOnDet, id, vLocalPos);

    // find all the clusters to put into the event (different interface)
    //const auto& clustersToPut = this->findAllNearbyClusters(clusterCollection, id, vLocalPos);

    for (const auto& cluster : clustersToPut) {
      spc.push_back(*cluster);
    }

    if (spc.empty())
      spc.abort();

  }  // loop on trajectory crossings

  iEvent.put(clusterPutToken_, std::move(outputClusters));
}

/*--------------------------------------------------------------------*/
const trajCrossings_t NearbyPixelClustersProducer::findAllTrajectoriesCrossings(
    const edm::Handle<TrajTrackAssociationCollection>& trajTrackCollectionHandle)
/*--------------------------------------------------------------------*/
{
  trajCrossings_t crossings;

  std::vector<uint32_t> treatedIds;

  for (const auto& pair : *trajTrackCollectionHandle) {
    const edm::Ref<std::vector<Trajectory>> traj = pair.key;

    for (const TrajectoryMeasurement& measurement : traj->measurements()) {
      //Check if the measurement infos can be read
      if (!measurement.updatedState().isValid())
        continue;

      const TransientTrackingRecHit::ConstRecHitPointer& recHit = measurement.recHit();

      // Only looking for pixel hits
      DetId recHitDetid = recHit->geographicalId();
      const auto& rawId = recHitDetid.rawId();

      if (!this->detidIsOnPixel(recHitDetid))
        continue;

      // Skipping hits with undeterminable positions
      TrajectoryStateOnSurface trajStateOnSurface = this->getTrajectoryStateOnSurface(measurement);

      if (!(trajStateOnSurface.isValid()))
        continue;

      // Position measurements
      // Looking for valid and missing hits
      LocalPoint localPosition = trajStateOnSurface.localPosition();

      if (std::find(treatedIds.begin(), treatedIds.end(), rawId) != treatedIds.end()) {
        crossings.at(rawId).push_back(localPosition);
      } else {
        crossings.insert(std::pair<uint32_t, std::vector<LocalPoint>>(rawId, {localPosition}));
        treatedIds.push_back(rawId);
      }
    }  // loop on measurements in trajectory
  }    // loop on trajectories

  return crossings;
}

/*--------------------------------------------------------------------*/
const std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> NearbyPixelClustersProducer::findAllNearbyClusters(
    const SiPixelClusterCollectionNew::const_iterator& clusterSet,
    const uint32_t rawId,
    const std::vector<LocalPoint>& vLocalPos)
/*--------------------------------------------------------------------*/
{
  std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> outputClusters;

  static constexpr unsigned int k_maxClustersInDet = 1024;

  // something funny is going on here ...
  if ((*clusterSet).size() > k_maxClustersInDet) {
    edm::LogWarning("NearbyPixelClustersProducer")
        << __func__ << "() number of clusters in det " << rawId /*(*clusterSet).id()*/ << " is " << (*clusterSet).size()
        << ", which is larger than maximum (1024).\n Something funny with the data is going on!" << std::endl;
    return outputClusters;
  }

  const PixelGeomDetUnit* pixdet = (const PixelGeomDetUnit*)trackerGeometry_->idToDetUnit(rawId);
  edmNew::DetSet<SiPixelCluster>::const_iterator itCluster = clusterSet->begin();

  // just copy straight the whole set of clusters in the detid
  if (dumpWholeDetId_) {
    for (; itCluster != clusterSet->end(); ++itCluster) {
      outputClusters.push_back(itCluster);
    }
    return outputClusters;
  }

  int count = 0;
  for (const auto& localPos : vLocalPos) {
    count++;
    //trajectory crossing local coordinates
    const auto& traj_lx = localPos.x();
    const auto& traj_ly = localPos.y();

    float minD = 10000.;
    edmNew::DetSet<SiPixelCluster>::const_iterator closest = nullptr;

    //std::cout<< rawId << " count: " << count << " n. clusters: " << (*clusterSet).size() << std::endl;
    LogDebug("NearbyPixelClustersProducer")
        << __func__ << rawId << " count: " << count << " n. clusters: " << (*clusterSet).size() << std::endl;

    for (; itCluster != clusterSet->end(); ++itCluster) {
      LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
      const auto& params = pixelCPE_->getParameters(*itCluster, *pixdet);
      lp = std::get<0>(params);

      float D = sqrt((lp.x() - traj_lx) * (lp.x() - traj_lx) + (lp.y() - traj_ly) * (lp.y() - traj_ly));
      if (D < minD) {
        closest = itCluster;
        minD = D;
      }
    }  // loop on cluster sets

    if (closest) {
      outputClusters.push_back(closest);
    }
  }  // loop on the trajectory crossings

  return outputClusters;
}

// overloaded method: use the whole DetSet
/*--------------------------------------------------------------------*/
const std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> NearbyPixelClustersProducer::findAllNearbyClusters(
    const SiPixelClusterCollectionNew& clusterCollection,
    const uint32_t rawId,
    const std::vector<LocalPoint>& vLocalPos)
/*--------------------------------------------------------------------*/
{
  std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> outputClusters;

  int count = 0;
  for (const auto& localPos : vLocalPos) {
    count++;

    //trajectory crossing local coordinates
    const auto& traj_lx = localPos.x();
    const auto& traj_ly = localPos.y();

    float minD = 10000.;

    SiPixelClusterCollectionNew::const_iterator itClusterSet = clusterCollection.begin();
    for (; itClusterSet != clusterCollection.end(); itClusterSet++) {
      DetId detId(itClusterSet->id());
      if (detId.rawId() != rawId)
        continue;

      unsigned int subDetId = detId.subdetId();
      if (subDetId != PixelSubdetector::PixelBarrel && subDetId != PixelSubdetector::PixelEndcap) {
        edm::LogError("NearByPixelClusterProducer")
            << "ERROR: not a pixel cluster!!!" << std::endl;  // should not happen
        continue;
      }

      edmNew::DetSet<SiPixelCluster>::const_iterator itCluster = itClusterSet->begin();

      // just copy straight the whole set of clusters in the detid
      if (dumpWholeDetId_ && count == 1) {
        for (; itCluster != itClusterSet->end(); ++itCluster) {
          outputClusters.push_back(itCluster);
        }
        return outputClusters;
      }

      //std::cout<< rawId << " count: " << count << " n. clusters: " << (*clusterSet).size() << std::endl;
      LogDebug("NearbyPixelClustersProducer")
          << __func__ << rawId << " count: " << count << " n. clusters: " << (*itClusterSet).size() << std::endl;

      const PixelGeomDetUnit* pixdet = (const PixelGeomDetUnit*)trackerGeometry_->idToDetUnit(rawId);

      edmNew::DetSet<SiPixelCluster>::const_iterator closest = nullptr;

      for (; itCluster != itClusterSet->end(); ++itCluster) {
        LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
        const auto& params = pixelCPE_->getParameters(*itCluster, *pixdet);
        lp = std::get<0>(params);

        float D = sqrt((lp.x() - traj_lx) * (lp.x() - traj_lx) + (lp.y() - traj_ly) * (lp.y() - traj_ly));
        if (D < minD) {
          closest = itCluster;
          minD = D;
        }
      }  // loop on cluster sets

      if (closest) {
        outputClusters.push_back(closest);
      }
    }  // loop on all clusters
  }    // loop on the trajectory crossings

  return outputClusters;
}

/*--------------------------------------------------------------------*/
bool NearbyPixelClustersProducer::detidIsOnPixel(const DetId& detid)
/*--------------------------------------------------------------------*/
{
  if (detid.det() != DetId::Tracker)
    return false;
  if (detid.subdetId() == PixelSubdetector::PixelBarrel)
    return true;
  if (detid.subdetId() == PixelSubdetector::PixelEndcap)
    return true;
  return false;
}

/*--------------------------------------------------------------------*/
TrajectoryStateOnSurface NearbyPixelClustersProducer::getTrajectoryStateOnSurface(
    const TrajectoryMeasurement& measurement)
/*--------------------------------------------------------------------*/
{
  const static TrajectoryStateCombiner trajStateCombiner;

  const auto& forwardPredictedState = measurement.forwardPredictedState();
  const auto& backwardPredictedState = measurement.backwardPredictedState();

  if (forwardPredictedState.isValid() && backwardPredictedState.isValid())
    return trajStateCombiner(forwardPredictedState, backwardPredictedState);

  else if (backwardPredictedState.isValid())
    return backwardPredictedState;

  else if (forwardPredictedState.isValid())
    return forwardPredictedState;

  edm::LogError("NearbyPixelClustersProducer") << "Error saving traj. measurement data."
                                               << " Trajectory state on surface cannot be determined." << std::endl;

  return TrajectoryStateOnSurface();
}

// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void NearbyPixelClustersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.setComment(
      "Produces the collection of SiPixelClusters closest, hit by hit, to the trajectory measurements of a given "
      "track");
  desc.add<bool>("throwBadComponents", false)
      ->setComment(
          "do not consider modules flagged as bad components. Careful, it changes the efficiency denominator!");
  desc.add<bool>("dumpWholeDetIds", false)
      ->setComment("put in the event all the pixel cluster on the impacted module, by default only the closest");
  ;
  desc.add<edm::InputTag>("clusterCollection", edm::InputTag("siPixelClusters"));
  desc.add<edm::InputTag>("trajectoryInput", edm::InputTag("myRefitter"));
  descriptions.addWithDefaultLabel(desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(NearbyPixelClustersProducer);