Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:11:21

0001 // Package:    SiPixelMonitorTrack
0002 // Class:      SiPixelHitEfficiencySource
0003 //
0004 // class SiPixelHitEfficiencyModule SiPixelHitEfficiencyModule.cc
0005 //       DQM/SiPixelMonitorTrack/src/SiPixelHitEfficiencyModule.cc
0006 //
0007 // Description: SiPixel hit efficiency data quality monitoring modules
0008 // Implementation: prototype -> improved -> never final - end of the 1st step
0009 //
0010 // Original Authors: Romain Rougny & Luca Mucibello
0011 //         Created: Mar Nov 10 13:29:00 CET 2009
0012 
0013 #include <cmath>
0014 #include <iostream>
0015 #include <map>
0016 #include <string>
0017 #include <utility>
0018 #include <vector>
0019 
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 
0023 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0024 #include "DataFormats/DetId/interface/DetId.h"
0025 #include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
0026 #include "DataFormats/SiPixelDetId/interface/PixelBarrelNameUpgrade.h"
0027 #include "DataFormats/TrackerCommon/interface/PixelEndcapName.h"
0028 #include "DataFormats/SiPixelDetId/interface/PixelEndcapNameUpgrade.h"
0029 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0030 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0031 #include "DataFormats/VertexReco/interface/Vertex.h"
0032 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0033 
0034 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0035 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0036 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0037 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0038 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0039 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0040 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0041 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0042 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0043 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0044 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0045 
0046 #include "DQM/SiPixelCommon/interface/SiPixelFolderOrganizer.h"
0047 #include "DQM/SiPixelMonitorTrack/interface/SiPixelHitEfficiencySource.h"
0048 
0049 using namespace std;
0050 using namespace edm;
0051 
0052 SiPixelHitEfficiencySource::SiPixelHitEfficiencySource(const edm::ParameterSet &pSet)
0053     : pSet_(pSet),
0054       modOn(pSet.getUntrackedParameter<bool>("modOn", true)),
0055       ladOn(pSet.getUntrackedParameter<bool>("ladOn", false)),
0056       layOn(pSet.getUntrackedParameter<bool>("layOn", false)),
0057       phiOn(pSet.getUntrackedParameter<bool>("phiOn", false)),
0058       ringOn(pSet.getUntrackedParameter<bool>("ringOn", false)),
0059       bladeOn(pSet.getUntrackedParameter<bool>("bladeOn", false)),
0060       diskOn(pSet.getUntrackedParameter<bool>("diskOn", false)),
0061       isUpgrade(pSet.getUntrackedParameter<bool>("isUpgrade", false))
0062 // updateEfficiencies(
0063 // pSet.getUntrackedParameter<bool>("updateEfficiencies",false) )
0064 {
0065   pSet_ = pSet;
0066   debug_ = pSet_.getUntrackedParameter<bool>("debug", false);
0067   applyEdgeCut_ = pSet_.getUntrackedParameter<bool>("applyEdgeCut");
0068   nSigma_EdgeCut_ = pSet_.getUntrackedParameter<double>("nSigma_EdgeCut");
0069   vtxsrc_ = pSet_.getUntrackedParameter<std::string>("vtxsrc", "offlinePrimaryVertices");
0070   vertexCollectionToken_ = consumes<reco::VertexCollection>(vtxsrc_);
0071   tracksrc_ = consumes<TrajTrackAssociationCollection>(pSet_.getParameter<edm::InputTag>("trajectoryInput"));
0072   clusterCollectionToken_ = consumes<edmNew::DetSetVector<SiPixelCluster>>(std::string("siPixelClusters"));
0073 
0074   measurementTrackerEventToken_ = consumes<MeasurementTrackerEvent>(std::string("MeasurementTrackerEvent"));
0075 
0076   trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0077   trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0078   measurementTrackerToken_ = esConsumes<MeasurementTracker, CkfComponentsRecord>();
0079   chi2MeasurementEstimatorBaseToken_ =
0080       esConsumes<Chi2MeasurementEstimatorBase, TrackingComponentsRecord>(edm::ESInputTag("", "Chi2"));
0081   propagatorToken_ = esConsumes<Propagator, TrackingComponentsRecord>(edm::ESInputTag("", "PropagatorWithMaterial"));
0082   pixelClusterParameterEstimatorToken_ =
0083       esConsumes<PixelClusterParameterEstimator, TkPixelCPERecord>(edm::ESInputTag("", "PixelCPEGeneric"));
0084   trackerTopoTokenBeginRun_ = esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>();
0085   trackerGeomTokenBeginRun_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>();
0086 
0087   firstRun = true;
0088 
0089   LogInfo("PixelDQM") << "SiPixelHitEfficiencySource constructor" << endl;
0090   LogInfo("PixelDQM") << "Mod/Lad/Lay/Phi " << modOn << "/" << ladOn << "/" << layOn << "/" << phiOn << std::endl;
0091   LogInfo("PixelDQM") << "Blade/Disk/Ring" << bladeOn << "/" << diskOn << "/" << ringOn << std::endl;
0092 }
0093 
0094 SiPixelHitEfficiencySource::~SiPixelHitEfficiencySource() {
0095   LogInfo("PixelDQM") << "SiPixelHitEfficiencySource destructor" << endl;
0096 
0097   std::map<uint32_t, SiPixelHitEfficiencyModule *>::iterator struct_iter;
0098   for (struct_iter = theSiPixelStructure.begin(); struct_iter != theSiPixelStructure.end(); struct_iter++) {
0099     delete struct_iter->second;
0100     struct_iter->second = nullptr;
0101   }
0102 }
0103 
0104 void SiPixelHitEfficiencySource::fillClusterProbability(int layer, int disk, bool plus, double probability) {
0105   // barrel
0106   if (layer != 0) {
0107     if (layer == 1) {
0108       if (plus)
0109         meClusterProbabilityL1_Plus_->Fill(probability);
0110       else
0111         meClusterProbabilityL1_Minus_->Fill(probability);
0112     }
0113 
0114     else if (layer == 2) {
0115       if (plus)
0116         meClusterProbabilityL2_Plus_->Fill(probability);
0117       else
0118         meClusterProbabilityL2_Minus_->Fill(probability);
0119     }
0120 
0121     else if (layer == 3) {
0122       if (plus)
0123         meClusterProbabilityL3_Plus_->Fill(probability);
0124       else
0125         meClusterProbabilityL3_Minus_->Fill(probability);
0126     }
0127   }
0128   // Endcap
0129   if (disk != 0) {
0130     if (disk == 1) {
0131       if (plus)
0132         meClusterProbabilityD1_Plus_->Fill(probability);
0133       else
0134         meClusterProbabilityD1_Minus_->Fill(probability);
0135     }
0136     if (disk == 2) {
0137       if (plus)
0138         meClusterProbabilityD2_Plus_->Fill(probability);
0139       else
0140         meClusterProbabilityD2_Minus_->Fill(probability);
0141     }
0142   }
0143 }
0144 
0145 void SiPixelHitEfficiencySource::dqmBeginRun(const edm::Run &r, edm::EventSetup const &iSetup) {
0146   LogInfo("PixelDQM") << "SiPixelHitEfficiencySource beginRun()" << endl;
0147 
0148   if (firstRun) {
0149     // retrieve TrackerGeometry for pixel dets
0150 
0151     nvalid = 0;
0152     nmissing = 0;
0153 
0154     firstRun = false;
0155   }
0156 
0157   edm::ESHandle<TrackerGeometry> TG = iSetup.getHandle(trackerGeomTokenBeginRun_);
0158   if (debug_)
0159     LogVerbatim("PixelDQM") << "TrackerGeometry " << &(*TG) << " size is " << TG->dets().size() << endl;
0160 
0161   // build theSiPixelStructure with the pixel barrel and endcap dets from
0162   // TrackerGeometry
0163   for (TrackerGeometry::DetContainer::const_iterator pxb = TG->detsPXB().begin(); pxb != TG->detsPXB().end(); pxb++) {
0164     if (dynamic_cast<PixelGeomDetUnit const *>((*pxb)) != nullptr) {
0165       SiPixelHitEfficiencyModule *module = new SiPixelHitEfficiencyModule((*pxb)->geographicalId().rawId());
0166       theSiPixelStructure.insert(
0167           pair<uint32_t, SiPixelHitEfficiencyModule *>((*pxb)->geographicalId().rawId(), module));
0168     }
0169   }
0170   for (TrackerGeometry::DetContainer::const_iterator pxf = TG->detsPXF().begin(); pxf != TG->detsPXF().end(); pxf++) {
0171     if (dynamic_cast<PixelGeomDetUnit const *>((*pxf)) != nullptr) {
0172       SiPixelHitEfficiencyModule *module = new SiPixelHitEfficiencyModule((*pxf)->geographicalId().rawId());
0173       theSiPixelStructure.insert(
0174           pair<uint32_t, SiPixelHitEfficiencyModule *>((*pxf)->geographicalId().rawId(), module));
0175     }
0176   }
0177   LogInfo("PixelDQM") << "SiPixelStructure size is " << theSiPixelStructure.size() << endl;
0178 }
0179 
0180 void SiPixelHitEfficiencySource::bookHistograms(DQMStore::IBooker &iBooker,
0181                                                 edm::Run const &iRun,
0182                                                 edm::EventSetup const &iSetup) {
0183   // book residual histograms in theSiPixelFolder - one (x,y) pair of histograms
0184   // per det
0185   SiPixelFolderOrganizer theSiPixelFolder(false);
0186 
0187   edm::ESHandle<TrackerTopology> tTopoHandle = iSetup.getHandle(trackerTopoTokenBeginRun_);
0188   const TrackerTopology *pTT = tTopoHandle.product();
0189 
0190   for (std::map<uint32_t, SiPixelHitEfficiencyModule *>::iterator pxd = theSiPixelStructure.begin();
0191        pxd != theSiPixelStructure.end();
0192        pxd++) {
0193     if (modOn) {
0194       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 0, isUpgrade))
0195         (*pxd).second->book(pSet_, pTT, iBooker, 0, isUpgrade);
0196       else
0197         throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Folder Creation Failed! ";
0198     }
0199     if (ladOn) {
0200       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 1, isUpgrade))
0201         (*pxd).second->book(pSet_, pTT, iBooker, 1, isUpgrade);
0202       else
0203         throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource ladder Folder Creation Failed! ";
0204     }
0205     if (layOn) {
0206       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 2, isUpgrade))
0207         (*pxd).second->book(pSet_, pTT, iBooker, 2, isUpgrade);
0208       else
0209         throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource layer Folder Creation Failed! ";
0210     }
0211     if (phiOn) {
0212       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 3, isUpgrade))
0213         (*pxd).second->book(pSet_, pTT, iBooker, 3, isUpgrade);
0214       else
0215         throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource phi Folder Creation Failed! ";
0216     }
0217     if (bladeOn) {
0218       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 4, isUpgrade))
0219         (*pxd).second->book(pSet_, pTT, iBooker, 4, isUpgrade);
0220       else
0221         throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Blade Folder Creation Failed! ";
0222     }
0223     if (diskOn) {
0224       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 5, isUpgrade))
0225         (*pxd).second->book(pSet_, pTT, iBooker, 5, isUpgrade);
0226       else
0227         throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Disk Folder Creation Failed! ";
0228     }
0229     if (ringOn) {
0230       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 6, isUpgrade))
0231         (*pxd).second->book(pSet_, pTT, iBooker, 6, isUpgrade);
0232       else
0233         throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Ring Folder Creation Failed! ";
0234     }
0235   }
0236 
0237   // book cluster probability histos for Barrel and Endcap
0238   iBooker.setCurrentFolder("Pixel/Barrel");
0239 
0240   meClusterProbabilityL1_Plus_ =
0241       iBooker.book1D("ClusterProbabilityXY_Layer1_Plus", "ClusterProbabilityXY_Layer1_Plus", 500, -14, 0.1);
0242   meClusterProbabilityL1_Plus_->setAxisTitle("Log(ClusterProbability)", 1);
0243 
0244   meClusterProbabilityL1_Minus_ =
0245       iBooker.book1D("ClusterProbabilityXY_Layer1_Minus", "ClusterProbabilityXY_Layer1_Minus", 500, -14, 0.1);
0246   meClusterProbabilityL1_Minus_->setAxisTitle("Log(ClusterProbability)", 1);
0247 
0248   meClusterProbabilityL2_Plus_ =
0249       iBooker.book1D("ClusterProbabilityXY_Layer2_Plus", "ClusterProbabilityXY_Layer2_Plus", 500, -14, 0.1);
0250   meClusterProbabilityL2_Plus_->setAxisTitle("Log(ClusterProbability)", 1);
0251 
0252   meClusterProbabilityL2_Minus_ =
0253       iBooker.book1D("ClusterProbabilityXY_Layer2_Minus", "ClusterProbabilityXY_Layer2_Minus", 500, -14, 0.1);
0254   meClusterProbabilityL2_Minus_->setAxisTitle("Log(ClusterProbability)", 1);
0255 
0256   meClusterProbabilityL3_Plus_ =
0257       iBooker.book1D("ClusterProbabilityXY_Layer3_Plus", "ClusterProbabilityXY_Layer3_Plus", 500, -14, 0.1);
0258   meClusterProbabilityL3_Plus_->setAxisTitle("Log(ClusterProbability)", 1);
0259 
0260   meClusterProbabilityL3_Minus_ =
0261       iBooker.book1D("ClusterProbabilityXY_Layer3_Minus", "ClusterProbabilityXY_Layer3_Minus", 500, -14, 0.1);
0262   meClusterProbabilityL3_Minus_->setAxisTitle("Log(ClusterProbability)", 1);
0263 
0264   iBooker.setCurrentFolder("Pixel/Endcap");
0265 
0266   meClusterProbabilityD1_Plus_ =
0267       iBooker.book1D("ClusterProbabilityXY_Disk1_Plus", "ClusterProbabilityXY_Disk1_Plus", 500, -14, 0.1);
0268   meClusterProbabilityD1_Plus_->setAxisTitle("Log(ClusterProbability)", 1);
0269 
0270   meClusterProbabilityD1_Minus_ =
0271       iBooker.book1D("ClusterProbabilityXY_Disk1_Minus", "ClusterProbabilityXY_Disk1_Minus", 500, -14, 0.1);
0272   meClusterProbabilityD1_Minus_->setAxisTitle("Log(ClusterProbability)", 1);
0273 
0274   meClusterProbabilityD2_Plus_ =
0275       iBooker.book1D("ClusterProbabilityXY_Disk2_Plus", "ClusterProbabilityXY_Disk2_Plus", 500, -14, 0.1);
0276   meClusterProbabilityD2_Plus_->setAxisTitle("Log(ClusterProbability)", 1);
0277 
0278   meClusterProbabilityD2_Minus_ =
0279       iBooker.book1D("ClusterProbabilityXY_Disk2_Minus", "ClusterProbabilityXY_Disk2_Minus", 500, -14, 0.1);
0280   meClusterProbabilityD2_Minus_->setAxisTitle("Log(ClusterProbability)", 1);
0281 }
0282 
0283 void SiPixelHitEfficiencySource::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0284   edm::ESHandle<TrackerTopology> tTopoHandle = iSetup.getHandle(trackerTopoToken_);
0285   const TrackerTopology *pTT = tTopoHandle.product();
0286 
0287   edm::Handle<reco::VertexCollection> vertexCollectionHandle;
0288   iEvent.getByToken(vertexCollectionToken_, vertexCollectionHandle);
0289   if (!vertexCollectionHandle.isValid())
0290     return;
0291   nvtx_ = 0;
0292   vtxntrk_ = -9999;
0293   vtxD0_ = -9999.;
0294   vtxX_ = -9999.;
0295   vtxY_ = -9999.;
0296   vtxZ_ = -9999.;
0297   vtxndof_ = -9999.;
0298   vtxchi2_ = -9999.;
0299   const reco::VertexCollection &vertices = *vertexCollectionHandle.product();
0300   reco::VertexCollection::const_iterator bestVtx = vertices.end();
0301   for (reco::VertexCollection::const_iterator it = vertices.begin(); it != vertices.end(); ++it) {
0302     if (!it->isValid())
0303       continue;
0304     if (vtxntrk_ == -9999 || vtxntrk_ < int(it->tracksSize()) ||
0305         (vtxntrk_ == int(it->tracksSize()) && fabs(vtxZ_) > fabs(it->z()))) {
0306       vtxntrk_ = it->tracksSize();
0307       vtxD0_ = it->position().rho();
0308       vtxX_ = it->x();
0309       vtxY_ = it->y();
0310       vtxZ_ = it->z();
0311       vtxndof_ = it->ndof();
0312       vtxchi2_ = it->chi2();
0313       bestVtx = it;
0314     }
0315     if (fabs(it->z()) <= 20. && fabs(it->position().rho()) <= 2. && it->ndof() > 4)
0316       nvtx_++;
0317   }
0318   if (nvtx_ < 1)
0319     return;
0320 
0321   // get the map
0322   edm::Handle<TrajTrackAssociationCollection> match;
0323   iEvent.getByToken(tracksrc_, match);
0324 
0325   if (debug_) {
0326     edm::EDConsumerBase::Labels labels;
0327     labelsForToken(tracksrc_, labels);
0328     std::cout << "Trajectories from " << labels.module << std::endl;
0329   }
0330 
0331   if (!match.isValid())
0332     return;
0333 
0334   const TrajTrackAssociationCollection ttac = *(match.product());
0335 
0336   if (debug_) {
0337     std::cout << "+++ NEW EVENT +++" << std::endl;
0338     std::cout << "Map entries \t : " << ttac.size() << std::endl;
0339   }
0340 
0341   std::set<SiPixelCluster> clusterSet;
0342   //  TrajectoryStateCombiner tsoscomb;
0343   // define variables for extrapolation
0344   int extrapolateFrom_ = 2;
0345   int extrapolateTo_ = 1;
0346   float maxlxmatch_ = 0.2;
0347   float maxlymatch_ = 0.2;
0348   bool keepOriginalMissingHit_ = true;
0349   edm::ESHandle<MeasurementTracker> measurementTrackerHandle = iSetup.getHandle(measurementTrackerToken_);
0350 
0351   edm::ESHandle<Chi2MeasurementEstimatorBase> est = iSetup.getHandle(chi2MeasurementEstimatorBaseToken_);
0352   edm::Handle<MeasurementTrackerEvent> measurementTrackerEventHandle;
0353   iEvent.getByToken(measurementTrackerEventToken_, measurementTrackerEventHandle);
0354   edm::ESHandle<Propagator> prop = iSetup.getHandle(propagatorToken_);
0355   Propagator *thePropagator = prop.product()->clone();
0356   // determines direction of the propagator => inward
0357   if (extrapolateFrom_ >= extrapolateTo_) {
0358     thePropagator->setPropagationDirection(oppositeToMomentum);
0359   }
0360   TrajectoryStateCombiner trajStateComb;
0361   bool debug_ = false;
0362 
0363   // Loop over track collection
0364   for (TrajTrackAssociationCollection::const_iterator it = ttac.begin(); it != ttac.end(); ++it) {
0365     // define vector to save extrapolated tracks
0366     std::vector<TrajectoryMeasurement> expTrajMeasurements;
0367 
0368     const edm::Ref<std::vector<Trajectory>> traj_iterator = it->key;
0369     // Trajectory Map, extract Trajectory for this track
0370     reco::TrackRef trackref = it->val;
0371     // tracks++;
0372 
0373     bool isBpixtrack = false, isFpixtrack = false;
0374     int nStripHits = 0;
0375     int L1hits = 0;
0376     int L2hits = 0;
0377     int L3hits = 0;
0378     int L4hits = 0;
0379     int D1hits = 0;
0380     int D2hits = 0;
0381     int D3hits = 0;
0382     std::vector<TrajectoryMeasurement> tmeasColl = traj_iterator->measurements();
0383     std::vector<TrajectoryMeasurement>::const_iterator tmeasIt;
0384     // loop on measurements to find out what kind of hits there are
0385     for (tmeasIt = tmeasColl.begin(); tmeasIt != tmeasColl.end(); tmeasIt++) {
0386       // if(! tmeasIt->updatedState().isValid()) continue; NOT NECESSARY (I
0387       // HOPE)
0388       TransientTrackingRecHit::ConstRecHitPointer testhit = tmeasIt->recHit();
0389       if (testhit->geographicalId().det() != DetId::Tracker)
0390         continue;
0391       uint testSubDetID = (testhit->geographicalId().subdetId());
0392       const DetId &hit_detId = testhit->geographicalId();
0393       int hit_layer = 0;
0394       int hit_ladder = 0;
0395       int hit_mod = 0;
0396       int hit_disk = 0;
0397 
0398       if (testSubDetID == PixelSubdetector::PixelBarrel) {
0399         isBpixtrack = true;
0400         hit_layer = PixelBarrelName(hit_detId, pTT, isUpgrade).layerName();
0401 
0402         hit_ladder = PXBDetId(hit_detId).ladder();
0403         hit_mod = PXBDetId(hit_detId).module();
0404 
0405         if (hit_layer == 1)
0406           L1hits++;
0407         if (hit_layer == 2)
0408           L2hits++;
0409         if (hit_layer == 3)
0410           L3hits++;
0411         if (hit_layer == 4)
0412           L4hits++;
0413       }
0414       if (testSubDetID == PixelSubdetector::PixelEndcap) {
0415         isFpixtrack = true;
0416         hit_disk = PixelEndcapName(hit_detId, pTT, isUpgrade).diskName();
0417 
0418         if (hit_disk == 1)
0419           D1hits++;
0420         if (hit_disk == 2)
0421           D2hits++;
0422         if (hit_disk == 3)
0423           D3hits++;
0424       }
0425       if (testSubDetID == StripSubdetector::TIB)
0426         nStripHits++;
0427       if (testSubDetID == StripSubdetector::TOB)
0428         nStripHits++;
0429       if (testSubDetID == StripSubdetector::TID)
0430         nStripHits++;
0431       if (testSubDetID == StripSubdetector::TEC)
0432         nStripHits++;
0433       // check if last valid hit is in Layer 2 or Disk 1
0434 
0435       bool lastValidL2 = false;
0436       if ((testSubDetID == PixelSubdetector::PixelBarrel && hit_layer == extrapolateFrom_) ||
0437           (testSubDetID == PixelSubdetector::PixelEndcap && hit_disk == 1)) {
0438         if (testhit->isValid()) {
0439           if (tmeasIt == tmeasColl.end() - 1) {
0440             lastValidL2 = true;
0441           } else {
0442             tmeasIt++;
0443             TransientTrackingRecHit::ConstRecHitPointer nextRecHit = tmeasIt->recHit();
0444             uint nextSubDetID = (nextRecHit->geographicalId().subdetId());
0445             int nextlayer = PixelBarrelName(nextRecHit->geographicalId()).layerName();
0446             if (nextSubDetID == PixelSubdetector::PixelBarrel && nextlayer == extrapolateTo_) {
0447               lastValidL2 = true;  //&& !nextRecHit->isValid()) lastValidL2=true;
0448             }
0449             tmeasIt--;
0450           }
0451         }
0452       }  // end check last valid layer
0453       if (lastValidL2) {
0454         std::vector<const BarrelDetLayer *> pxbLayers =
0455             measurementTrackerHandle->geometricSearchTracker()->pixelBarrelLayers();
0456         const DetLayer *pxb1 = pxbLayers[extrapolateTo_ - 1];
0457         const MeasurementEstimator *estimator = est.product();
0458         const LayerMeasurements *theLayerMeasurements =
0459             new LayerMeasurements(*measurementTrackerHandle, *measurementTrackerEventHandle);
0460         const TrajectoryStateOnSurface tsosPXB2 = tmeasIt->updatedState();
0461         expTrajMeasurements = theLayerMeasurements->measurements(*pxb1, tsosPXB2, *thePropagator, *estimator);
0462         delete theLayerMeasurements;
0463         if (!expTrajMeasurements.empty()) {
0464           for (uint p = 0; p < expTrajMeasurements.size(); p++) {
0465             TrajectoryMeasurement pxb1TM(expTrajMeasurements[p]);
0466             const auto &pxb1Hit = pxb1TM.recHit();
0467             // remove hits with rawID == 0
0468             if (pxb1Hit->geographicalId().rawId() == 0) {
0469               expTrajMeasurements.erase(expTrajMeasurements.begin() + p);
0470               continue;
0471             }
0472           }
0473         }
0474         //
0475       }
0476       // check if extrapolated hit to layer 1 one matches the original hit
0477       TrajectoryStateOnSurface chkPredTrajState =
0478           trajStateComb(tmeasIt->forwardPredictedState(), tmeasIt->backwardPredictedState());
0479       if (!chkPredTrajState.isValid())
0480         continue;
0481       float chkx = chkPredTrajState.globalPosition().x();
0482       float chky = chkPredTrajState.globalPosition().y();
0483       float chkz = chkPredTrajState.globalPosition().z();
0484       LocalPoint chklp = chkPredTrajState.localPosition();
0485       if (testSubDetID == PixelSubdetector::PixelBarrel && hit_layer == extrapolateTo_) {
0486         // Here we will drop the extrapolated hits if there is a hit and use
0487         // that hit
0488         vector<int> imatches;
0489         size_t imatch = 0;
0490         float glmatch = 9999.;
0491         for (size_t iexp = 0; iexp < expTrajMeasurements.size(); iexp++) {
0492           const DetId &exphit_detId = expTrajMeasurements[iexp].recHit()->geographicalId();
0493           int exphit_ladder = PXBDetId(exphit_detId).ladder();
0494           int exphit_mod = PXBDetId(exphit_detId).module();
0495           int dladder = abs(exphit_ladder - hit_ladder);
0496           if (dladder > 10)
0497             dladder = 20 - dladder;
0498           int dmodule = abs(exphit_mod - hit_mod);
0499           if (dladder != 0 || dmodule != 0) {
0500             continue;
0501           }
0502 
0503           TrajectoryStateOnSurface predTrajState = expTrajMeasurements[iexp].updatedState();
0504           float x = predTrajState.globalPosition().x();
0505           float y = predTrajState.globalPosition().y();
0506           float z = predTrajState.globalPosition().z();
0507           float dxyz = sqrt((chkx - x) * (chkx - x) + (chky - y) * (chky - y) + (chkz - z) * (chkz - z));
0508 
0509           if (dxyz <= glmatch) {
0510             glmatch = dxyz;
0511             imatch = iexp;
0512             imatches.push_back(int(imatch));
0513           }
0514 
0515         }  // found the propagated traj best matching the hit in data
0516 
0517         float lxmatch = 9999.0;
0518         float lymatch = 9999.0;
0519         if (!expTrajMeasurements.empty()) {
0520           if (glmatch < 9999.) {  // if there is any propagated trajectory for this hit
0521             const DetId &matchhit_detId = expTrajMeasurements[imatch].recHit()->geographicalId();
0522 
0523             int matchhit_ladder = PXBDetId(matchhit_detId).ladder();
0524             int dladder = abs(matchhit_ladder - hit_ladder);
0525             if (dladder > 10)
0526               dladder = 20 - dladder;
0527             LocalPoint lp = expTrajMeasurements[imatch].updatedState().localPosition();
0528             lxmatch = fabs(lp.x() - chklp.x());
0529             lymatch = fabs(lp.y() - chklp.y());
0530           }
0531           if (lxmatch < maxlxmatch_ && lymatch < maxlymatch_) {
0532             if (testhit->getType() != TrackingRecHit::missing || keepOriginalMissingHit_) {
0533               expTrajMeasurements.erase(expTrajMeasurements.begin() + imatch);
0534             }
0535           }
0536 
0537         }  // expected trajectory measurment not empty
0538       }
0539     }  // loop on trajectory measurments tmeasColl
0540 
0541     // if an extrapolated hit was found but not matched to an exisitng L1 hit
0542     // then push the hit back into the collection now keep the first one that is
0543     // left
0544     if (!expTrajMeasurements.empty()) {
0545       for (size_t f = 0; f < expTrajMeasurements.size(); f++) {
0546         TrajectoryMeasurement AddHit = expTrajMeasurements[f];
0547         if (AddHit.recHit()->getType() == TrackingRecHit::missing) {
0548           tmeasColl.push_back(AddHit);
0549           isBpixtrack = true;
0550         }
0551       }
0552     }
0553 
0554     if (isBpixtrack || isFpixtrack) {
0555       if (trackref->pt() < 0.6 || nStripHits < 8 || fabs(trackref->dxy(bestVtx->position())) > 0.05 ||
0556           fabs(trackref->dz(bestVtx->position())) > 0.5)
0557         continue;
0558 
0559       if (debug_) {
0560         std::cout << "isBpixtrack : " << isBpixtrack << std::endl;
0561         std::cout << "isFpixtrack : " << isFpixtrack << std::endl;
0562       }
0563       // std::cout<<"This tracks has so many hits:
0564       // "<<tmeasColl.size()<<std::endl;
0565       for (std::vector<TrajectoryMeasurement>::const_iterator tmeasIt = tmeasColl.begin(); tmeasIt != tmeasColl.end();
0566            tmeasIt++) {
0567         // if(! tmeasIt->updatedState().isValid()) continue;
0568         TrajectoryStateOnSurface tsos = tmeasIt->updatedState();
0569 
0570         TransientTrackingRecHit::ConstRecHitPointer hit = tmeasIt->recHit();
0571         if (hit->geographicalId().det() != DetId::Tracker)
0572           continue;
0573         else {
0574           //      //residual
0575           const DetId &hit_detId = hit->geographicalId();
0576           // uint IntRawDetID = (hit_detId.rawId());
0577           uint IntSubDetID = (hit_detId.subdetId());
0578 
0579           if (IntSubDetID == 0) {
0580             if (debug_)
0581               std::cout << "NO IntSubDetID\n";
0582             continue;
0583           }
0584           if (IntSubDetID != PixelSubdetector::PixelBarrel && IntSubDetID != PixelSubdetector::PixelEndcap)
0585             continue;
0586 
0587           int disk = 0;
0588           int layer = 0;
0589           int panel = 0;
0590           int module = 0;
0591           bool isHalfModule = false;
0592 
0593           const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit *>(hit->hit());
0594 
0595           if (IntSubDetID == PixelSubdetector::PixelBarrel) {  // it's a BPIX hit
0596             layer = PixelBarrelName(hit_detId, pTT, isUpgrade).layerName();
0597             isHalfModule = PixelBarrelName(hit_detId, pTT, isUpgrade).isHalfModule();
0598 
0599             if (hit->isValid()) {  // fill the cluster probability in barrel
0600               bool plus = true;
0601               if ((PixelBarrelName(hit_detId, pTT, isUpgrade).shell() == PixelBarrelName::Shell::mO) ||
0602                   (PixelBarrelName(hit_detId, pTT, isUpgrade).shell() == PixelBarrelName::Shell::mI))
0603                 plus = false;
0604               double clusterProbability = pixhit->clusterProbability(0);
0605               if (clusterProbability != 0)
0606                 fillClusterProbability(layer, 0, plus, log10(clusterProbability));
0607             }
0608 
0609           } else if (IntSubDetID == PixelSubdetector::PixelEndcap) {  // it's an FPIX hit
0610             disk = PixelEndcapName(hit_detId, pTT, isUpgrade).diskName();
0611             panel = PixelEndcapName(hit_detId, pTT, isUpgrade).pannelName();
0612             module = PixelEndcapName(hit_detId, pTT, isUpgrade).plaquetteName();
0613 
0614             if (hit->isValid()) {
0615               bool plus = true;
0616               if ((PixelEndcapName(hit_detId, pTT, isUpgrade).halfCylinder() == PixelEndcapName::HalfCylinder::mO) ||
0617                   (PixelEndcapName(hit_detId, pTT, isUpgrade).halfCylinder() == PixelEndcapName::HalfCylinder::mI))
0618                 plus = false;
0619               double clusterProbability = pixhit->clusterProbability(0);
0620               if (clusterProbability != 0)
0621                 fillClusterProbability(0, disk, plus, log10(clusterProbability));
0622             }
0623           }
0624 
0625           if (layer == 1) {
0626             if (fabs(trackref->dxy(bestVtx->position())) > 0.01 || fabs(trackref->dz(bestVtx->position())) > 0.1)
0627               continue;
0628             if (!(L2hits > 0 && L3hits > 0) && !(L2hits > 0 && D1hits > 0) && !(D1hits > 0 && D2hits > 0))
0629               continue;
0630           } else if (layer == 2) {
0631             if (fabs(trackref->dxy(bestVtx->position())) > 0.02 || fabs(trackref->dz(bestVtx->position())) > 0.1)
0632               continue;
0633             if (!(L1hits > 0 && L3hits > 0) && !(L1hits > 0 && D1hits > 0))
0634               continue;
0635           } else if (layer == 3) {
0636             if (fabs(trackref->dxy(bestVtx->position())) > 0.02 || fabs(trackref->dz(bestVtx->position())) > 0.1)
0637               continue;
0638             if (!(L1hits > 0 && L2hits > 0))
0639               continue;
0640           } else if (layer == 4) {
0641             if (fabs(trackref->dxy(bestVtx->position())) > 0.02 || fabs(trackref->dz(bestVtx->position())) > 0.1)
0642               continue;
0643           } else if (disk == 1) {
0644             if (fabs(trackref->dxy(bestVtx->position())) > 0.05 || fabs(trackref->dz(bestVtx->position())) > 0.5)
0645               continue;
0646             if (!(L1hits > 0 && D2hits > 0) && !(L2hits > 0 && D2hits > 0))
0647               continue;
0648           } else if (disk == 2) {
0649             if (fabs(trackref->dxy(bestVtx->position())) > 0.05 || fabs(trackref->dz(bestVtx->position())) > 0.5)
0650               continue;
0651             if (!(L1hits > 0 && D1hits > 0))
0652               continue;
0653           } else if (disk == 3) {
0654             if (fabs(trackref->dxy(bestVtx->position())) > 0.05 || fabs(trackref->dz(bestVtx->position())) > 0.5)
0655               continue;
0656           }
0657 
0658           // check whether hit is valid or missing using track algo flag
0659           bool isHitValid = hit->hit()->getType() == TrackingRecHit::valid;
0660           bool isHitMissing = hit->hit()->getType() == TrackingRecHit::missing;
0661 
0662           if (debug_)
0663             std::cout << "the hit is persistent\n";
0664 
0665           std::map<uint32_t, SiPixelHitEfficiencyModule *>::iterator pxd =
0666               theSiPixelStructure.find((*hit).geographicalId().rawId());
0667 
0668           // calculate alpha and beta from cluster position
0669           LocalTrajectoryParameters ltp = tsos.localParameters();
0670           // LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
0671 
0672           //*************** Edge cut ********************
0673           double lx = tsos.localPosition().x();
0674           double ly = tsos.localPosition().y();
0675 
0676           if (fabs(lx) > 0.55 || fabs(ly) > 3.0)
0677             continue;
0678 
0679           bool passedFiducial = true;
0680 
0681           // Module fiducials:
0682           if (IntSubDetID == PixelSubdetector::PixelBarrel && fabs(ly) >= 3.1)
0683             passedFiducial = false;
0684           if (IntSubDetID == PixelSubdetector::PixelEndcap &&
0685               !((panel == 1 &&
0686                  ((module == 1 && fabs(ly) < 0.7) ||
0687                   ((module == 2 && fabs(ly) < 1.1) && !(disk == -1 && ly > 0.8 && lx > 0.2) &&
0688                    !(disk == 1 && ly < -0.7 && lx > 0.2) && !(disk == 2 && ly < -0.8)) ||
0689                   ((module == 3 && fabs(ly) < 1.5) && !(disk == -2 && lx > 0.1 && ly > 1.0) &&
0690                    !(disk == 2 && lx > 0.1 && ly < -1.0)) ||
0691                   ((module == 4 && fabs(ly) < 1.9) && !(disk == -2 && ly > 1.5) && !(disk == 2 && ly < -1.5)))) ||
0692                 (panel == 2 &&
0693                  ((module == 1 && fabs(ly) < 0.7) ||
0694                   (module == 2 && fabs(ly) < 1.2 && !(disk > 0 && ly > 1.1) && !(disk < 0 && ly < -1.1)) ||
0695                   (module == 3 && fabs(ly) < 1.6 && !(disk > 0 && ly > 1.5) && !(disk < 0 && ly < -1.5))))))
0696             passedFiducial = false;
0697           if (IntSubDetID == PixelSubdetector::PixelEndcap &&
0698               ((panel == 1 && (module == 1 || (module >= 3 && abs(disk) == 1))) ||
0699                (panel == 2 && ((module == 1 && abs(disk) == 2) || (module == 3 && abs(disk) == 1)))))
0700             passedFiducial = false;
0701           // ROC fiducials:
0702           double ly_mod = fabs(ly);
0703           if (IntSubDetID == PixelSubdetector::PixelEndcap && (panel + module) % 2 == 1)
0704             ly_mod = ly_mod + 0.405;
0705           float d_rocedge = fabs(fmod(ly_mod, 0.81) - 0.405);
0706           if (d_rocedge <= 0.0625)
0707             passedFiducial = false;
0708           if (!((IntSubDetID == PixelSubdetector::PixelBarrel &&
0709                  ((!isHalfModule && fabs(lx) < 0.6) || (isHalfModule && lx > -0.3 && lx < 0.2))) ||
0710                 (IntSubDetID == PixelSubdetector::PixelEndcap &&
0711                  ((panel == 1 &&
0712                    ((module == 1 && fabs(lx) < 0.2) ||
0713                     (module == 2 && ((fabs(lx) < 0.55 && abs(disk) == 1) || (lx > -0.5 && lx < 0.2 && disk == -2) ||
0714                                      (lx > -0.5 && lx < 0.0 && disk == 2))) ||
0715                     (module == 3 && lx > -0.6 && lx < 0.5) || (module == 4 && lx > -0.3 && lx < 0.15))) ||
0716                   (panel == 2 && ((module == 1 && fabs(lx) < 0.6) ||
0717                                   (module == 2 && ((fabs(lx) < 0.55 && abs(disk) == 1) ||
0718                                                    (lx > -0.6 && lx < 0.5 && abs(disk) == 2))) ||
0719                                   (module == 3 && fabs(lx) < 0.5)))))))
0720             passedFiducial = false;
0721           if (((IntSubDetID == PixelSubdetector::PixelBarrel && !isHalfModule) ||
0722                (IntSubDetID == PixelSubdetector::PixelEndcap && !(panel == 1 && (module == 1 || module == 4)))) &&
0723               fabs(lx) < 0.06)
0724             passedFiducial = false;
0725 
0726           //*************** find closest clusters ********************
0727           float dx_cl[2];
0728           float dy_cl[2];
0729           dx_cl[0] = dx_cl[1] = dy_cl[0] = dy_cl[1] = -9999.;
0730           edm::ESHandle<PixelClusterParameterEstimator> cpEstimator =
0731               iSetup.getHandle(pixelClusterParameterEstimatorToken_);
0732           if (cpEstimator.isValid()) {
0733             const PixelClusterParameterEstimator &cpe(*cpEstimator);
0734             edm::ESHandle<TrackerGeometry> tracker = iSetup.getHandle(trackerGeomToken_);
0735             if (tracker.isValid()) {
0736               const TrackerGeometry *tkgeom = &(*tracker);
0737               edm::Handle<edmNew::DetSetVector<SiPixelCluster>> clusterCollectionHandle;
0738               iEvent.getByToken(clusterCollectionToken_, clusterCollectionHandle);
0739               if (clusterCollectionHandle.isValid()) {
0740                 const edmNew::DetSetVector<SiPixelCluster> &clusterCollection = *clusterCollectionHandle;
0741                 edmNew::DetSetVector<SiPixelCluster>::const_iterator itClusterSet = clusterCollection.begin();
0742                 float minD[2];
0743                 minD[0] = minD[1] = 10000.;
0744                 for (; itClusterSet != clusterCollection.end(); itClusterSet++) {
0745                   DetId detId(itClusterSet->id());
0746                   if (detId.rawId() != hit->geographicalId().rawId())
0747                     continue;
0748                   // unsigned int sdId=detId.subdetId();
0749                   const PixelGeomDetUnit *pixdet = (const PixelGeomDetUnit *)tkgeom->idToDetUnit(detId);
0750                   edmNew::DetSet<SiPixelCluster>::const_iterator itCluster = itClusterSet->begin();
0751                   for (; itCluster != itClusterSet->end(); ++itCluster) {
0752                     LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
0753                     PixelClusterParameterEstimator::ReturnType params = cpe.getParameters(*itCluster, *pixdet);
0754                     lp = std::get<0>(params);
0755                     float D = sqrt((lp.x() - lx) * (lp.x() - lx) + (lp.y() - ly) * (lp.y() - ly));
0756                     if (D < minD[0]) {
0757                       minD[1] = minD[0];
0758                       dx_cl[1] = dx_cl[0];
0759                       dy_cl[1] = dy_cl[0];
0760                       minD[0] = D;
0761                       dx_cl[0] = lp.x();
0762                       dy_cl[0] = lp.y();
0763                     } else if (D < minD[1]) {
0764                       minD[1] = D;
0765                       dx_cl[1] = lp.x();
0766                       dy_cl[1] = lp.y();
0767                     }
0768                   }  // loop on clusterSets
0769                 }    // loop on clusterCollection
0770                 for (size_t i = 0; i < 2; i++) {
0771                   if (minD[i] < 9999.) {
0772                     dx_cl[i] = fabs(dx_cl[i] - lx);
0773                     dy_cl[i] = fabs(dy_cl[i] - ly);
0774                   }
0775                 }
0776               }  // valid clusterCollectionHandle
0777             }    // valid tracker
0778           }      // valid cpEstimator
0779           // distance of hit from closest cluster!
0780           float d_cl[2];
0781           d_cl[0] = d_cl[1] = -9999.;
0782           if (dx_cl[0] != -9999. && dy_cl[0] != -9999.)
0783             d_cl[0] = sqrt(dx_cl[0] * dx_cl[0] + dy_cl[0] * dy_cl[0]);
0784           if (dx_cl[1] != -9999. && dy_cl[1] != -9999.)
0785             d_cl[1] = sqrt(dx_cl[1] * dx_cl[1] + dy_cl[1] * dy_cl[1]);
0786           if (isHitMissing && (d_cl[0] < 0.05 || d_cl[1] < 0.05)) {
0787             isHitMissing = false;
0788             isHitValid = true;
0789           }
0790 
0791           if (debug_) {
0792             std::cout << "Ready to add hit in histogram:\n";
0793             // std::cout << "detid: "<<hit_detId<<std::endl;
0794             std::cout << "isHitValid: " << isHitValid << std::endl;
0795             std::cout << "isHitMissing: " << isHitMissing << std::endl;
0796             // std::cout << "passedEdgeCut: "<<passedFiducial<<std::endl;
0797           }
0798 
0799           if (nStripHits < 11)
0800             continue;  // Efficiency plots are filled with hits on tracks that
0801                        // have at least 11 Strip hits
0802 
0803           if (pxd != theSiPixelStructure.end() && isHitValid && passedFiducial)
0804             ++nvalid;
0805           if (pxd != theSiPixelStructure.end() && isHitMissing && passedFiducial)
0806             ++nmissing;
0807 
0808           if (pxd != theSiPixelStructure.end() && passedFiducial && (isHitValid || isHitMissing))
0809             (*pxd).second->fill(pTT, ltp, isHitValid, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
0810 
0811           //}//end if (persistent hit exists and is pixel hit)
0812 
0813         }  // end of else
0814 
0815       }  // end for (all traj measurements of pixeltrack)
0816     }    // end if (is pixeltrack)
0817     else if (debug_)
0818       std::cout << "no pixeltrack:\n";
0819 
0820   }  // end loop on map entries
0821 }
0822 
0823 DEFINE_FWK_MODULE(SiPixelHitEfficiencySource);  // define this as a plug-in