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_;
0092     std::vector<SiPixelTemplateStore> thePixelTemp_;
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   };
0105
0106   SiPixelPhase1TrackClusters::SiPixelPhase1TrackClusters(const edm::ParameterSet& iConfig)
0107       : SiPixelPhase1Base(iConfig), applyVertexCut_(iConfig.getUntrackedParameter<bool>("VertexCut", true)) {
0108     tracksToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0109
0110     offlinePrimaryVerticesToken_ =
0111         applyVertexCut_ ? consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))
0112                         : edm::EDGetTokenT<reco::VertexCollection>();
0113
0114     pixelClusterShapeCacheToken_ =
0115         consumes<SiPixelClusterShapeCache>(iConfig.getParameter<edm::InputTag>("clusterShapeCache"));
0116
0117     trackerTopoToken_ = esConsumes<edm::Transition::BeginRun>();
0118     trackerGeomToken_ = esConsumes<edm::Transition::BeginRun>();
0119     clusterShapeHitFilterToken_ =
0120         esConsumes<ClusterShapeHitFilter, CkfComponentsRecord>(edm::ESInputTag("", "ClusterShapeHitFilter"));
0121     templateDBobjectToken_ = esConsumes<edm::Transition::BeginRun>();
0122   }
0123
0124   void SiPixelPhase1TrackClusters::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0125     //get topology
0126     tkTpl = &iSetup.getData(trackerTopoToken_);
0127
0128     // Initialize 1D templates
0129     templateDBobject_ = &iSetup.getData(templateDBobjectToken_);
0130     if (!SiPixelTemplate::pushfile(*templateDBobject_, thePixelTemp_))
0131       edm::LogError("SiPixelPhase1TrackClusters")
0132           << "Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
0133           << (*templateDBobject_).version() << std::endl;
0134   }
0135
0136   void SiPixelPhase1TrackClusters::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0137     if (!checktrigger(iEvent, iSetup, DCS))
0138       return;
0139
0140     if (histo.size() != ENUM_SIZE) {
0141       edm::LogError("SiPixelPhase1TrackClusters")
0142           << "incompatible configuration " << histo.size() << "!=" << ENUM_SIZE << std::endl;
0143       return;
0144     }
0145
0146     edm::ESHandle<ClusterShapeHitFilter> shapeFilterH = iSetup.getHandle(clusterShapeHitFilterToken_);
0147     auto const& shapeFilter = *shapeFilterH;
0148
0149     edm::Handle<reco::VertexCollection> vertices;
0150     if (applyVertexCut_) {
0151       iEvent.getByToken(offlinePrimaryVerticesToken_, vertices);
0152       if (!vertices.isValid() || vertices->empty())
0153         return;
0154     }
0155
0156     //get the map
0157     edm::Handle<reco::TrackCollection> tracks;
0158     iEvent.getByToken(tracksToken_, tracks);
0159
0160     if (!tracks.isValid()) {
0161       edm::LogWarning("SiPixelPhase1TrackClusters") << "track collection is not valid";
0162       return;
0163     }
0164
0165     SiPixelTemplate templ(thePixelTemp_);
0166
0167     edm::Handle<SiPixelClusterShapeCache> pixelClusterShapeCacheH;
0168     iEvent.getByToken(pixelClusterShapeCacheToken_, pixelClusterShapeCacheH);
0169     if (!pixelClusterShapeCacheH.isValid()) {
0170       edm::LogWarning("SiPixelPhase1TrackClusters") << "PixelClusterShapeCache collection is not valid";
0171       return;
0172     }
0173     auto const& pixelClusterShapeCache = *pixelClusterShapeCacheH;
0174
0175     for (auto const& track : *tracks) {
0176       if (applyVertexCut_ &&
0177           (track.pt() < 0.75 || std::abs(track.dxy((*vertices)[0].position())) > 5 * track.dxyError()))
0178         continue;
0179
0180       bool isBpixtrack = false, isFpixtrack = false, crossesPixVol = false;
0181
0182       // find out whether track crosses pixel fiducial volume (for cosmic tracks)
0183       auto d0 = track.d0(), dz = track.dz();
0184       if (std::abs(d0) < 16 && std::abs(dz) < 50)
0185         crossesPixVol = true;
0186
0187       auto etatk = track.eta();
0188
0189       auto const& trajParams = track.extra()->trajParams();
0190       assert(trajParams.size() == track.recHitsSize());
0191       auto hb = track.recHitsBegin();
0192
0193       for (unsigned int h = 0; h < track.recHitsSize(); h++) {
0194         auto hit = *(hb + h);
0195         if (!hit->isValid())
0196           continue;
0197         auto id = hit->geographicalId();
0198
0199         // check that we are in the pixel
0200         auto subdetid = (id.subdetId());
0201         if (subdetid == PixelSubdetector::PixelBarrel)
0202           isBpixtrack = true;
0203         if (subdetid == PixelSubdetector::PixelEndcap)
0204           isFpixtrack = true;
0205         if (subdetid != PixelSubdetector::PixelBarrel && subdetid != PixelSubdetector::PixelEndcap)
0206           continue;
0207         bool iAmBarrel = subdetid == PixelSubdetector::PixelBarrel;
0208
0209         // PXB_L4 IS IN THE OTHER WAY
0210         // CAN BE XORed BUT LETS KEEP THINGS SIMPLE
0211         bool iAmOuter = ((tkTpl->pxbLadder(id) % 2 == 1) && tkTpl->pxbLayer(id) != 4) ||
0212                         ((tkTpl->pxbLadder(id) % 2 != 1) && tkTpl->pxbLayer(id) == 4);
0213
0214         auto pixhit = dynamic_cast<const SiPixelRecHit*>(hit->hit());
0215         if (!pixhit)
0216           continue;
0217
0218         auto geomdetunit = dynamic_cast<const PixelGeomDetUnit*>(pixhit->detUnit());
0219         auto const& topol = geomdetunit->specificTopology();
0220
0221         // get the cluster
0222         auto clustp = pixhit->cluster();
0223         if (clustp.isNull())
0224           continue;
0225         auto const& cluster = *clustp;
0226         const std::vector<SiPixelCluster::Pixel> pixelsVec = cluster.pixels();
0227         for (unsigned int i = 0; i < pixelsVec.size(); ++i) {
0228           float pixx = pixelsVec[i].x;  // index as float=iteger, row index
0229           float pixy = pixelsVec[i].y;  // same, col index
0230
0231           bool bigInX = topol.isItBigPixelInX(int(pixx));
0232           bool bigInY = topol.isItBigPixelInY(int(pixy));
0233           float pixel_charge = pixelsVec[i].adc;
0234
0235           if (bigInX == true || bigInY == true) {
0236             histo[ON_TRACK_BIGPIXELCHARGE].fill(pixel_charge, id, &iEvent);
0237           } else {
0238             histo[ON_TRACK_NOTBIGPIXELCHARGE].fill(pixel_charge, id, &iEvent);
0239           }
0240         }  // End loop over pixels
0241         auto const& ltp = trajParams[h];
0242
0243         auto localDir = ltp.momentum() / ltp.momentum().mag();
0244
0245         // correct charge for track impact angle
0246         auto charge = cluster.charge() * ltp.absdz();
0247         //Correct charge with Template1D
0248         float cotAlpha = ltp.dxdz();
0249         float cotBeta = ltp.dydz();
0250         float locBx = 1.;
0251         if (cotBeta < 0.)
0252           locBx = -1.;
0253         float locBz = locBx;
0254         if (cotAlpha < 0.)
0255           locBz = -locBx;
0256         templ.interpolate(templateDBobject_->getTemplateID(id), cotAlpha, cotBeta, locBz, locBx);
0257         auto charge_cor = (charge * templ.qscale()) / templ.r_qMeas_qTrue();
0258         auto tmpl = templ.qscale() / templ.r_qMeas_qTrue();
0259
0260         auto clustgp = pixhit->globalPosition();  // from rechit
0261
0262         int part;
0263         ClusterData::ArrayType meas;
0264         std::pair<float, float> pred;
0265         if (shapeFilter.getSizes(*pixhit, localDir, pixelClusterShapeCache, part, meas, pred)) {
0266           auto shape = shapeFilter.isCompatible(*pixhit, localDir, pixelClusterShapeCache);
0267           unsigned shapeVal = (shape ? 1 : 0);
0268
0269           if (iAmBarrel) {
0270             if (iAmOuter) {
0271               histo[ON_TRACK_SIZE_X_OUTER].fill(pred.first, cluster.sizeX(), id, &iEvent);
0272               histo[ON_TRACK_SIZE_Y_OUTER].fill(pred.second, cluster.sizeY(), id, &iEvent);
0273               histo[ON_TRACK_SIZE_XY_OUTER].fill(cluster.sizeY(), cluster.sizeX(), id, &iEvent);
0274
0275               histo[ON_TRACK_SHAPE_OUTER].fill(shapeVal, id, &iEvent);
0276             } else {
0277               histo[ON_TRACK_SIZE_X_INNER].fill(pred.first, cluster.sizeX(), id, &iEvent);
0278               histo[ON_TRACK_SIZE_Y_INNER].fill(pred.second, cluster.sizeY(), id, &iEvent);
0279               histo[ON_TRACK_SIZE_XY_INNER].fill(cluster.sizeY(), cluster.sizeX(), id, &iEvent);
0280
0281               histo[ON_TRACK_SHAPE_INNER].fill(shapeVal, id, &iEvent);
0282             }
0283           } else {
0284             histo[ON_TRACK_SIZE_X_F].fill(pred.first, cluster.sizeX(), id, &iEvent);
0285             histo[ON_TRACK_SIZE_Y_F].fill(pred.second, cluster.sizeY(), id, &iEvent);
0286             histo[ON_TRACK_SIZE_XY_F].fill(cluster.sizeY(), cluster.sizeX(), id, &iEvent);
0287           }
0288           histo[ON_TRACK_SHAPE].fill(shapeVal, id, &iEvent);
0289         }
0290
0291         for (int i = 0; i < cluster.size(); i++) {
0292           SiPixelCluster::Pixel const& vecipxl = cluster.pixel(i);
0293           histo[DIGIS_HITMAP_ON_TRACK].fill(id, &iEvent, vecipxl.y, vecipxl.x);
0294           histo[ON_TRACK_NDIGIS].fill(id, &iEvent);
0295         }
0296
0297         histo[ON_TRACK_NCLUSTERS].fill(id, &iEvent);
0298         histo[ON_TRACK_CHARGE].fill(charge, id, &iEvent);
0299         histo[ON_TRACK_CORRECTEDCHARGE].fill(charge_cor, id, &iEvent);
0300         histo[TEMPLATE_CORRECTION].fill(tmpl, id, &iEvent);
0301         histo[ON_TRACK_SIZE].fill(cluster.size(), id, &iEvent);
0302
0303         histo[ON_TRACK_POSITIONB].fill(clustgp.z(), clustgp.phi(), id, &iEvent);
0304         histo[ON_TRACK_POSITIONF].fill(clustgp.x(), clustgp.y(), id, &iEvent);
0305
0306         histo[CHARGE_VS_SIZE_ON_TRACK].fill(cluster.size(), charge, id, &iEvent);
0307
0308         if (iAmBarrel)  // Avoid mistakes even if specification < should > handle it
0309         {
0310           if (iAmOuter) {
0311             histo[SIZE_VS_ETA_ON_TRACK_OUTER].fill(etatk, cluster.sizeY(), id, &iEvent);
0312             histo[ON_TRACK_CHARGE_OUTER].fill(charge, id, &iEvent);
0313             histo[ON_TRACK_CORRECTEDCHARGE_OUTER].fill(charge_cor, id, &iEvent);
0314             histo[TEMPLATE_CORRECTION_OUTER].fill(tmpl, id, &iEvent);
0315           } else {
0316             histo[SIZE_VS_ETA_ON_TRACK_INNER].fill(etatk, cluster.sizeY(), id, &iEvent);
0317             histo[ON_TRACK_CHARGE_INNER].fill(charge, id, &iEvent);
0318             histo[ON_TRACK_CORRECTEDCHARGE_INNER].fill(charge_cor, id, &iEvent);
0319             histo[TEMPLATE_CORRECTION_INNER].fill(tmpl, id, &iEvent);
0320           }
0321         }
0322       }
0323
0324       // statistics on tracks
0325       histo[NTRACKS].fill(1, DetId(0), &iEvent);
0326       if (isBpixtrack || isFpixtrack)
0327         histo[NTRACKS].fill(2, DetId(0), &iEvent);
0328       if (isBpixtrack)
0329         histo[NTRACKS].fill(3, DetId(0), &iEvent);
0330       if (isFpixtrack)
0331         histo[NTRACKS].fill(4, DetId(0), &iEvent);
0332
0333       if (crossesPixVol) {
0334         if (isBpixtrack || isFpixtrack)
0335           histo[NTRACKS_INVOLUME].fill(1, DetId(0), &iEvent);
0336         else
0337           histo[NTRACKS_INVOLUME].fill(0, DetId(0), &iEvent);
0338       }
0339     }
0340
0341     histo[ON_TRACK_NCLUSTERS].executePerEventHarvesting(&iEvent);
0342     histo[ON_TRACK_NDIGIS].executePerEventHarvesting(&iEvent);
0343   }
0344
0345 }  // namespace
0346
0347 #include "FWCore/Framework/interface/MakerMacros.h"
0348 DEFINE_FWK_MODULE(SiPixelPhase1TrackClusters);