SiPixelPhase1Track/plugins/SiPixelPhase1TrackClusters.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     SiPixelPhase1TrackClusters
0004 // Class  :     SiPixelPhase1TrackClusters
0005 //
0006
0007 // Original Author: Marcel Schneider
0008
0009 #include "DQM/SiPixelPhase1Common/interface/SiPixelPhase1Base.h"
0010 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0011
0012 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0013 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016
0017 #include "DataFormats/SiPixelCluster/interface/SiPixelClusterShapeCache.h"
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0022 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0023 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0024 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0025 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0026 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0027 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0028
0029 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0030
0031 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0032 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h"
0033
0034 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h"
0035 #include "CalibTracker/Records/interface/SiPixelTemplateDBObjectESProducerRcd.h"
0036
0037 namespace {
0038
0039   class SiPixelPhase1TrackClusters final : public SiPixelPhase1Base {
0040     enum {
0041       ON_TRACK_CHARGE,
0042       ON_TRACK_CORRECTEDCHARGE,
0043       TEMPLATE_CORRECTION,
0044       ON_TRACK_BIGPIXELCHARGE,
0045       ON_TRACK_NOTBIGPIXELCHARGE,
0046       ON_TRACK_SIZE,
0047       ON_TRACK_SHAPE,
0048       ON_TRACK_NCLUSTERS,
0049       ON_TRACK_POSITIONB,
0050       ON_TRACK_POSITIONF,
0051       DIGIS_HITMAP_ON_TRACK,
0052       ON_TRACK_NDIGIS,
0053
0054       NTRACKS,
0055       NTRACKS_INVOLUME,
0056
0057       SIZE_VS_ETA_ON_TRACK_OUTER,
0058       SIZE_VS_ETA_ON_TRACK_INNER,
0059       ON_TRACK_CHARGE_OUTER,
0060       ON_TRACK_CHARGE_INNER,
0061       ON_TRACK_CORRECTEDCHARGE_OUTER,
0062       ON_TRACK_CORRECTEDCHARGE_INNER,
0063       TEMPLATE_CORRECTION_OUTER,
0064       TEMPLATE_CORRECTION_INNER,
0065
0066       ON_TRACK_SHAPE_OUTER,
0067       ON_TRACK_SHAPE_INNER,
0068
0069       ON_TRACK_SIZE_X_OUTER,
0070       ON_TRACK_SIZE_X_INNER,
0071       ON_TRACK_SIZE_X_F,
0072       ON_TRACK_SIZE_Y_OUTER,
0073       ON_TRACK_SIZE_Y_INNER,
0074       ON_TRACK_SIZE_Y_F,
0075
0076       ON_TRACK_SIZE_XY_OUTER,
0077       ON_TRACK_SIZE_XY_INNER,
0078       ON_TRACK_SIZE_XY_F,
0079       CHARGE_VS_SIZE_ON_TRACK,
0080
0081       ENUM_SIZE
0082     };
0083
0084   public:
0085     explicit SiPixelPhase1TrackClusters(const edm::ParameterSet& conf);
0086     void dqmBeginRun(const edm::Run&, const edm::EventSetup&) override;
0087     void analyze(const edm::Event&, const edm::EventSetup&) override;
0088
0089   private:
0090     const bool applyVertexCut_;
0091     const SiPixelTemplateDBObject* templateDBobject_ = nullptr;
0092     const std::vector<SiPixelTemplateStore>* thePixelTemp_ = nullptr;
0093     const TrackerTopology* tkTpl = nullptr;
0094
0095     edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0096     edm::EDGetTokenT<reco::VertexCollection> offlinePrimaryVerticesToken_;
0097     edm::EDGetTokenT<SiPixelClusterShapeCache> pixelClusterShapeCacheToken_;
0098
0099     edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken_;
0100     edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken_;
0101     edm::ESGetToken<ClusterShapeHitFilter, CkfComponentsRecord> clusterShapeHitFilterToken_;
0102
0103     edm::ESGetToken<SiPixelTemplateDBObject, SiPixelTemplateDBObjectESProducerRcd> templateDBobjectToken_;
0104     edm::ESGetToken<std::vector<SiPixelTemplateStore>, SiPixelTemplateDBObjectESProducerRcd> templateStoreToken_;
0105   };
0106
0107   SiPixelPhase1TrackClusters::SiPixelPhase1TrackClusters(const edm::ParameterSet& iConfig)
0108       : SiPixelPhase1Base(iConfig), applyVertexCut_(iConfig.getUntrackedParameter<bool>("VertexCut", true)) {
0109     tracksToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0110
0111     offlinePrimaryVerticesToken_ =
0112         applyVertexCut_ ? consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))
0113                         : edm::EDGetTokenT<reco::VertexCollection>();
0114
0115     pixelClusterShapeCacheToken_ =
0116         consumes<SiPixelClusterShapeCache>(iConfig.getParameter<edm::InputTag>("clusterShapeCache"));
0117
0118     trackerTopoToken_ = esConsumes<edm::Transition::BeginRun>();
0119     trackerGeomToken_ = esConsumes<edm::Transition::BeginRun>();
0120     clusterShapeHitFilterToken_ =
0121         esConsumes<ClusterShapeHitFilter, CkfComponentsRecord>(edm::ESInputTag("", "ClusterShapeHitFilter"));
0122     templateStoreToken_ = esConsumes<edm::Transition::BeginRun>();
0123     templateDBobjectToken_ = esConsumes<edm::Transition::BeginRun>();
0124   }
0125
0126   void SiPixelPhase1TrackClusters::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0127     //get topology
0128     tkTpl = &iSetup.getData(trackerTopoToken_);
0129
0130     // Initialize 1D templates
0131     templateDBobject_ = &iSetup.getData(templateDBobjectToken_);
0132     thePixelTemp_ = &iSetup.getData(templateStoreToken_);
0133   }
0134
0135   void SiPixelPhase1TrackClusters::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0136     if (!checktrigger(iEvent, iSetup, DCS))
0137       return;
0138
0139     if (histo.size() != ENUM_SIZE) {
0140       edm::LogError("SiPixelPhase1TrackClusters")
0141           << "incompatible configuration " << histo.size() << "!=" << ENUM_SIZE << std::endl;
0142       return;
0143     }
0144
0145     edm::ESHandle<ClusterShapeHitFilter> shapeFilterH = iSetup.getHandle(clusterShapeHitFilterToken_);
0146     auto const& shapeFilter = *shapeFilterH;
0147
0148     edm::Handle<reco::VertexCollection> vertices;
0149     if (applyVertexCut_) {
0150       iEvent.getByToken(offlinePrimaryVerticesToken_, vertices);
0151       if (!vertices.isValid() || vertices->empty())
0152         return;
0153     }
0154
0155     //get the map
0156     edm::Handle<reco::TrackCollection> tracks;
0157     iEvent.getByToken(tracksToken_, tracks);
0158
0159     if (!tracks.isValid()) {
0160       edm::LogWarning("SiPixelPhase1TrackClusters") << "track collection is not valid";
0161       return;
0162     }
0163
0164     SiPixelTemplate templ(*thePixelTemp_);
0165
0166     edm::Handle<SiPixelClusterShapeCache> pixelClusterShapeCacheH;
0167     iEvent.getByToken(pixelClusterShapeCacheToken_, pixelClusterShapeCacheH);
0168     if (!pixelClusterShapeCacheH.isValid()) {
0169       edm::LogWarning("SiPixelPhase1TrackClusters") << "PixelClusterShapeCache collection is not valid";
0170       return;
0171     }
0172     auto const& pixelClusterShapeCache = *pixelClusterShapeCacheH;
0173
0174     for (auto const& track : *tracks) {
0175       if (applyVertexCut_ &&
0176           (track.pt() < 0.75 || std::abs(track.dxy((*vertices)[0].position())) > 5 * track.dxyError()))
0177         continue;
0178
0179       bool isBpixtrack = false, isFpixtrack = false, crossesPixVol = false;
0180
0181       // find out whether track crosses pixel fiducial volume (for cosmic tracks)
0182       auto d0 = track.d0(), dz = track.dz();
0183       if (std::abs(d0) < 16 && std::abs(dz) < 50)
0184         crossesPixVol = true;
0185
0186       auto etatk = track.eta();
0187
0188       auto const& trajParams = track.extra()->trajParams();
0189       assert(trajParams.size() == track.recHitsSize());
0190       auto hb = track.recHitsBegin();
0191
0192       for (unsigned int h = 0; h < track.recHitsSize(); h++) {
0193         auto hit = *(hb + h);
0194         if (!hit->isValid())
0195           continue;
0196         auto id = hit->geographicalId();
0197
0198         // check that we are in the pixel
0199         auto subdetid = (id.subdetId());
0200         if (subdetid == PixelSubdetector::PixelBarrel)
0201           isBpixtrack = true;
0202         if (subdetid == PixelSubdetector::PixelEndcap)
0203           isFpixtrack = true;
0204         if (subdetid != PixelSubdetector::PixelBarrel && subdetid != PixelSubdetector::PixelEndcap)
0205           continue;
0206         bool iAmBarrel = subdetid == PixelSubdetector::PixelBarrel;
0207
0208         // PXB_L4 IS IN THE OTHER WAY
0209         // CAN BE XORed BUT LETS KEEP THINGS SIMPLE
0210         bool iAmOuter = ((tkTpl->pxbLadder(id) % 2 == 1) && tkTpl->pxbLayer(id) != 4) ||
0211                         ((tkTpl->pxbLadder(id) % 2 != 1) && tkTpl->pxbLayer(id) == 4);
0212
0213         auto pixhit = dynamic_cast<const SiPixelRecHit*>(hit->hit());
0214         if (!pixhit)
0215           continue;
0216
0217         auto geomdetunit = dynamic_cast<const PixelGeomDetUnit*>(pixhit->detUnit());
0218         auto const& topol = geomdetunit->specificTopology();
0219
0220         // get the cluster
0221         auto clustp = pixhit->cluster();
0222         if (clustp.isNull())
0223           continue;
0224         auto const& cluster = *clustp;
0225         const std::vector<SiPixelCluster::Pixel> pixelsVec = cluster.pixels();
0226         for (unsigned int i = 0; i < pixelsVec.size(); ++i) {
0227           float pixx = pixelsVec[i].x;  // index as float=iteger, row index
0228           float pixy = pixelsVec[i].y;  // same, col index
0229
0230           bool bigInX = topol.isItBigPixelInX(int(pixx));
0231           bool bigInY = topol.isItBigPixelInY(int(pixy));
0232           float pixel_charge = pixelsVec[i].adc;
0233
0234           if (bigInX == true || bigInY == true) {
0235             histo[ON_TRACK_BIGPIXELCHARGE].fill(pixel_charge, id, &iEvent);
0236           } else {
0237             histo[ON_TRACK_NOTBIGPIXELCHARGE].fill(pixel_charge, id, &iEvent);
0238           }
0239         }  // End loop over pixels
0240         auto const& ltp = trajParams[h];
0241
0242         auto localDir = ltp.momentum() / ltp.momentum().mag();
0243
0244         // correct charge for track impact angle
0245         auto charge = cluster.charge() * ltp.absdz();
0246         //Correct charge with Template1D
0247         float cotAlpha = ltp.dxdz();
0248         float cotBeta = ltp.dydz();
0249         float locBx = 1.;
0250         if (cotBeta < 0.)
0251           locBx = -1.;
0252         float locBz = locBx;
0253         if (cotAlpha < 0.)
0254           locBz = -locBx;
0255         templ.interpolate(templateDBobject_->getTemplateID(id), cotAlpha, cotBeta, locBz, locBx);
0256         auto charge_cor = (charge * templ.qscale()) / templ.r_qMeas_qTrue();
0257         auto tmpl = templ.qscale() / templ.r_qMeas_qTrue();
0258
0259         auto clustgp = pixhit->globalPosition();  // from rechit
0260
0261         int part;
0262         ClusterData::ArrayType meas;
0263         std::pair<float, float> pred;
0264         if (shapeFilter.getSizes(*pixhit, localDir, pixelClusterShapeCache, part, meas, pred)) {
0265           auto shape = shapeFilter.isCompatible(*pixhit, localDir, pixelClusterShapeCache);
0266           unsigned shapeVal = (shape ? 1 : 0);
0267
0268           if (iAmBarrel) {
0269             if (iAmOuter) {
0270               histo[ON_TRACK_SIZE_X_OUTER].fill(pred.first, cluster.sizeX(), id, &iEvent);
0271               histo[ON_TRACK_SIZE_Y_OUTER].fill(pred.second, cluster.sizeY(), id, &iEvent);
0272               histo[ON_TRACK_SIZE_XY_OUTER].fill(cluster.sizeY(), cluster.sizeX(), id, &iEvent);
0273
0274               histo[ON_TRACK_SHAPE_OUTER].fill(shapeVal, id, &iEvent);
0275             } else {
0276               histo[ON_TRACK_SIZE_X_INNER].fill(pred.first, cluster.sizeX(), id, &iEvent);
0277               histo[ON_TRACK_SIZE_Y_INNER].fill(pred.second, cluster.sizeY(), id, &iEvent);
0278               histo[ON_TRACK_SIZE_XY_INNER].fill(cluster.sizeY(), cluster.sizeX(), id, &iEvent);
0279
0280               histo[ON_TRACK_SHAPE_INNER].fill(shapeVal, id, &iEvent);
0281             }
0282           } else {
0283             histo[ON_TRACK_SIZE_X_F].fill(pred.first, cluster.sizeX(), id, &iEvent);
0284             histo[ON_TRACK_SIZE_Y_F].fill(pred.second, cluster.sizeY(), id, &iEvent);
0285             histo[ON_TRACK_SIZE_XY_F].fill(cluster.sizeY(), cluster.sizeX(), id, &iEvent);
0286           }
0287           histo[ON_TRACK_SHAPE].fill(shapeVal, id, &iEvent);
0288         }
0289
0290         for (int i = 0; i < cluster.size(); i++) {
0291           SiPixelCluster::Pixel const& vecipxl = cluster.pixel(i);
0292           histo[DIGIS_HITMAP_ON_TRACK].fill(id, &iEvent, vecipxl.y, vecipxl.x);
0293           histo[ON_TRACK_NDIGIS].fill(id, &iEvent);
0294         }
0295
0296         histo[ON_TRACK_NCLUSTERS].fill(id, &iEvent);
0297         histo[ON_TRACK_CHARGE].fill(charge, id, &iEvent);
0298         histo[ON_TRACK_CORRECTEDCHARGE].fill(charge_cor, id, &iEvent);
0299         histo[TEMPLATE_CORRECTION].fill(tmpl, id, &iEvent);
0300         histo[ON_TRACK_SIZE].fill(cluster.size(), id, &iEvent);
0301
0302         histo[ON_TRACK_POSITIONB].fill(clustgp.z(), clustgp.phi(), id, &iEvent);
0303         histo[ON_TRACK_POSITIONF].fill(clustgp.x(), clustgp.y(), id, &iEvent);
0304
0305         histo[CHARGE_VS_SIZE_ON_TRACK].fill(cluster.size(), charge, id, &iEvent);
0306
0307         if (iAmBarrel)  // Avoid mistakes even if specification < should > handle it
0308         {
0309           if (iAmOuter) {
0310             histo[SIZE_VS_ETA_ON_TRACK_OUTER].fill(etatk, cluster.sizeY(), id, &iEvent);
0311             histo[ON_TRACK_CHARGE_OUTER].fill(charge, id, &iEvent);
0312             histo[ON_TRACK_CORRECTEDCHARGE_OUTER].fill(charge_cor, id, &iEvent);
0313             histo[TEMPLATE_CORRECTION_OUTER].fill(tmpl, id, &iEvent);
0314           } else {
0315             histo[SIZE_VS_ETA_ON_TRACK_INNER].fill(etatk, cluster.sizeY(), id, &iEvent);
0316             histo[ON_TRACK_CHARGE_INNER].fill(charge, id, &iEvent);
0317             histo[ON_TRACK_CORRECTEDCHARGE_INNER].fill(charge_cor, id, &iEvent);
0318             histo[TEMPLATE_CORRECTION_INNER].fill(tmpl, id, &iEvent);
0319           }
0320         }
0321       }
0322
0323       // statistics on tracks
0324       histo[NTRACKS].fill(1, DetId(0), &iEvent);
0325       if (isBpixtrack || isFpixtrack)
0326         histo[NTRACKS].fill(2, DetId(0), &iEvent);
0327       if (isBpixtrack)
0328         histo[NTRACKS].fill(3, DetId(0), &iEvent);
0329       if (isFpixtrack)
0330         histo[NTRACKS].fill(4, DetId(0), &iEvent);
0331
0332       if (crossesPixVol) {
0333         if (isBpixtrack || isFpixtrack)
0334           histo[NTRACKS_INVOLUME].fill(1, DetId(0), &iEvent);
0335         else
0336           histo[NTRACKS_INVOLUME].fill(0, DetId(0), &iEvent);
0337       }
0338     }
0339
0340     histo[ON_TRACK_NCLUSTERS].executePerEventHarvesting(&iEvent);
0341     histo[ON_TRACK_NDIGIS].executePerEventHarvesting(&iEvent);
0342   }
0343
0344 }  // namespace
0345
0346 #include "FWCore/Framework/interface/MakerMacros.h"
0347 DEFINE_FWK_MODULE(SiPixelPhase1TrackClusters);