Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-11 04:32:37 

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 D1hits = 0;
0379     int D2hits = 0;
0380     std::vector<TrajectoryMeasurement> tmeasColl = traj_iterator->measurements();
0381     std::vector<TrajectoryMeasurement>::const_iterator tmeasIt;
0382     // loop on measurements to find out what kind of hits there are
0383     for (tmeasIt = tmeasColl.begin(); tmeasIt != tmeasColl.end(); tmeasIt++) {
0384       // if(! tmeasIt->updatedState().isValid()) continue; NOT NECESSARY (I
0385       // HOPE)
0386       TransientTrackingRecHit::ConstRecHitPointer testhit = tmeasIt->recHit();
0387       if (testhit->geographicalId().det() != DetId::Tracker)
0388         continue;
0389       uint testSubDetID = (testhit->geographicalId().subdetId());
0390       const DetId &hit_detId = testhit->geographicalId();
0391       int hit_layer = 0;
0392       int hit_ladder = 0;
0393       int hit_mod = 0;
0394       int hit_disk = 0;
0395 
0396       if (testSubDetID == PixelSubdetector::PixelBarrel) {
0397         isBpixtrack = true;
0398         hit_layer = PixelBarrelName(hit_detId, pTT, isUpgrade).layerName();
0399 
0400         hit_ladder = PXBDetId(hit_detId).ladder();
0401         hit_mod = PXBDetId(hit_detId).module();
0402 
0403         if (hit_layer == 1)
0404           L1hits++;
0405         if (hit_layer == 2)
0406           L2hits++;
0407         if (hit_layer == 3)
0408           L3hits++;
0409       }
0410       if (testSubDetID == PixelSubdetector::PixelEndcap) {
0411         isFpixtrack = true;
0412         hit_disk = PixelEndcapName(hit_detId, pTT, isUpgrade).diskName();
0413 
0414         if (hit_disk == 1)
0415           D1hits++;
0416         if (hit_disk == 2)
0417           D2hits++;
0418       }
0419       if (testSubDetID == StripSubdetector::TIB)
0420         nStripHits++;
0421       if (testSubDetID == StripSubdetector::TOB)
0422         nStripHits++;
0423       if (testSubDetID == StripSubdetector::TID)
0424         nStripHits++;
0425       if (testSubDetID == StripSubdetector::TEC)
0426         nStripHits++;
0427       // check if last valid hit is in Layer 2 or Disk 1
0428 
0429       bool lastValidL2 = false;
0430       if ((testSubDetID == PixelSubdetector::PixelBarrel && hit_layer == extrapolateFrom_) ||
0431           (testSubDetID == PixelSubdetector::PixelEndcap && hit_disk == 1)) {
0432         if (testhit->isValid()) {
0433           if (tmeasIt == tmeasColl.end() - 1) {
0434             lastValidL2 = true;
0435           } else {
0436             tmeasIt++;
0437             TransientTrackingRecHit::ConstRecHitPointer nextRecHit = tmeasIt->recHit();
0438             uint nextSubDetID = (nextRecHit->geographicalId().subdetId());
0439             int nextlayer = PixelBarrelName(nextRecHit->geographicalId()).layerName();
0440             if (nextSubDetID == PixelSubdetector::PixelBarrel && nextlayer == extrapolateTo_) {
0441               lastValidL2 = true;  //&& !nextRecHit->isValid()) lastValidL2=true;
0442             }
0443             tmeasIt--;
0444           }
0445         }
0446       }  // end check last valid layer
0447       if (lastValidL2) {
0448         std::vector<const BarrelDetLayer *> pxbLayers =
0449             measurementTrackerHandle->geometricSearchTracker()->pixelBarrelLayers();
0450         const DetLayer *pxb1 = pxbLayers[extrapolateTo_ - 1];
0451         const MeasurementEstimator *estimator = est.product();
0452         const LayerMeasurements *theLayerMeasurements =
0453             new LayerMeasurements(*measurementTrackerHandle, *measurementTrackerEventHandle);
0454         const TrajectoryStateOnSurface tsosPXB2 = tmeasIt->updatedState();
0455         expTrajMeasurements = theLayerMeasurements->measurements(*pxb1, tsosPXB2, *thePropagator, *estimator);
0456         delete theLayerMeasurements;
0457         if (!expTrajMeasurements.empty()) {
0458           for (uint p = 0; p < expTrajMeasurements.size(); p++) {
0459             TrajectoryMeasurement pxb1TM(expTrajMeasurements[p]);
0460             const auto &pxb1Hit = pxb1TM.recHit();
0461             // remove hits with rawID == 0
0462             if (pxb1Hit->geographicalId().rawId() == 0) {
0463               expTrajMeasurements.erase(expTrajMeasurements.begin() + p);
0464               continue;
0465             }
0466           }
0467         }
0468         //
0469       }
0470       // check if extrapolated hit to layer 1 one matches the original hit
0471       TrajectoryStateOnSurface chkPredTrajState =
0472           trajStateComb(tmeasIt->forwardPredictedState(), tmeasIt->backwardPredictedState());
0473       if (!chkPredTrajState.isValid())
0474         continue;
0475       float chkx = chkPredTrajState.globalPosition().x();
0476       float chky = chkPredTrajState.globalPosition().y();
0477       float chkz = chkPredTrajState.globalPosition().z();
0478       LocalPoint chklp = chkPredTrajState.localPosition();
0479       if (testSubDetID == PixelSubdetector::PixelBarrel && hit_layer == extrapolateTo_) {
0480         // Here we will drop the extrapolated hits if there is a hit and use
0481         // that hit
0482         vector<int> imatches;
0483         size_t imatch = 0;
0484         float glmatch = 9999.;
0485         for (size_t iexp = 0; iexp < expTrajMeasurements.size(); iexp++) {
0486           const DetId &exphit_detId = expTrajMeasurements[iexp].recHit()->geographicalId();
0487           int exphit_ladder = PXBDetId(exphit_detId).ladder();
0488           int exphit_mod = PXBDetId(exphit_detId).module();
0489           int dladder = abs(exphit_ladder - hit_ladder);
0490           if (dladder > 10)
0491             dladder = 20 - dladder;
0492           int dmodule = abs(exphit_mod - hit_mod);
0493           if (dladder != 0 || dmodule != 0) {
0494             continue;
0495           }
0496 
0497           TrajectoryStateOnSurface predTrajState = expTrajMeasurements[iexp].updatedState();
0498           float x = predTrajState.globalPosition().x();
0499           float y = predTrajState.globalPosition().y();
0500           float z = predTrajState.globalPosition().z();
0501           float dxyz = sqrt((chkx - x) * (chkx - x) + (chky - y) * (chky - y) + (chkz - z) * (chkz - z));
0502 
0503           if (dxyz <= glmatch) {
0504             glmatch = dxyz;
0505             imatch = iexp;
0506             imatches.push_back(int(imatch));
0507           }
0508 
0509         }  // found the propagated traj best matching the hit in data
0510 
0511         float lxmatch = 9999.0;
0512         float lymatch = 9999.0;
0513         if (!expTrajMeasurements.empty()) {
0514           if (glmatch < 9999.) {  // if there is any propagated trajectory for this hit
0515             const DetId &matchhit_detId = expTrajMeasurements[imatch].recHit()->geographicalId();
0516 
0517             int matchhit_ladder = PXBDetId(matchhit_detId).ladder();
0518             int dladder = abs(matchhit_ladder - hit_ladder);
0519             if (dladder > 10)
0520               dladder = 20 - dladder;
0521             LocalPoint lp = expTrajMeasurements[imatch].updatedState().localPosition();
0522             lxmatch = fabs(lp.x() - chklp.x());
0523             lymatch = fabs(lp.y() - chklp.y());
0524           }
0525           if (lxmatch < maxlxmatch_ && lymatch < maxlymatch_) {
0526             if (testhit->getType() != TrackingRecHit::missing || keepOriginalMissingHit_) {
0527               expTrajMeasurements.erase(expTrajMeasurements.begin() + imatch);
0528             }
0529           }
0530 
0531         }  // expected trajectory measurment not empty
0532       }
0533     }  // loop on trajectory measurments tmeasColl
0534 
0535     // if an extrapolated hit was found but not matched to an exisitng L1 hit
0536     // then push the hit back into the collection now keep the first one that is
0537     // left
0538     if (!expTrajMeasurements.empty()) {
0539       for (size_t f = 0; f < expTrajMeasurements.size(); f++) {
0540         TrajectoryMeasurement AddHit = expTrajMeasurements[f];
0541         if (AddHit.recHit()->getType() == TrackingRecHit::missing) {
0542           tmeasColl.push_back(AddHit);
0543           isBpixtrack = true;
0544         }
0545       }
0546     }
0547 
0548     if (isBpixtrack || isFpixtrack) {
0549       if (trackref->pt() < 0.6 || nStripHits < 8 || fabs(trackref->dxy(bestVtx->position())) > 0.05 ||
0550           fabs(trackref->dz(bestVtx->position())) > 0.5)
0551         continue;
0552 
0553       if (debug_) {
0554         std::cout << "isBpixtrack : " << isBpixtrack << std::endl;
0555         std::cout << "isFpixtrack : " << isFpixtrack << std::endl;
0556       }
0557       // std::cout<<"This tracks has so many hits:
0558       // "<<tmeasColl.size()<<std::endl;
0559       for (std::vector<TrajectoryMeasurement>::const_iterator tmeasIt = tmeasColl.begin(); tmeasIt != tmeasColl.end();
0560            tmeasIt++) {
0561         // if(! tmeasIt->updatedState().isValid()) continue;
0562         TrajectoryStateOnSurface tsos = tmeasIt->updatedState();
0563 
0564         TransientTrackingRecHit::ConstRecHitPointer hit = tmeasIt->recHit();
0565         if (hit->geographicalId().det() != DetId::Tracker)
0566           continue;
0567         else {
0568           //      //residual
0569           const DetId &hit_detId = hit->geographicalId();
0570           // uint IntRawDetID = (hit_detId.rawId());
0571           uint IntSubDetID = (hit_detId.subdetId());
0572 
0573           if (IntSubDetID == 0) {
0574             if (debug_)
0575               std::cout << "NO IntSubDetID\n";
0576             continue;
0577           }
0578           if (IntSubDetID != PixelSubdetector::PixelBarrel && IntSubDetID != PixelSubdetector::PixelEndcap)
0579             continue;
0580 
0581           int disk = 0;
0582           int layer = 0;
0583           int panel = 0;
0584           int module = 0;
0585           bool isHalfModule = false;
0586 
0587           const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit *>(hit->hit());
0588 
0589           if (IntSubDetID == PixelSubdetector::PixelBarrel) {  // it's a BPIX hit
0590             layer = PixelBarrelName(hit_detId, pTT, isUpgrade).layerName();
0591             isHalfModule = PixelBarrelName(hit_detId, pTT, isUpgrade).isHalfModule();
0592 
0593             if (hit->isValid()) {  // fill the cluster probability in barrel
0594               bool plus = true;
0595               if ((PixelBarrelName(hit_detId, pTT, isUpgrade).shell() == PixelBarrelName::Shell::mO) ||
0596                   (PixelBarrelName(hit_detId, pTT, isUpgrade).shell() == PixelBarrelName::Shell::mI))
0597                 plus = false;
0598               double clusterProbability = pixhit->clusterProbability(0);
0599               if (clusterProbability != 0)
0600                 fillClusterProbability(layer, 0, plus, log10(clusterProbability));
0601             }
0602 
0603           } else if (IntSubDetID == PixelSubdetector::PixelEndcap) {  // it's an FPIX hit
0604             disk = PixelEndcapName(hit_detId, pTT, isUpgrade).diskName();
0605             panel = PixelEndcapName(hit_detId, pTT, isUpgrade).pannelName();
0606             module = PixelEndcapName(hit_detId, pTT, isUpgrade).plaquetteName();
0607 
0608             if (hit->isValid()) {
0609               bool plus = true;
0610               if ((PixelEndcapName(hit_detId, pTT, isUpgrade).halfCylinder() == PixelEndcapName::HalfCylinder::mO) ||
0611                   (PixelEndcapName(hit_detId, pTT, isUpgrade).halfCylinder() == PixelEndcapName::HalfCylinder::mI))
0612                 plus = false;
0613               double clusterProbability = pixhit->clusterProbability(0);
0614               if (clusterProbability != 0)
0615                 fillClusterProbability(0, disk, plus, log10(clusterProbability));
0616             }
0617           }
0618 
0619           if (layer == 1) {
0620             if (fabs(trackref->dxy(bestVtx->position())) > 0.01 || fabs(trackref->dz(bestVtx->position())) > 0.1)
0621               continue;
0622             if (!(L2hits > 0 && L3hits > 0) && !(L2hits > 0 && D1hits > 0) && !(D1hits > 0 && D2hits > 0))
0623               continue;
0624           } else if (layer == 2) {
0625             if (fabs(trackref->dxy(bestVtx->position())) > 0.02 || fabs(trackref->dz(bestVtx->position())) > 0.1)
0626               continue;
0627             if (!(L1hits > 0 && L3hits > 0) && !(L1hits > 0 && D1hits > 0))
0628               continue;
0629           } else if (layer == 3) {
0630             if (fabs(trackref->dxy(bestVtx->position())) > 0.02 || fabs(trackref->dz(bestVtx->position())) > 0.1)
0631               continue;
0632             if (!(L1hits > 0 && L2hits > 0))
0633               continue;
0634           } else if (layer == 4) {
0635             if (fabs(trackref->dxy(bestVtx->position())) > 0.02 || fabs(trackref->dz(bestVtx->position())) > 0.1)
0636               continue;
0637           } else if (disk == 1) {
0638             if (fabs(trackref->dxy(bestVtx->position())) > 0.05 || fabs(trackref->dz(bestVtx->position())) > 0.5)
0639               continue;
0640             if (!(L1hits > 0 && D2hits > 0) && !(L2hits > 0 && D2hits > 0))
0641               continue;
0642           } else if (disk == 2) {
0643             if (fabs(trackref->dxy(bestVtx->position())) > 0.05 || fabs(trackref->dz(bestVtx->position())) > 0.5)
0644               continue;
0645             if (!(L1hits > 0 && D1hits > 0))
0646               continue;
0647           } else if (disk == 3) {
0648             if (fabs(trackref->dxy(bestVtx->position())) > 0.05 || fabs(trackref->dz(bestVtx->position())) > 0.5)
0649               continue;
0650           }
0651 
0652           // check whether hit is valid or missing using track algo flag
0653           bool isHitValid = hit->hit()->getType() == TrackingRecHit::valid;
0654           bool isHitMissing = hit->hit()->getType() == TrackingRecHit::missing;
0655 
0656           if (debug_)
0657             std::cout << "the hit is persistent\n";
0658 
0659           std::map<uint32_t, SiPixelHitEfficiencyModule *>::iterator pxd =
0660               theSiPixelStructure.find((*hit).geographicalId().rawId());
0661 
0662           // calculate alpha and beta from cluster position
0663           LocalTrajectoryParameters ltp = tsos.localParameters();
0664           // LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
0665 
0666           //*************** Edge cut ********************
0667           double lx = tsos.localPosition().x();
0668           double ly = tsos.localPosition().y();
0669 
0670           if (fabs(lx) > 0.55 || fabs(ly) > 3.0)
0671             continue;
0672 
0673           bool passedFiducial = true;
0674 
0675           // Module fiducials:
0676           if (IntSubDetID == PixelSubdetector::PixelBarrel && fabs(ly) >= 3.1)
0677             passedFiducial = false;
0678           if (IntSubDetID == PixelSubdetector::PixelEndcap &&
0679               !((panel == 1 &&
0680                  ((module == 1 && fabs(ly) < 0.7) ||
0681                   ((module == 2 && fabs(ly) < 1.1) && !(disk == -1 && ly > 0.8 && lx > 0.2) &&
0682                    !(disk == 1 && ly < -0.7 && lx > 0.2) && !(disk == 2 && ly < -0.8)) ||
0683                   ((module == 3 && fabs(ly) < 1.5) && !(disk == -2 && lx > 0.1 && ly > 1.0) &&
0684                    !(disk == 2 && lx > 0.1 && ly < -1.0)) ||
0685                   ((module == 4 && fabs(ly) < 1.9) && !(disk == -2 && ly > 1.5) && !(disk == 2 && ly < -1.5)))) ||
0686                 (panel == 2 &&
0687                  ((module == 1 && fabs(ly) < 0.7) ||
0688                   (module == 2 && fabs(ly) < 1.2 && !(disk > 0 && ly > 1.1) && !(disk < 0 && ly < -1.1)) ||
0689                   (module == 3 && fabs(ly) < 1.6 && !(disk > 0 && ly > 1.5) && !(disk < 0 && ly < -1.5))))))
0690             passedFiducial = false;
0691           if (IntSubDetID == PixelSubdetector::PixelEndcap &&
0692               ((panel == 1 && (module == 1 || (module >= 3 && abs(disk) == 1))) ||
0693                (panel == 2 && ((module == 1 && abs(disk) == 2) || (module == 3 && abs(disk) == 1)))))
0694             passedFiducial = false;
0695           // ROC fiducials:
0696           double ly_mod = fabs(ly);
0697           if (IntSubDetID == PixelSubdetector::PixelEndcap && (panel + module) % 2 == 1)
0698             ly_mod = ly_mod + 0.405;
0699           float d_rocedge = fabs(fmod(ly_mod, 0.81) - 0.405);
0700           if (d_rocedge <= 0.0625)
0701             passedFiducial = false;
0702           if (!((IntSubDetID == PixelSubdetector::PixelBarrel &&
0703                  ((!isHalfModule && fabs(lx) < 0.6) || (isHalfModule && lx > -0.3 && lx < 0.2))) ||
0704                 (IntSubDetID == PixelSubdetector::PixelEndcap &&
0705                  ((panel == 1 &&
0706                    ((module == 1 && fabs(lx) < 0.2) ||
0707                     (module == 2 && ((fabs(lx) < 0.55 && abs(disk) == 1) || (lx > -0.5 && lx < 0.2 && disk == -2) ||
0708                                      (lx > -0.5 && lx < 0.0 && disk == 2))) ||
0709                     (module == 3 && lx > -0.6 && lx < 0.5) || (module == 4 && lx > -0.3 && lx < 0.15))) ||
0710                   (panel == 2 && ((module == 1 && fabs(lx) < 0.6) ||
0711                                   (module == 2 && ((fabs(lx) < 0.55 && abs(disk) == 1) ||
0712                                                    (lx > -0.6 && lx < 0.5 && abs(disk) == 2))) ||
0713                                   (module == 3 && fabs(lx) < 0.5)))))))
0714             passedFiducial = false;
0715           if (((IntSubDetID == PixelSubdetector::PixelBarrel && !isHalfModule) ||
0716                (IntSubDetID == PixelSubdetector::PixelEndcap && !(panel == 1 && (module == 1 || module == 4)))) &&
0717               fabs(lx) < 0.06)
0718             passedFiducial = false;
0719 
0720           //*************** find closest clusters ********************
0721           float dx_cl[2];
0722           float dy_cl[2];
0723           dx_cl[0] = dx_cl[1] = dy_cl[0] = dy_cl[1] = -9999.;
0724           edm::ESHandle<PixelClusterParameterEstimator> cpEstimator =
0725               iSetup.getHandle(pixelClusterParameterEstimatorToken_);
0726           if (cpEstimator.isValid()) {
0727             const PixelClusterParameterEstimator &cpe(*cpEstimator);
0728             edm::ESHandle<TrackerGeometry> tracker = iSetup.getHandle(trackerGeomToken_);
0729             if (tracker.isValid()) {
0730               const TrackerGeometry *tkgeom = &(*tracker);
0731               edm::Handle<edmNew::DetSetVector<SiPixelCluster>> clusterCollectionHandle;
0732               iEvent.getByToken(clusterCollectionToken_, clusterCollectionHandle);
0733               if (clusterCollectionHandle.isValid()) {
0734                 const edmNew::DetSetVector<SiPixelCluster> &clusterCollection = *clusterCollectionHandle;
0735                 edmNew::DetSetVector<SiPixelCluster>::const_iterator itClusterSet = clusterCollection.begin();
0736                 float minD[2];
0737                 minD[0] = minD[1] = 10000.;
0738                 for (; itClusterSet != clusterCollection.end(); itClusterSet++) {
0739                   DetId detId(itClusterSet->id());
0740                   if (detId.rawId() != hit->geographicalId().rawId())
0741                     continue;
0742                   // unsigned int sdId=detId.subdetId();
0743                   const PixelGeomDetUnit *pixdet = (const PixelGeomDetUnit *)tkgeom->idToDetUnit(detId);
0744                   edmNew::DetSet<SiPixelCluster>::const_iterator itCluster = itClusterSet->begin();
0745                   for (; itCluster != itClusterSet->end(); ++itCluster) {
0746                     LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
0747                     PixelClusterParameterEstimator::ReturnType params = cpe.getParameters(*itCluster, *pixdet);
0748                     lp = std::get<0>(params);
0749                     float D = sqrt((lp.x() - lx) * (lp.x() - lx) + (lp.y() - ly) * (lp.y() - ly));
0750                     if (D < minD[0]) {
0751                       minD[1] = minD[0];
0752                       dx_cl[1] = dx_cl[0];
0753                       dy_cl[1] = dy_cl[0];
0754                       minD[0] = D;
0755                       dx_cl[0] = lp.x();
0756                       dy_cl[0] = lp.y();
0757                     } else if (D < minD[1]) {
0758                       minD[1] = D;
0759                       dx_cl[1] = lp.x();
0760                       dy_cl[1] = lp.y();
0761                     }
0762                   }  // loop on clusterSets
0763                 }  // loop on clusterCollection
0764                 for (size_t i = 0; i < 2; i++) {
0765                   if (minD[i] < 9999.) {
0766                     dx_cl[i] = fabs(dx_cl[i] - lx);
0767                     dy_cl[i] = fabs(dy_cl[i] - ly);
0768                   }
0769                 }
0770               }  // valid clusterCollectionHandle
0771             }  // valid tracker
0772           }  // valid cpEstimator
0773           // distance of hit from closest cluster!
0774           float d_cl[2];
0775           d_cl[0] = d_cl[1] = -9999.;
0776           if (dx_cl[0] != -9999. && dy_cl[0] != -9999.)
0777             d_cl[0] = sqrt(dx_cl[0] * dx_cl[0] + dy_cl[0] * dy_cl[0]);
0778           if (dx_cl[1] != -9999. && dy_cl[1] != -9999.)
0779             d_cl[1] = sqrt(dx_cl[1] * dx_cl[1] + dy_cl[1] * dy_cl[1]);
0780           if (isHitMissing && (d_cl[0] < 0.05 || d_cl[1] < 0.05)) {
0781             isHitMissing = false;
0782             isHitValid = true;
0783           }
0784 
0785           if (debug_) {
0786             std::cout << "Ready to add hit in histogram:\n";
0787             // std::cout << "detid: "<<hit_detId<<std::endl;
0788             std::cout << "isHitValid: " << isHitValid << std::endl;
0789             std::cout << "isHitMissing: " << isHitMissing << std::endl;
0790             // std::cout << "passedEdgeCut: "<<passedFiducial<<std::endl;
0791           }
0792 
0793           if (nStripHits < 11)
0794             continue;  // Efficiency plots are filled with hits on tracks that
0795                        // have at least 11 Strip hits
0796 
0797           if (pxd != theSiPixelStructure.end() && isHitValid && passedFiducial)
0798             ++nvalid;
0799           if (pxd != theSiPixelStructure.end() && isHitMissing && passedFiducial)
0800             ++nmissing;
0801 
0802           if (pxd != theSiPixelStructure.end() && passedFiducial && (isHitValid || isHitMissing))
0803             (*pxd).second->fill(pTT, ltp, isHitValid, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
0804 
0805           //}//end if (persistent hit exists and is pixel hit)
0806 
0807         }  // end of else
0808 
0809       }  // end for (all traj measurements of pixeltrack)
0810     }  // end if (is pixeltrack)
0811     else if (debug_)
0812       std::cout << "no pixeltrack:\n";
0813 
0814   }  // end loop on map entries
0815 }
0816 
0817 DEFINE_FWK_MODULE(SiPixelHitEfficiencySource);  // define this as a plug-in
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!