CalibTracker/plugins/SiPixelCalSingleMuonAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    DQMOffline/CalibTracker
0004 // Class:      SiPixelCalSingleMuonAnalyzer
0005 //
0006 /**\class SiPixelCalSingleMuonAnalyzer SiPixelCalSingleMuonAnalyzer.cc DQMOffline/CalibTracker/plugins/SiPixelCalSingleMuonAnalyzer.cc
0007
0008  Description: Analysis of the closebyPixelClusters collections for Pixel Hit Efficiency mearurements purposes
0009 */
0010 //
0011 // Original Author:  Marco Musich
0012 //         Created:  Tue, 30 Mar 2021 09:22:07 GMT
0013 //
0014 //
0015
0016 // system include files
0017 #include <memory>
0018
0019 // user include files
0020 #include "CalibTracker/SiPixelESProducers/interface/SiPixelDetInfoFileReader.h"
0021 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0022 #include "DataFormats/Common/interface/DetSetVector.h"
0023 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0024 #include "DataFormats/DetId/interface/DetId.h"
0025 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0026 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0027 #include "DataFormats/TrackReco/interface/Track.h"
0028 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0029 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0030 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "FWCore/Framework/interface/Frameworkfwd.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0036 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0037 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0038 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0039 #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
0040 #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
0041 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0042 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0043 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0044 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0045 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0046
0047 //
0048 // class declaration
0049 //
0050 class SiPixelCalSingleMuonAnalyzer : public DQMEDAnalyzer {
0051 public:
0052   explicit SiPixelCalSingleMuonAnalyzer(const edm::ParameterSet&);
0053   ~SiPixelCalSingleMuonAnalyzer() override = default;
0054   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0055
0056 private:
0057   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0058   void analyze(const edm::Event&, const edm::EventSetup&) override;
0059   void countClusters(const edm::Handle<SiPixelClusterCollectionNew>& handle, unsigned int& nClusGlobal);
0060   const bool detidIsOnPixel(const DetId& detid);
0061
0062   TrajectoryStateOnSurface getTrajectoryStateOnSurface(const TrajectoryMeasurement& measurement);
0063   std::pair<float, float> findClosestCluster(const edm::Handle<SiPixelClusterCollectionNew>& handle,
0064                                              const PixelClusterParameterEstimator* pixelCPE_,
0065                                              const TrackerGeometry* trackerGeometry_,
0066                                              uint32_t rawId,
0067                                              float traj_lx,
0068                                              float traj_ly);
0069
0070   // ----------member data ---------------------------
0071   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
0072   const edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> pixelCPEEsToken_;
0073   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsTokenBR_;
0074   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoEsTokenBR_;
0075
0076   const edm::EDGetTokenT<SiPixelClusterCollectionNew> clustersToken_;
0077   const edm::EDGetTokenT<SiPixelClusterCollectionNew> nearByClustersToken_;
0078   const edm::EDGetTokenT<TrajTrackAssociationCollection> trajTrackCollectionToken_;
0079   const edm::EDGetTokenT<edm::ValueMap<std::vector<float>>> distanceToken_;
0080   const edm::EDGetTokenT<edm::View<reco::Track>> muonTracksToken_;
0081
0082   const std::string dqm_path_;
0083   SiPixelDetInfoFileReader reader_;
0084
0085   typedef dqm::reco::MonitorElement MonitorElement;
0086   MonitorElement* h_nALCARECOClusters;
0087   MonitorElement* h_nCloseByClusters;
0088   MonitorElement* h_distClosestValid;
0089   MonitorElement* h_distClosestMissing;
0090   MonitorElement* h_distClosestInactive;
0091   MonitorElement* h_distClosestTrack;
0092   bool phase_;
0093 };
0094
0095 //
0096 // constructors and destructor
0097 //
0098 SiPixelCalSingleMuonAnalyzer::SiPixelCalSingleMuonAnalyzer(const edm::ParameterSet& iConfig)
0099     : geomEsToken_(esConsumes()),
0100       pixelCPEEsToken_(esConsumes(edm::ESInputTag("", "PixelCPEGeneric"))),
0101       geomEsTokenBR_(esConsumes<edm::Transition::BeginRun>()),
0102       topoEsTokenBR_(esConsumes<edm::Transition::BeginRun>()),
0103       clustersToken_(consumes<SiPixelClusterCollectionNew>(iConfig.getParameter<edm::InputTag>("clusterCollection"))),
0104       nearByClustersToken_(
0105           consumes<SiPixelClusterCollectionNew>(iConfig.getParameter<edm::InputTag>("nearByClusterCollection"))),
0106       trajTrackCollectionToken_(
0107           consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("trajectoryInput"))),
0108       distanceToken_(consumes<edm::ValueMap<std::vector<float>>>(iConfig.getParameter<edm::InputTag>("distToTrack"))),
0109       muonTracksToken_(consumes<edm::View<reco::Track>>(iConfig.getParameter<edm::InputTag>("muonTracks"))),
0110       dqm_path_(iConfig.getParameter<std::string>("dqmPath")),
0111       reader_(edm::FileInPath(iConfig.getParameter<std::string>("skimmedGeometryPath")).fullPath()) {}
0112
0113 //
0114 // member functions
0115 //
0116
0117 // ------------ method called for each event  ------------
0118 void SiPixelCalSingleMuonAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0119   using namespace edm;
0120
0121   // get the Tracker geometry from event setup
0122   const TrackerGeometry* trackerGeometry_ = &iSetup.getData(geomEsToken_);
0123
0124   // get the Pixel CPE from event setup
0125   const PixelClusterParameterEstimator* pixelCPE_ = &iSetup.getData(pixelCPEEsToken_);
0126
0127   // get the muon track collection
0128   const auto& muonTrackCollectionHandle = iEvent.getHandle(muonTracksToken_);
0129   auto const& muonTracks = *muonTrackCollectionHandle;
0130
0131   // get the track distances
0132   const auto& distancesToTrack = iEvent.getHandle(distanceToken_);
0133
0134   unsigned int nMuons = muonTracks.size();
0135   for (unsigned int ij = 0; ij < nMuons; ij++) {
0136     auto muon = muonTrackCollectionHandle->ptrAt(ij);
0137     edm::RefToBase<reco::Track> trackRef = muonTrackCollectionHandle->refAt(ij);
0138     const auto& distances = (*distancesToTrack)[trackRef];
0139
0140     LogDebug("SiPixelCalSingleMuonAnalyzer") << "distances size: " << distances.size() << std::endl;
0141
0142     unsigned counter = 0;
0143     double closestDR = 999.;
0144     for (const auto& distance : distances) {
0145       counter++;
0146       LogDebug("SiPixelCalSingleMuonAnalyzer")
0147           << "track: " << counter << " distance:" << std::sqrt(distance) << std::endl;
0148       if (distance < closestDR && distance > 0) {
0149         closestDR = distance;
0150       }
0151     }
0152
0153     h_distClosestTrack->Fill(std::sqrt(closestDR));
0154   }
0155
0156   // Get cluster collection
0157   const auto& clusterCollectionHandle = iEvent.getHandle(clustersToken_);
0158
0159   unsigned int nClusGlobal = 0;
0160   countClusters(clusterCollectionHandle, nClusGlobal);
0161
0162   h_nALCARECOClusters->Fill(nClusGlobal);
0163   LogTrace("SiPixelCalSingleMuonAnalyzer") << "total ALCARECO clusters: " << nClusGlobal << std::endl;
0164
0165   // Get nearby cluster collection
0166   const auto& nearByClusterCollectionHandle = iEvent.getHandle(nearByClustersToken_);
0167
0168   unsigned int nNearByClusGlobal = 0;
0169   countClusters(nearByClusterCollectionHandle, nNearByClusGlobal);
0170
0171   h_nCloseByClusters->Fill(nNearByClusGlobal);
0172   LogTrace("SiPixelCalSingleMuonAnalyzer") << "total close-by clusters: " << nNearByClusGlobal << std::endl;
0173
0174   // Get Traj-Track Collection
0175   const auto& trajTrackCollectionHandle = iEvent.getHandle(trajTrackCollectionToken_);
0176
0177   if (!trajTrackCollectionHandle.isValid())
0178     return;
0179
0180   for (const auto& pair : *trajTrackCollectionHandle) {
0181     const edm::Ref<std::vector<Trajectory>> traj = pair.key;
0182     const reco::TrackRef track = pair.val;
0183
0184     for (const TrajectoryMeasurement& measurement : traj->measurements()) {
0185       if (!measurement.updatedState().isValid())
0186         return;
0187
0188       const TransientTrackingRecHit::ConstRecHitPointer& recHit = measurement.recHit();
0189
0190       // Only looking for pixel hits
0191       DetId r_rawId = recHit->geographicalId();
0192
0193       if (!this->detidIsOnPixel(r_rawId))
0194         continue;
0195
0196       // Skipping hits with undeterminable positions
0197       TrajectoryStateOnSurface trajStateOnSurface = this->getTrajectoryStateOnSurface(measurement);
0198
0199       if (!(trajStateOnSurface.isValid()))
0200         continue;
0201
0202       // Position measurements
0203       // Looking for valid and missing hits
0204       LocalPoint localPosition = trajStateOnSurface.localPosition();
0205
0206       const auto& traj_lx = localPosition.x();
0207       const auto& traj_ly = localPosition.y();
0208
0209       const auto loc = this->findClosestCluster(
0210           nearByClusterCollectionHandle, pixelCPE_, trackerGeometry_, r_rawId.rawId(), traj_lx, traj_ly);
0211
0212       float dist = (loc.first != -999.) ? std::sqrt(loc.first * loc.first + loc.second * loc.second) : -0.1;
0213
0214       if (recHit->getType() == TrackingRecHit::valid) {
0215         LogTrace("SiPixelCalSingleMuonAnalyzer")
0216             << "RawID:" << r_rawId.rawId() << " (valid hit), distance: " << dist << std::endl;
0217         h_distClosestValid->Fill(dist);
0218       }
0219
0220       if (recHit->getType() == TrackingRecHit::missing) {
0221         LogTrace("SiPixelCalSingleMuonAnalyzer")
0222             << "RawID:" << r_rawId.rawId() << " (missing hit), distance: " << dist << std::endl;
0223         h_distClosestMissing->Fill(dist);
0224       }
0225
0226       if (recHit->getType() == TrackingRecHit::inactive) {
0227         LogTrace("SiPixelCalSingleMuonAnalyzer")
0228             << "RawID:" << r_rawId.rawId() << " (inactive hit), distance: " << dist << std::endl;
0229         h_distClosestInactive->Fill(dist);
0230       }
0231     }
0232   }
0233 }
0234
0235 // ------------ method called once each job just before starting event loop  ------------
0236 void SiPixelCalSingleMuonAnalyzer::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const&, edm::EventSetup const&) {
0237   // book the overall event count and event types histograms
0238   iBooker.setCurrentFolder(dqm_path_ + "/ClusterCounts");
0239   h_nALCARECOClusters = iBooker.book1I(
0240       "h_nALCARECOClusters", "Number of Pixel clusters per event (ALCARECO) ;N_{clusters};events", 20, 0, 20);
0241   h_nCloseByClusters = iBooker.book1I(
0242       "h_nCloseByClusters", "Number of Pixel clusters per event (close-by) ;N_{clusters};events", 20, 0, 20);
0243
0244   iBooker.setCurrentFolder(dqm_path_ + "/TrajDistance");
0245   h_distClosestValid =
0246       iBooker.book1D("h_distClosestValid",
0247                      "Distance of Closest cluster to trajectory (valid);distance (cm); valid trajectory measurements",
0248                      110,
0249                      -0.105,
0250                      0.995);
0251   h_distClosestMissing = iBooker.book1D(
0252       "h_distClosestMissing",
0253       "Distance of Closest cluster to trajectory (missing);distance (cm);missing trajectory measurements",
0254       110,
0255       -0.105,
0256       0.995);
0257   h_distClosestInactive = iBooker.book1D(
0258       "h_distClosestInactive",
0259       "Distance of Closest cluster to trajectory (inactive);distance (cm);inactive trajectory measurements",
0260       110,
0261       -0.105,
0262       0.995);
0263
0264   iBooker.setCurrentFolder(dqm_path_ + "/OtherTrackDistance");
0265   h_distClosestTrack =
0266       iBooker.book1D("h_distClosestTrack",
0267                      "#DeltaR Distance of Closest track to the muon trajectory;#DeltaR distance; muon trajectories",
0268                      100,
0269                      0.,
0270                      5.);
0271 }
0272
0273 /*--------------------------------------------------------------------*/
0274 const bool SiPixelCalSingleMuonAnalyzer::detidIsOnPixel(const DetId& detid)
0275 /*--------------------------------------------------------------------*/
0276 {
0277   if (detid.det() != DetId::Tracker)
0278     return false;
0279   if (detid.subdetId() == PixelSubdetector::PixelBarrel)
0280     return true;
0281   if (detid.subdetId() == PixelSubdetector::PixelEndcap)
0282     return true;
0283   return false;
0284 }
0285
0286 /*--------------------------------------------------------------------*/
0287 TrajectoryStateOnSurface SiPixelCalSingleMuonAnalyzer::getTrajectoryStateOnSurface(
0288     const TrajectoryMeasurement& measurement)
0289 /*--------------------------------------------------------------------*/
0290 {
0291   const static TrajectoryStateCombiner trajStateCombiner;
0292
0293   const auto& forwardPredictedState = measurement.forwardPredictedState();
0294   const auto& backwardPredictedState = measurement.backwardPredictedState();
0295
0296   if (forwardPredictedState.isValid() && backwardPredictedState.isValid())
0297     return trajStateCombiner(forwardPredictedState, backwardPredictedState);
0298
0299   else if (backwardPredictedState.isValid())
0300     return backwardPredictedState;
0301
0302   else if (forwardPredictedState.isValid())
0303     return forwardPredictedState;
0304
0305   edm::LogError("NearbyPixelClusterProducer") << "Error saving traj. measurement data."
0306                                               << " Trajectory state on surface cannot be determined." << std::endl;
0307
0308   return TrajectoryStateOnSurface();
0309 }
0310
0311 /*--------------------------------------------------------------------*/
0312 void SiPixelCalSingleMuonAnalyzer::countClusters(const edm::Handle<SiPixelClusterCollectionNew>& handle,
0313                                                  //const TrackerGeometry* tkGeom_,
0314                                                  unsigned int& nClusGlobal)
0315 /*--------------------------------------------------------------------*/
0316 {
0317   for (const auto& DSVItr : *handle) {
0318     uint32_t rawid(DSVItr.detId());
0319     DetId detId(rawid);
0320     LogDebug("SiPixelCalSingleMuonAnalyzer") << "DetId: " << detId.rawId() << " size: " << DSVItr.size() << std::endl;
0321     nClusGlobal += DSVItr.size();
0322   }
0323 }
0324
0325 /*--------------------------------------------------------------------*/
0326 std::pair<float, float> SiPixelCalSingleMuonAnalyzer::findClosestCluster(
0327     const edm::Handle<SiPixelClusterCollectionNew>& handle,
0328     const PixelClusterParameterEstimator* pixelCPE_,
0329     const TrackerGeometry* trackerGeometry_,
0330     uint32_t rawId,
0331     float traj_lx,
0332     float traj_ly)
0333 /*--------------------------------------------------------------------*/
0334 {
0335   const SiPixelClusterCollectionNew& clusterCollection = *handle;
0336   SiPixelClusterCollectionNew::const_iterator itClusterSet = clusterCollection.begin();
0337
0338   float minD = 10e7;
0339
0340   auto loc = std::make_pair(-999., -999.);
0341
0342   for (; itClusterSet != clusterCollection.end(); itClusterSet++) {
0343     DetId detId(itClusterSet->id());
0344     if (detId.rawId() != rawId)
0345       continue;
0346
0347     unsigned int subDetId = detId.subdetId();
0348     if (subDetId != PixelSubdetector::PixelBarrel && subDetId != PixelSubdetector::PixelEndcap) {
0349       edm::LogError("NearByPixelClustersAnalyzer")
0350           << "ERROR: not a pixel cluster!!!" << std::endl;  // should not happen
0351       continue;
0352     }
0353
0354     const PixelGeomDetUnit* pixdet = (const PixelGeomDetUnit*)trackerGeometry_->idToDetUnit(detId);
0355     edmNew::DetSet<SiPixelCluster>::const_iterator itCluster = itClusterSet->begin();
0356     for (; itCluster != itClusterSet->end(); ++itCluster) {
0357       LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
0358       PixelClusterParameterEstimator::ReturnType params = pixelCPE_->getParameters(*itCluster, *pixdet);
0359       lp = std::get<0>(params);
0360
0361       float D = (lp.x() - traj_lx) * (lp.x() - traj_lx) + (lp.y() - traj_ly) * (lp.y() - traj_ly);
0362       if (D < minD) {
0363         minD = D;
0364         loc.first = (lp.x() - traj_lx);
0365         loc.second = (lp.y() - traj_ly);
0366       }
0367     }  // loop on cluster sets
0368   }
0369   return loc;
0370 }
0371
0372 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0373 void SiPixelCalSingleMuonAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0374   edm::ParameterSetDescription desc;
0375   desc.setComment("Analysis of the closebyPixelClusters collections");
0376   desc.add<edm::InputTag>("clusterCollection", edm::InputTag("ALCARECOSiPixelCalSingleMuonTight"));
0377   desc.add<edm::InputTag>("nearByClusterCollection", edm::InputTag("closebyPixelClusters"));
0378   desc.add<edm::InputTag>("trajectoryInput", edm::InputTag("refittedTracks"));
0379   desc.add<edm::InputTag>("muonTracks", edm::InputTag("ALCARECOSiPixelCalSingleMuonTight"));
0380   desc.add<edm::InputTag>("distToTrack", edm::InputTag("trackDistances"));
0381   desc.add<std::string>("dqmPath", "SiPixelCalSingleMuonTight");
0382   desc.add<std::string>("skimmedGeometryPath",
0383                         "SLHCUpgradeSimulations/Geometry/data/PhaseI/PixelSkimmedGeometry_phase1.txt");
0384   descriptions.addWithDefaultLabel(desc);
0385 }
0386
0387 //define this as a plug-in
0388 DEFINE_FWK_MODULE(SiPixelCalSingleMuonAnalyzer);