TkAlCaRecoProducers/plugins/NearbyPixelClustersProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Calibration/TkAlCaRecoProducers
0004 // Class:      NearbyPixelClustersProducer
0005 //
0006 /**\class NearbyPixelClustersProducer NearbyPixelClustersProducer.cc Calibration/TkAlCaRecoProducers/plugins/NearbyPixelClustersProducer.cc
0007
0008  Description: Class to produce the collection of SiPixelClusters closest, hit by hit, to the trajectory measurements of a given track
0009
0010  Implementation:
0011      Implementation of this class is heavily endebted to https://github.com/jkarancs/PhaseIPixelNtuplizer/blob/master/plugins/PhaseIPixelNtuplizer.h
0012
0013 */
0014 //
0015 // Original Author:  Marco Musich
0016 //         Created:  Mon, 29 Mar 2021 12:29:30 GMT
0017 //
0018 //
0019
0020 // system include files
0021 #include <memory>
0022 #include <map>
0023
0024 // user include files
0025 #include "DataFormats/Common/interface/DetSetVector.h"
0026 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0027 #include "DataFormats/DetId/interface/DetId.h"
0028 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0029 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0030 #include "DataFormats/TrackReco/interface/Track.h"
0031 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0032 #include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
0033 #include "DataFormats/TrackerCommon/interface/PixelEndcapName.h"
0034 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0035 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "FWCore/Framework/interface/Frameworkfwd.h"
0038 #include "FWCore/Framework/interface/MakerMacros.h"
0039 #include "FWCore/Framework/interface/stream/EDProducer.h"
0040 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0041 #include "FWCore/Utilities/interface/StreamID.h"
0042 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0043 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0044 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0045 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0046 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0047 #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
0048 #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
0049 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0050 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0051 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0052 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0053 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0054 #include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
0055 #include "CondFormats/DataRecord/interface/SiPixelQualityFromDbRcd.h"
0056
0057 using trajCrossings_t = std::map<uint32_t, std::vector<LocalPoint>>;
0058
0059 //
0060 // class declaration
0061 //
0062
0063 class NearbyPixelClustersProducer : public edm::stream::EDProducer<> {
0064 public:
0065   explicit NearbyPixelClustersProducer(const edm::ParameterSet&);
0066   ~NearbyPixelClustersProducer() override = default;
0067
0068   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0069
0070 private:
0071   void produce(edm::Event&, const edm::EventSetup&) override;
0072   const trajCrossings_t findAllTrajectoriesCrossings(
0073       const edm::Handle<TrajTrackAssociationCollection>& trajTrackCollectionHandle);
0074
0075   const std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> findAllNearbyClusters(
0076       const SiPixelClusterCollectionNew::const_iterator& clusterSet,
0077       const uint32_t rawId,
0078       const std::vector<LocalPoint>& vLocalPos);
0079
0080   const std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> findAllNearbyClusters(
0081       const SiPixelClusterCollectionNew& clusterSet, const uint32_t rawId, const std::vector<LocalPoint>& vLocalPos);
0082
0083   TrajectoryStateOnSurface getTrajectoryStateOnSurface(const TrajectoryMeasurement& measurement);
0084   bool detidIsOnPixel(const DetId& detid);
0085
0086   // ----------member data ---------------------------
0087
0088   // switches
0089   const bool throwBadComponents_;
0090   const bool dumpWholeDetId_;
0091
0092   // esTokens
0093   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
0094   const edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> pixelCPEEsToken_;
0095   const edm::ESGetToken<SiPixelQuality, SiPixelQualityFromDbRcd> badModuleToken_;
0096
0097   // edToken
0098   edm::EDGetTokenT<SiPixelClusterCollectionNew> clustersToken_;
0099   edm::EDGetTokenT<TrajTrackAssociationCollection> trajTrackCollectionToken_;
0100
0101   // putToken
0102   edm::EDPutTokenT<SiPixelClusterCollectionNew> clusterPutToken_;
0103
0104   // event setup
0105   const TrackerGeometry* trackerGeometry_;
0106   const PixelClusterParameterEstimator* pixelCPE_;
0107 };
0108
0109 //
0110 // constructors and destructor
0111 //
0112 NearbyPixelClustersProducer::NearbyPixelClustersProducer(const edm::ParameterSet& iConfig)
0113     : throwBadComponents_(iConfig.getParameter<bool>("throwBadComponents")),
0114       dumpWholeDetId_(iConfig.getParameter<bool>("dumpWholeDetIds")),
0115       geomEsToken_(esConsumes()),
0116       pixelCPEEsToken_(esConsumes(edm::ESInputTag("", "PixelCPEGeneric"))),
0117       badModuleToken_(esConsumes()),
0118       clustersToken_(consumes<SiPixelClusterCollectionNew>(iConfig.getParameter<edm::InputTag>("clusterCollection"))),
0119       trajTrackCollectionToken_(
0120           consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("trajectoryInput"))),
0121       clusterPutToken_(produces<SiPixelClusterCollectionNew>()) {
0122   if (dumpWholeDetId_) {
0123     edm::LogPrint("NearbyPixelClustersProducer") << "WARNING: selected to dump the whole DetId's worth of clusters.\n "
0124                                                     "This will have consequences on the event size!"
0125                                                  << std::endl;
0126   }
0127 }
0128
0129 //
0130 // member functions
0131 //
0132
0133 // ------------ method called to produce the data  ------------
0134 void NearbyPixelClustersProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0135   using namespace edm;
0136
0137   auto outputClusters = std::make_unique<SiPixelClusterCollectionNew>();
0138
0139   // get the Tracker geometry from event setup
0140   trackerGeometry_ = &iSetup.getData(geomEsToken_);
0141
0142   // get the Pixel CPE from event setup
0143   pixelCPE_ = &iSetup.getData(pixelCPEEsToken_);
0144
0145   const auto& SiPixelBadModule_ = &iSetup.getData(badModuleToken_);
0146
0147   // get cluster collection
0148   const auto& clusterCollectionHandle = iEvent.getHandle(clustersToken_);
0149   const SiPixelClusterCollectionNew& clusterCollection = *clusterCollectionHandle;
0150
0151   // get Traj-Track Collection
0152   const auto& trajTrackCollectionHandle = iEvent.getHandle(trajTrackCollectionToken_);
0153   if (!trajTrackCollectionHandle.isValid())
0154     return;
0155
0156   // find all trajectory crossings in the event
0157   const auto& allCrossings = this->findAllTrajectoriesCrossings(trajTrackCollectionHandle);
0158
0159   LogDebug("NearbyPixelClustersProducer") << allCrossings.size() << std::endl;
0160
0161   // now find all nearby clusters
0162   for (const auto& [id, vLocalPos] : allCrossings) {
0163     // prepare the filler
0164     SiPixelClusterCollectionNew::FastFiller spc(*outputClusters, id);
0165
0166     // retrieve the clusters of the right detId
0167     const auto& clustersOnDet = clusterCollection.find(DetId(id));
0168
0169     if (clustersOnDet == clusterCollection.end())
0170       continue;
0171
0172     // if the cluster DetSet is not valid, move on
0173     if (!(*clustersOnDet).isValid())
0174       continue;
0175
0176     // if the module is bad continue
0177     if (throwBadComponents_ && SiPixelBadModule_->IsModuleBad(id))
0178       continue;
0179
0180     // find all the clusters to put into the event
0181     const auto& clustersToPut = this->findAllNearbyClusters(clustersOnDet, id, vLocalPos);
0182
0183     // find all the clusters to put into the event (different interface)
0184     //const auto& clustersToPut = this->findAllNearbyClusters(clusterCollection, id, vLocalPos);
0185
0186     for (const auto& cluster : clustersToPut) {
0187       spc.push_back(*cluster);
0188     }
0189
0190     if (spc.empty())
0191       spc.abort();
0192
0193   }  // loop on trajectory crossings
0194
0195   iEvent.put(clusterPutToken_, std::move(outputClusters));
0196 }
0197
0198 /*--------------------------------------------------------------------*/
0199 const trajCrossings_t NearbyPixelClustersProducer::findAllTrajectoriesCrossings(
0200     const edm::Handle<TrajTrackAssociationCollection>& trajTrackCollectionHandle)
0201 /*--------------------------------------------------------------------*/
0202 {
0203   trajCrossings_t crossings;
0204
0205   std::vector<uint32_t> treatedIds;
0206
0207   for (const auto& pair : *trajTrackCollectionHandle) {
0208     const edm::Ref<std::vector<Trajectory>> traj = pair.key;
0209
0210     for (const TrajectoryMeasurement& measurement : traj->measurements()) {
0211       //Check if the measurement infos can be read
0212       if (!measurement.updatedState().isValid())
0213         continue;
0214
0215       const TransientTrackingRecHit::ConstRecHitPointer& recHit = measurement.recHit();
0216
0217       // Only looking for pixel hits
0218       DetId recHitDetid = recHit->geographicalId();
0219       const auto& rawId = recHitDetid.rawId();
0220
0221       if (!this->detidIsOnPixel(recHitDetid))
0222         continue;
0223
0224       // Skipping hits with undeterminable positions
0225       TrajectoryStateOnSurface trajStateOnSurface = this->getTrajectoryStateOnSurface(measurement);
0226
0227       if (!(trajStateOnSurface.isValid()))
0228         continue;
0229
0230       // Position measurements
0231       // Looking for valid and missing hits
0232       LocalPoint localPosition = trajStateOnSurface.localPosition();
0233
0234       if (std::find(treatedIds.begin(), treatedIds.end(), rawId) != treatedIds.end()) {
0235         crossings.at(rawId).push_back(localPosition);
0236       } else {
0237         crossings.insert(std::pair<uint32_t, std::vector<LocalPoint>>(rawId, {localPosition}));
0238         treatedIds.push_back(rawId);
0239       }
0240     }  // loop on measurements in trajectory
0241   }    // loop on trajectories
0242
0243   return crossings;
0244 }
0245
0246 /*--------------------------------------------------------------------*/
0247 const std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> NearbyPixelClustersProducer::findAllNearbyClusters(
0248     const SiPixelClusterCollectionNew::const_iterator& clusterSet,
0249     const uint32_t rawId,
0250     const std::vector<LocalPoint>& vLocalPos)
0251 /*--------------------------------------------------------------------*/
0252 {
0253   std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> outputClusters;
0254
0255   static constexpr unsigned int k_maxClustersInDet = 1024;
0256
0257   // something funny is going on here ...
0258   if ((*clusterSet).size() > k_maxClustersInDet) {
0259     edm::LogWarning("NearbyPixelClustersProducer")
0260         << __func__ << "() number of clusters in det " << rawId /*(*clusterSet).id()*/ << " is " << (*clusterSet).size()
0261         << ", which is larger than maximum (1024).\n Something funny with the data is going on!" << std::endl;
0262     return outputClusters;
0263   }
0264
0265   const PixelGeomDetUnit* pixdet = (const PixelGeomDetUnit*)trackerGeometry_->idToDetUnit(rawId);
0266   edmNew::DetSet<SiPixelCluster>::const_iterator itCluster = clusterSet->begin();
0267
0268   // just copy straight the whole set of clusters in the detid
0269   if (dumpWholeDetId_) {
0270     for (; itCluster != clusterSet->end(); ++itCluster) {
0271       outputClusters.push_back(itCluster);
0272     }
0273     return outputClusters;
0274   }
0275
0276   int count = 0;
0277   for (const auto& localPos : vLocalPos) {
0278     count++;
0279     //trajectory crossing local coordinates
0280     const auto& traj_lx = localPos.x();
0281     const auto& traj_ly = localPos.y();
0282
0283     float minD = 10000.;
0284     edmNew::DetSet<SiPixelCluster>::const_iterator closest = nullptr;
0285
0286     //std::cout<< rawId << " count: " << count << " n. clusters: " << (*clusterSet).size() << std::endl;
0287     LogDebug("NearbyPixelClustersProducer")
0288         << __func__ << rawId << " count: " << count << " n. clusters: " << (*clusterSet).size() << std::endl;
0289
0290     for (; itCluster != clusterSet->end(); ++itCluster) {
0291       LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
0292       const auto& params = pixelCPE_->getParameters(*itCluster, *pixdet);
0293       lp = std::get<0>(params);
0294
0295       float D = sqrt((lp.x() - traj_lx) * (lp.x() - traj_lx) + (lp.y() - traj_ly) * (lp.y() - traj_ly));
0296       if (D < minD) {
0297         closest = itCluster;
0298         minD = D;
0299       }
0300     }  // loop on cluster sets
0301
0302     if (closest) {
0303       outputClusters.push_back(closest);
0304     }
0305   }  // loop on the trajectory crossings
0306
0307   return outputClusters;
0308 }
0309
0310 // overloaded method: use the whole DetSet
0311 /*--------------------------------------------------------------------*/
0312 const std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> NearbyPixelClustersProducer::findAllNearbyClusters(
0313     const SiPixelClusterCollectionNew& clusterCollection,
0314     const uint32_t rawId,
0315     const std::vector<LocalPoint>& vLocalPos)
0316 /*--------------------------------------------------------------------*/
0317 {
0318   std::vector<edmNew::DetSet<SiPixelCluster>::const_iterator> outputClusters;
0319
0320   int count = 0;
0321   for (const auto& localPos : vLocalPos) {
0322     count++;
0323
0324     //trajectory crossing local coordinates
0325     const auto& traj_lx = localPos.x();
0326     const auto& traj_ly = localPos.y();
0327
0328     float minD = 10000.;
0329
0330     SiPixelClusterCollectionNew::const_iterator itClusterSet = clusterCollection.begin();
0331     for (; itClusterSet != clusterCollection.end(); itClusterSet++) {
0332       DetId detId(itClusterSet->id());
0333       if (detId.rawId() != rawId)
0334         continue;
0335
0336       unsigned int subDetId = detId.subdetId();
0337       if (subDetId != PixelSubdetector::PixelBarrel && subDetId != PixelSubdetector::PixelEndcap) {
0338         edm::LogError("NearByPixelClusterProducer")
0339             << "ERROR: not a pixel cluster!!!" << std::endl;  // should not happen
0340         continue;
0341       }
0342
0343       edmNew::DetSet<SiPixelCluster>::const_iterator itCluster = itClusterSet->begin();
0344
0345       // just copy straight the whole set of clusters in the detid
0346       if (dumpWholeDetId_ && count == 1) {
0347         for (; itCluster != itClusterSet->end(); ++itCluster) {
0348           outputClusters.push_back(itCluster);
0349         }
0350         return outputClusters;
0351       }
0352
0353       //std::cout<< rawId << " count: " << count << " n. clusters: " << (*clusterSet).size() << std::endl;
0354       LogDebug("NearbyPixelClustersProducer")
0355           << __func__ << rawId << " count: " << count << " n. clusters: " << (*itClusterSet).size() << std::endl;
0356
0357       const PixelGeomDetUnit* pixdet = (const PixelGeomDetUnit*)trackerGeometry_->idToDetUnit(rawId);
0358
0359       edmNew::DetSet<SiPixelCluster>::const_iterator closest = nullptr;
0360
0361       for (; itCluster != itClusterSet->end(); ++itCluster) {
0362         LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
0363         const auto& params = pixelCPE_->getParameters(*itCluster, *pixdet);
0364         lp = std::get<0>(params);
0365
0366         float D = sqrt((lp.x() - traj_lx) * (lp.x() - traj_lx) + (lp.y() - traj_ly) * (lp.y() - traj_ly));
0367         if (D < minD) {
0368           closest = itCluster;
0369           minD = D;
0370         }
0371       }  // loop on cluster sets
0372
0373       if (closest) {
0374         outputClusters.push_back(closest);
0375       }
0376     }  // loop on all clusters
0377   }    // loop on the trajectory crossings
0378
0379   return outputClusters;
0380 }
0381
0382 /*--------------------------------------------------------------------*/
0383 bool NearbyPixelClustersProducer::detidIsOnPixel(const DetId& detid)
0384 /*--------------------------------------------------------------------*/
0385 {
0386   if (detid.det() != DetId::Tracker)
0387     return false;
0388   if (detid.subdetId() == PixelSubdetector::PixelBarrel)
0389     return true;
0390   if (detid.subdetId() == PixelSubdetector::PixelEndcap)
0391     return true;
0392   return false;
0393 }
0394
0395 /*--------------------------------------------------------------------*/
0396 TrajectoryStateOnSurface NearbyPixelClustersProducer::getTrajectoryStateOnSurface(
0397     const TrajectoryMeasurement& measurement)
0398 /*--------------------------------------------------------------------*/
0399 {
0400   const static TrajectoryStateCombiner trajStateCombiner;
0401
0402   const auto& forwardPredictedState = measurement.forwardPredictedState();
0403   const auto& backwardPredictedState = measurement.backwardPredictedState();
0404
0405   if (forwardPredictedState.isValid() && backwardPredictedState.isValid())
0406     return trajStateCombiner(forwardPredictedState, backwardPredictedState);
0407
0408   else if (backwardPredictedState.isValid())
0409     return backwardPredictedState;
0410
0411   else if (forwardPredictedState.isValid())
0412     return forwardPredictedState;
0413
0414   edm::LogError("NearbyPixelClustersProducer") << "Error saving traj. measurement data."
0415                                                << " Trajectory state on surface cannot be determined." << std::endl;
0416
0417   return TrajectoryStateOnSurface();
0418 }
0419
0420 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0421 void NearbyPixelClustersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0422   edm::ParameterSetDescription desc;
0423   desc.setComment(
0424       "Produces the collection of SiPixelClusters closest, hit by hit, to the trajectory measurements of a given "
0425       "track");
0426   desc.add<bool>("throwBadComponents", false)
0427       ->setComment(
0428           "do not consider modules flagged as bad components. Careful, it changes the efficiency denominator!");
0429   desc.add<bool>("dumpWholeDetIds", false)
0430       ->setComment("put in the event all the pixel cluster on the impacted module, by default only the closest");
0431   ;
0432   desc.add<edm::InputTag>("clusterCollection", edm::InputTag("siPixelClusters"));
0433   desc.add<edm::InputTag>("trajectoryInput", edm::InputTag("myRefitter"));
0434   descriptions.addWithDefaultLabel(desc);
0435 }
0436
0437 //define this as a plug-in
0438 DEFINE_FWK_MODULE(NearbyPixelClustersProducer);