SiPixelPhase1Track/plugins/SiPixelPhase1TrackEfficiency.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     SiPixelPhase1TrackEfficiency
0004 // Class:       SiPixelPhase1TrackEfficiency
0005 //
0006
0007 // Original Author: Marcel Schneider
0008
0009 #include "DQM/SiPixelPhase1Common/interface/SiPixelPhase1Base.h"
0010 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0011 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0012 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014
0015 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0016 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0017 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0018 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0019 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0020 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0023 #include "DataFormats/VertexReco/interface/Vertex.h"
0024 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0025 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0026 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0027 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0028 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0029 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0030
0031 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0032 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0033 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0034 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0035 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0036 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0037 #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
0038
0039 ///commnet
0040
0041 namespace {
0042
0043   class SiPixelPhase1TrackEfficiency final : public SiPixelPhase1Base {
0044     enum { VALID, MISSING, INACTIVE, EFFICIENCY, VERTICES };
0045
0046   public:
0047     explicit SiPixelPhase1TrackEfficiency(const edm::ParameterSet& conf);
0048     void analyze(const edm::Event&, const edm::EventSetup&) override;
0049
0050   private:
0051     edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster>> clustersToken_;
0052     edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0053     edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0054     edm::EDGetTokenT<TrajTrackAssociationCollection> trajTrackCollectionToken_;
0055     edm::EDGetTokenT<MeasurementTrackerEvent> tracker_;  //new
0056     bool applyVertexCut_;
0057
0058     edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken_;
0059     edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken_;
0060     edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorToken_;
0061     edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> chi2MeasurementEstimatorBaseToken_;
0062     edm::ESGetToken<MeasurementTracker, CkfComponentsRecord> measurementTrackerToken_;
0063     edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> pixelClusterParameterEstimatorToken_;
0064
0065     const TrackerTopology* trackerTopology_;
0066     const Propagator* trackerPropagator_;
0067     const MeasurementEstimator* chi2MeasurementEstimator_;
0068   };
0069
0070   SiPixelPhase1TrackEfficiency::SiPixelPhase1TrackEfficiency(const edm::ParameterSet& iConfig)
0071       : SiPixelPhase1Base(iConfig)  //,
0072   {
0073     tracker_ = consumes<MeasurementTrackerEvent>(iConfig.getParameter<edm::InputTag>("tracker"));
0074     tracksToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0075     vtxToken_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryvertices"));
0076     applyVertexCut_ = iConfig.getUntrackedParameter<bool>("VertexCut", true);
0077     trajTrackCollectionToken_ =
0078         consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("trajectoryInput"));
0079     clustersToken_ = consumes<edmNew::DetSetVector<SiPixelCluster>>(iConfig.getParameter<edm::InputTag>("clusters"));
0080
0081     trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0082     trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0083     propagatorToken_ = esConsumes<Propagator, TrackingComponentsRecord>(edm::ESInputTag("", "PropagatorWithMaterial"));
0084     chi2MeasurementEstimatorBaseToken_ =
0085         esConsumes<Chi2MeasurementEstimatorBase, TrackingComponentsRecord>(edm::ESInputTag("", "Chi2"));
0086     measurementTrackerToken_ = esConsumes<MeasurementTracker, CkfComponentsRecord>();
0087     pixelClusterParameterEstimatorToken_ =
0088         esConsumes<PixelClusterParameterEstimator, TkPixelCPERecord>(edm::ESInputTag("", "PixelCPEGeneric"));
0089   }
0090
0091   void SiPixelPhase1TrackEfficiency::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0092     if (!checktrigger(iEvent, iSetup, DCS))
0093       return;
0094
0095     // get geometry
0096     edm::ESHandle<TrackerGeometry> tracker = iSetup.getHandle(trackerGeomToken_);
0097     assert(tracker.isValid());
0098
0099     // get primary vertex
0100     edm::Handle<reco::VertexCollection> vertices;
0101     iEvent.getByToken(vtxToken_, vertices);
0102
0103     // TrackerTopology for module informations
0104     edm::ESHandle<TrackerTopology> trackerTopologyHandle = iSetup.getHandle(trackerTopoToken_);
0105     trackerTopology_ = trackerTopologyHandle.product();
0106
0107     // Tracker propagator for propagating tracks to other layers
0108     edm::ESHandle<Propagator> propagatorHandle = iSetup.getHandle(propagatorToken_);
0109     std::unique_ptr<Propagator> propagatorUniquePtr(propagatorHandle.product()->clone());
0110     trackerPropagator_ = propagatorUniquePtr.get();
0111     const_cast<Propagator*>(trackerPropagator_)->setPropagationDirection(oppositeToMomentum);
0112
0113     // Measurement estimator
0114     edm::ESHandle<Chi2MeasurementEstimatorBase> chi2MeasurementEstimatorHandle =
0115         iSetup.getHandle(chi2MeasurementEstimatorBaseToken_);
0116     chi2MeasurementEstimator_ = chi2MeasurementEstimatorHandle.product();
0117
0118     //Tracker
0119     edm::Handle<MeasurementTrackerEvent> trackerMeas;
0120     iEvent.getByToken(tracker_, trackerMeas);
0121
0122     edm::ESHandle<MeasurementTracker> measurementTrackerHandle = iSetup.getHandle(measurementTrackerToken_);
0123
0124     //vertices
0125     if (!vertices.isValid())
0126       return;
0127     histo[VERTICES].fill(vertices->size(), DetId(0), &iEvent);
0128     if (applyVertexCut_ && vertices->empty())
0129       return;
0130
0131     // should be used for weird cuts
0132     //const auto primaryVertex = vertices->at(0);
0133
0134     // get the map
0135     edm::Handle<reco::TrackCollection> tracks;
0136     iEvent.getByToken(tracksToken_, tracks);
0137     if (!tracks.isValid())
0138       return;
0139
0140     //new
0141     edm::Handle<TrajTrackAssociationCollection> trajTrackCollectionHandle;
0142     iEvent.getByToken(trajTrackCollectionToken_, trajTrackCollectionHandle);
0143     if (!trajTrackCollectionHandle.isValid())
0144       return;
0145
0146     //Access Pixel Clusters
0147     edm::Handle<edmNew::DetSetVector<SiPixelCluster>> siPixelClusters;
0148     iEvent.getByToken(clustersToken_, siPixelClusters);
0149     if (!siPixelClusters.isValid())
0150       return;
0151     //
0152
0153     edm::ESHandle<PixelClusterParameterEstimator> cpEstimator = iSetup.getHandle(pixelClusterParameterEstimatorToken_);
0154     if (!cpEstimator.isValid())
0155       return;
0156
0157     const PixelClusterParameterEstimator& cpe(*cpEstimator);
0158     const TrackerGeometry* tkgeom = &(*tracker);
0159
0160     //////////////////////////////////////////////////////////////////////////////////////////
0161
0162     // Hp cut
0163     static constexpr int TRACK_QUALITY_HIGH_PURITY_BIT = 2;
0164     static constexpr int TRACK_QUALITY_HIGH_PURITY_MASK = 1 << TRACK_QUALITY_HIGH_PURITY_BIT;
0165
0166     // Pt cut
0167     static constexpr float TRACK_PT_CUT_VAL = 1.0f;
0168
0169     // Nstrip cut
0170     static constexpr int TRACK_NSTRIP_CUT_VAL = 10;
0171
0172     //D0
0173     static constexpr std::array<float, 4> TRACK_D0_CUT_BARREL_VAL = {{0.01f, 0.02f, 0.02f, 0.02f}};
0174     static constexpr float TRACK_D0_CUT_FORWARD_VAL = 0.05f;
0175
0176     //Dz
0177     static constexpr float TRACK_DZ_CUT_BARREL_VAL = 0.01f;
0178     static constexpr float TRACK_DZ_CUT_FORWARD_VAL = 0.5f;
0179
0180     TrajectoryStateOnSurface tsosPXB2;
0181     bool valid_layerFrom = false;
0182
0183     const GeometricSearchTracker* gst_ = trackerMeas->geometricSearchTracker();
0184     const auto* pxbLayer1_ = gst_->pixelBarrelLayers().front();
0185     const LayerMeasurements* theLayerMeasurements_ = new LayerMeasurements(*measurementTrackerHandle, *trackerMeas);
0186
0187     std::vector<TrajectoryMeasurement> expTrajMeasurements;
0188     std::vector<std::pair<int, bool[3]>> eff_pxb1_vector;
0189
0190     for (const auto& pair : *trajTrackCollectionHandle) {
0191       const edm::Ref<std::vector<Trajectory>> traj = pair.key;
0192       const reco::TrackRef track = pair.val;
0193
0194       expTrajMeasurements.clear();
0195       eff_pxb1_vector.clear();
0196       //this cut is needed to be consisten with residuals calculation
0197       if (applyVertexCut_ &&
0198           (track->pt() < 0.75 || std::abs(track->dxy(vertices->at(0).position())) > 5 * track->dxyError()))
0199         continue;
0200
0201       bool isBpixtrack = false;
0202       bool isFpixtrack = false;
0203       int nStripHits = 0;
0204       int nBpixL1Hits = 0;
0205       int nBpixL2Hits = 0;
0206       int nBpixL3Hits = 0;
0207       int nBpixL4Hits = 0;
0208       int nFpixD1Hits = 0;
0209       int nFpixD2Hits = 0;
0210       int nFpixD3Hits = 0;
0211       bool passcuts = true;
0212
0213       // first, look at the full track to see whether it is good
0214       // auto const & trajParams = track.extra()->trajParams();
0215
0216       auto hb = track->recHitsBegin();
0217       for (unsigned int h = 0; h < track->recHitsSize(); h++) {
0218         auto hit = *(hb + h);
0219         if (!hit->isValid())
0220           continue;
0221
0222         DetId id = hit->geographicalId();
0223         uint32_t subdetid = (id.subdetId());
0224
0225         //Check the location of valid hit
0226         if (subdetid == PixelSubdetector::PixelBarrel && hit->isValid()) {
0227           isBpixtrack = true;
0228           if (trackerTopology_->pxbLayer(id) == 1)
0229             nBpixL1Hits++;
0230           else if (trackerTopology_->pxbLayer(id) == 2)
0231             nBpixL2Hits++;
0232           else if (trackerTopology_->pxbLayer(id) == 3)
0233             nBpixL3Hits++;
0234           else if (trackerTopology_->pxbLayer(id) == 4)
0235             nBpixL4Hits++;
0236         } else if (subdetid == PixelSubdetector::PixelEndcap && hit->isValid()) {
0237           isFpixtrack = true;
0238           if (trackerTopology_->pxfDisk(id) == 1)
0239             nFpixD1Hits++;
0240           else if (trackerTopology_->pxfDisk(id) == 2)
0241             nFpixD2Hits++;
0242           else if (trackerTopology_->pxfDisk(id) == 3)
0243             nFpixD3Hits++;
0244         }
0245
0246         // count strip hits
0247         if (subdetid == StripSubdetector::TIB || subdetid == StripSubdetector::TOB ||
0248             subdetid == StripSubdetector::TID || subdetid == StripSubdetector::TEC)
0249           nStripHits++;
0250       }
0251
0252       if (!isBpixtrack && !isFpixtrack)
0253         continue;
0254
0255       // Hp cut
0256       if (!((track->qualityMask() & TRACK_QUALITY_HIGH_PURITY_MASK) >> TRACK_QUALITY_HIGH_PURITY_BIT))
0257         passcuts = false;
0258
0259       // Pt cut
0260       if (!(TRACK_PT_CUT_VAL < track->pt()))
0261         passcuts = false;
0262
0263       // Nstrip cut
0264       if (!(TRACK_NSTRIP_CUT_VAL < nStripHits))
0265         passcuts = false;
0266
0267       // then, look at each hit
0268       for (unsigned int h = 0; h < track->recHitsSize(); h++) {
0269         bool passcuts_hit = true;
0270         auto hit = *(hb + h);
0271
0272         DetId id = hit->geographicalId();
0273         uint32_t subdetid = (id.subdetId());
0274         if (subdetid != PixelSubdetector::PixelBarrel && subdetid != PixelSubdetector::PixelEndcap)
0275           continue;
0276
0277         bool isHitValid = hit->getType() == TrackingRecHit::valid;
0278         bool isHitMissing = hit->getType() == TrackingRecHit::missing;
0279         bool isHitInactive = hit->getType() == TrackingRecHit::inactive;
0280
0281         //D0
0282         if (subdetid == PixelSubdetector::PixelBarrel) {
0283           if (!((std::abs(track->dxy(vertices->at(0).position())) * -1.0) <
0284                 TRACK_D0_CUT_BARREL_VAL[trackerTopology_->pxbLayer(id) - 1]))
0285             passcuts_hit = false;
0286         } else if (subdetid == PixelSubdetector::PixelEndcap) {
0287           if (!((std::abs(track->dxy(vertices->at(0).position())) * -1.0) < TRACK_D0_CUT_FORWARD_VAL))
0288             passcuts_hit = false;
0289         }
0290
0291         //Dz
0292         if (subdetid == PixelSubdetector::PixelBarrel) {
0293           if (!(std::abs(track->dz(vertices->at(0).position())) < TRACK_DZ_CUT_BARREL_VAL))
0294             passcuts_hit = false;
0295         } else if (subdetid == PixelSubdetector::PixelEndcap) {
0296           if (!(std::abs(track->dz(vertices->at(0).position())) < TRACK_DZ_CUT_FORWARD_VAL))
0297             passcuts_hit = false;
0298         }
0299
0300         // Pixhit cut
0301         if (subdetid == PixelSubdetector::PixelBarrel) {
0302           if (trackerTopology_->pxbLayer(id) == 1) {
0303             if (!((nBpixL2Hits > 0 && nBpixL3Hits > 0 && nBpixL4Hits > 0) ||
0304                   (nBpixL2Hits > 0 && nBpixL3Hits > 0 && nFpixD1Hits > 0) ||
0305                   (nBpixL2Hits > 0 && nFpixD1Hits > 0 && nFpixD2Hits > 0) ||
0306                   (nFpixD1Hits > 0 && nFpixD2Hits > 0 && nFpixD3Hits > 0)))
0307               passcuts_hit = false;
0308           } else if (trackerTopology_->pxbLayer(id) == 2) {
0309             if (!((nBpixL1Hits > 0 && nBpixL3Hits > 0 && nBpixL4Hits > 0) ||
0310                   (nBpixL1Hits > 0 && nBpixL3Hits > 0 && nFpixD1Hits > 0) ||
0311                   (nBpixL1Hits > 0 && nFpixD1Hits > 0 && nFpixD2Hits > 0)))
0312               passcuts_hit = false;
0313           } else if (trackerTopology_->pxbLayer(id) == 3) {
0314             if (!((nBpixL1Hits > 0 && nBpixL2Hits > 0 && nBpixL4Hits > 0) ||
0315                   (nBpixL1Hits > 0 && nBpixL2Hits > 0 && nFpixD1Hits > 0)))
0316               passcuts_hit = false;
0317           } else if (trackerTopology_->pxbLayer(id) == 4)
0318             if (!((nBpixL1Hits > 0 && nBpixL2Hits > 0 && nBpixL3Hits > 0)))
0319               passcuts_hit = false;
0320         } else if (subdetid == PixelSubdetector::PixelEndcap) {
0321           if (trackerTopology_->pxfDisk(id) == 1) {
0322             if (!((nBpixL1Hits > 0 && nBpixL2Hits > 0 && nBpixL3Hits > 0) ||
0323                   (nBpixL1Hits > 0 && nBpixL2Hits > 0 && nFpixD2Hits > 0) ||
0324                   (nBpixL1Hits > 0 && nFpixD2Hits > 0 && nFpixD3Hits > 0)))
0325               passcuts_hit = false;
0326           } else if (trackerTopology_->pxfDisk(id) == 2) {
0327             if (!((nBpixL1Hits > 0 && nBpixL2Hits > 0 && nFpixD1Hits > 0) ||
0328                   (nBpixL1Hits > 0 && nFpixD1Hits > 0 && nFpixD3Hits > 0)))
0329               passcuts_hit = false;
0330           } else if (trackerTopology_->pxfDisk(id) == 3) {
0331             if (!((nBpixL1Hits > 0 && nFpixD1Hits > 0 && nFpixD2Hits > 0)))
0332               passcuts_hit = false;
0333           }
0334         }
0335         /*
0336       //Fiducial Cut - will work on it later
0337       const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(hit);
0338       const PixelGeomDetUnit* geomdetunit = dynamic_cast<const PixelGeomDetUnit*> ( tracker->idToDet(id) );
0339       const PixelTopology& topol = geomdetunit->specificTopology();
0340
0341       LocalPoint lp;
0342       if (pixhit) {
0343     lp = pixhit->localPosition();
0344       }
0345
0346       MeasurementPoint mp = topol.measurementPosition(lp);
0347       int row = (int) mp.x() % 80;
0348       int col = (int) mp.y() % 52;
0349
0350       int centerrow = 40;
0351       int centercol = 26;
0352
0353       if (!((col < (centercol + 10)) && (col > (centercol - 10)) && (row < (centerrow + 10)) && (row > (centerrow -10 )))) passcuts_hit = false;
0354       */
0355
0356         if (passcuts_hit && passcuts) {
0357           if (!(subdetid == PixelSubdetector::PixelBarrel && trackerTopology_->pxbLayer(id) == 1)) {
0358             if (isHitValid) {
0359               histo[VALID].fill(id, &iEvent);
0360               histo[EFFICIENCY].fill(1, id, &iEvent);
0361             }
0362             if (isHitMissing) {
0363               histo[MISSING].fill(id, &iEvent);
0364               histo[EFFICIENCY].fill(0, id, &iEvent);
0365             }
0366             if (isHitInactive) {
0367               histo[INACTIVE].fill(id, &iEvent);
0368             }
0369           }
0370         }
0371       }
0372
0373       ///////////////////////////////////////////////layer 1 specific here/////////////////////////////////////////////////////////////////////
0374       valid_layerFrom = false;
0375
0376       //propagation only from PXB2 and PXD1, more cuts later
0377       for (const auto& tm : traj->measurements()) {
0378         if (tm.recHit().get() && tm.recHitR().isValid()) {
0379           DetId where = tm.recHitR().geographicalId();
0380           int source_det = where.subdetId();
0381
0382           if (source_det == PixelSubdetector::SubDetector::PixelBarrel) {
0383             int source_layer = trackerTopology_->pxbLayer(where);
0384             if (source_layer == 2) {
0385               if (tm.updatedState().isValid()) {
0386                 tsosPXB2 = tm.updatedState();
0387                 valid_layerFrom = true;
0388               }
0389             }
0390           }
0391
0392           if (source_det == PixelSubdetector::SubDetector::PixelEndcap) {
0393             int source_layer = trackerTopology_->pxfDisk(where);
0394             if (source_layer == 1) {
0395               if (tm.updatedState().isValid()) {
0396                 tsosPXB2 = tm.updatedState();
0397                 valid_layerFrom = true;
0398               }
0399             }
0400           }
0401         }
0402       }  //uodated tsosPXB2 here
0403
0404       if (!valid_layerFrom)
0405         continue;
0406       if (!tsosPXB2.isValid())
0407         continue;
0408
0409       //propagation A: Calculate the efficiency by the distance to the closest cluster
0410       expTrajMeasurements =
0411           theLayerMeasurements_->measurements(*pxbLayer1_, tsosPXB2, *trackerPropagator_, *chi2MeasurementEstimator_);
0412       auto compDets = pxbLayer1_->compatibleDets(tsosPXB2, *trackerPropagator_, *chi2MeasurementEstimator_);
0413       std::pair<int, bool[3]> eff_map;
0414
0415       //Fiducial Cut, only calculate the efficiency of the central pixels
0416       for (uint p = 0; p < expTrajMeasurements.size(); p++) {
0417         bool valid = false;
0418         bool missing = false;
0419         bool passcuts_hit = true;
0420         TrajectoryMeasurement pxb1TM(expTrajMeasurements[p]);
0421         const auto& pxb1Hit = pxb1TM.recHit();
0422         bool inactive = (pxb1Hit->getType() == TrackingRecHit::inactive);
0423         int detidHit = pxb1Hit->geographicalId();
0424         if (detidHit == 0)
0425           continue;
0426
0427         const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(pxb1Hit->hit());
0428         const PixelGeomDetUnit* geomdetunit = dynamic_cast<const PixelGeomDetUnit*>(tracker->idToDet(detidHit));
0429         const PixelTopology& topol = geomdetunit->specificTopology();
0430
0431         if (!pixhit)
0432           continue;
0433
0434         LocalPoint lp = pixhit->localPosition();
0435         MeasurementPoint mp = topol.measurementPosition(lp);
0436         const int nRows = topol.rowsperroc();
0437         const int nColumns = topol.colsperroc();
0438         int row = (int)mp.x() % nRows;
0439         int col = (int)mp.y() % nColumns;
0440
0441         int centerrow = nRows / 2;
0442         int centercol = nColumns / 2;
0443
0444         if (!((col < (centercol + 10)) && (col > (centercol - 10)) && (row < (centerrow + 10)) &&
0445               (row > (centerrow - 10))))
0446           passcuts_hit = false;
0447
0448         //Access the distance to the closest cluster
0449         for (const auto& detAndState : compDets) {
0450           const auto& pXb1_lpos = detAndState.second.localPosition();
0451           if (pxb1Hit->geographicalId().rawId() != detAndState.first->geographicalId().rawId())
0452             continue;
0453           int detid = detAndState.first->geographicalId().rawId();
0454
0455           for (edmNew::DetSetVector<SiPixelCluster>::const_iterator iter_cl = siPixelClusters->begin();
0456                iter_cl != siPixelClusters->end();
0457                iter_cl++) {
0458             DetId detId(iter_cl->id());
0459             float minD[2], minDist = 10000;
0460             minD[0] = minD[1] = 10000.;
0461             if (detId.rawId() != detAndState.first->geographicalId().rawId())
0462               continue;
0463
0464             const PixelGeomDetUnit* pixdet = (const PixelGeomDetUnit*)tkgeom->idToDetUnit(detId);
0465             edmNew::DetSet<SiPixelCluster>::const_iterator itCluster = iter_cl->begin();
0466             if (passcuts_hit) {
0467               for (; itCluster != iter_cl->end(); ++itCluster) {
0468                 LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
0469                 PixelClusterParameterEstimator::ReturnType params = cpe.getParameters(*itCluster, *pixdet);
0470                 lp = std::get<0>(params);
0471
0472                 float Xdist = abs(lp.x() - pXb1_lpos.x());
0473                 float Ydist = abs(lp.y() - pXb1_lpos.y());
0474                 float dist = sqrt(Xdist * Xdist + Ydist * Ydist);
0475                 if (dist < minDist) {
0476                   minDist = dist;
0477                   minD[0] = Xdist;
0478                   minD[1] = Ydist;
0479                 }
0480               }
0481
0482               if ((minD[0] < 0.02) && (minD[1] < 0.02)) {
0483                 valid = true;
0484                 missing = false;
0485                 inactive = false;
0486               } else if (inactive) {
0487                 valid = false;
0488                 missing = false;
0489               } else {
0490                 missing = true;
0491                 valid = false;
0492               }
0493             }
0494           }
0495
0496           //cuts: exactly the same as for other hits but assuming PXB1
0497
0498           //D0
0499           if (!((std::abs(track->dxy(vertices->at(0).position())) * -1.0) <
0500                 TRACK_D0_CUT_BARREL_VAL[trackerTopology_->pxbLayer(detid) - 1]))
0501             passcuts_hit = false;
0502           //Dz
0503           if (!(std::abs(track->dz(vertices->at(0).position())) < TRACK_DZ_CUT_BARREL_VAL))
0504             passcuts_hit = false;
0505           // Pixhit cut
0506           if (!((nBpixL2Hits > 0 && nBpixL3Hits > 0 && nBpixL4Hits > 0) ||
0507                 (nBpixL2Hits > 0 && nBpixL3Hits > 0 && nFpixD1Hits > 0) ||
0508                 (nBpixL2Hits > 0 && nFpixD1Hits > 0 && nFpixD2Hits > 0) ||
0509                 (nFpixD1Hits > 0 && nFpixD2Hits > 0 && nFpixD3Hits > 0)))
0510             passcuts_hit = false;
0511           bool found_det = false;
0512
0513           if (passcuts && passcuts_hit) {
0514             for (unsigned int i_eff = 0; i_eff < eff_pxb1_vector.size(); i_eff++) {
0515               //in case found hit in the same det, take only the valid hit
0516               if (eff_pxb1_vector[i_eff].first == detid) {
0517                 found_det = true;
0518                 if (eff_pxb1_vector[i_eff].second[0] == false && valid == true) {
0519                   eff_pxb1_vector[i_eff].second[0] = valid;
0520                   eff_pxb1_vector[i_eff].second[1] = missing;
0521                   eff_pxb1_vector[i_eff].second[2] = inactive;
0522                 }
0523               }
0524             }
0525             if (!found_det) {
0526               eff_map.first = detid;
0527               eff_map.second[0] = valid;
0528               eff_map.second[1] = missing;
0529               eff_map.second[2] = inactive;
0530               eff_pxb1_vector.push_back(eff_map);
0531             }
0532           }
0533         }
0534       }
0535
0536       if (eff_pxb1_vector.size() == 1) {
0537         //eff map is filled -> decide what to do for double hits, ie eff_pxb1_vector.size>1 ... if 1 just use MISSING and VALID as usual
0538
0539         if (eff_pxb1_vector[0].second[0]) {
0540           histo[VALID].fill(eff_pxb1_vector[0].first, &iEvent);
0541           histo[EFFICIENCY].fill(1, eff_pxb1_vector[0].first, &iEvent);
0542         }
0543         if (eff_pxb1_vector[0].second[1]) {
0544           histo[MISSING].fill(eff_pxb1_vector[0].first, &iEvent);
0545           histo[EFFICIENCY].fill(0, eff_pxb1_vector[0].first, &iEvent);
0546         }
0547
0548         if (eff_pxb1_vector[0].second[2]) {
0549           histo[INACTIVE].fill(eff_pxb1_vector[0].first, &iEvent);
0550         }
0551       }
0552
0553     }  //trajTrackHandle
0554
0555     histo[VALID].executePerEventHarvesting(&iEvent);
0556     histo[MISSING].executePerEventHarvesting(&iEvent);
0557     histo[INACTIVE].executePerEventHarvesting(&iEvent);
0558   }
0559
0560 }  // namespace
0561
0562 DEFINE_FWK_MODULE(SiPixelPhase1TrackEfficiency);