DeDx/plugins/DeDxHitInfoProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    DeDxHitInfoProducer
0004 // Class:      DeDxHitInfoProducer
0005 //
0006 /**\class DeDxHitInfoProducer DeDxHitInfoProducer.cc RecoTracker/DeDx/plugins/DeDxHitInfoProducer.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  loic Quertenmont (querten)
0015 //         Created:  Mon Nov 21 14:09:02 CEST 2014
0016 // Modifications: Tamas Almos Vami (2022)
0017
0018 // system include files
0019 #include <memory>
0020
0021 #include "DataFormats/Common/interface/ValueMap.h"
0022 #include "DataFormats/TrackReco/interface/DeDxData.h"
0023 #include "DataFormats/TrackReco/interface/DeDxHit.h"
0024 #include "DataFormats/TrackReco/interface/DeDxHitInfo.h"
0025 #include "DataFormats/TrackReco/interface/Track.h"
0026 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
0027 #include "FWCore/Framework/interface/ESHandle.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/Frameworkfwd.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 #include "FWCore/Framework/interface/stream/EDProducer.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 #include "FWCore/Utilities/interface/ESGetToken.h"
0034 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0035 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0036 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0037 #include "RecoTracker/DeDx/interface/DeDxTools.h"
0038 #include "RecoTracker/DeDx/interface/GenericTruncatedAverageDeDxEstimator.h"
0039 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShape.h"
0040 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h"
0041 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0042
0043 using namespace reco;
0044 using namespace std;
0045 using namespace edm;
0046
0047 class DeDxHitInfoProducer : public edm::stream::EDProducer<> {
0048 public:
0049   explicit DeDxHitInfoProducer(const edm::ParameterSet&);
0050   ~DeDxHitInfoProducer() override;
0051
0052 private:
0053   void beginRun(edm::Run const& run, const edm::EventSetup&) override;
0054   void produce(edm::Event&, const edm::EventSetup&) override;
0055
0056   void makeCalibrationMap(const TrackerGeometry& tkGeom_);
0057   void processRec(reco::DeDxHitInfo&, const SiStripRecHit2D&, const LocalPoint&, const LocalVector&, const float&);
0058   void processHit(const TrackingRecHit* recHit,
0059                   const float trackMomentum,
0060                   const LocalVector& trackDirection,
0061                   reco::DeDxHitInfo& hitDeDxInfo,
0062                   std::vector<float>& hitMomentum,
0063                   const LocalPoint& hitLocalPos);
0064
0065   // ----------member data ---------------------------
0066   const bool usePixel_;
0067   const bool useStrip_;
0068   const float theMeVperADCPixel_;
0069   const float theMeVperADCStrip_;
0070
0071   const unsigned int minTrackHits_;
0072   const float minTrackPt_;
0073   const float minTrackPtPrescale_;
0074   const float maxTrackEta_;
0075
0076   const std::string calibrationPath_;
0077   const bool useCalibration_;
0078   const bool doShapeTest_;
0079   const bool usePixelShape_;
0080   const bool storeMomentumAtHit_;
0081
0082   const unsigned int lowPtTracksPrescalePass_, lowPtTracksPrescaleFail_;
0083   GenericTruncatedAverageDeDxEstimator lowPtTracksEstimator_;
0084   const float lowPtTracksDeDxThreshold_;
0085   const bool usePixelForPrescales_;
0086
0087   const edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0088   const edm::EDGetTokenT<SiPixelClusterShapeCache> pixShapeCacheToken_;
0089   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0090   const edm::ESGetToken<ClusterShapeHitFilter, CkfComponentsRecord> clShapeToken_;
0091   edm::ESHandle<TrackerGeometry> tkGeom_;
0092   edm::ESHandle<ClusterShapeHitFilter> clShape_;
0093   edm::Handle<SiPixelClusterShapeCache> pixShapeCache_;
0094
0095   std::vector<std::vector<float>> calibGains_;
0096   unsigned int offsetDU_;
0097
0098   uint64_t xorshift128p(uint64_t state[2]) {
0099     uint64_t x = state[0];
0100     uint64_t const y = state[1];
0101     state[0] = y;
0102     x ^= x << 23;                              // a
0103     state[1] = x ^ y ^ (x >> 17) ^ (y >> 26);  // b, c
0104     return state[1] + y;
0105   }
0106 };
0107
0108 DeDxHitInfoProducer::DeDxHitInfoProducer(const edm::ParameterSet& iConfig)
0109     : usePixel_(iConfig.getParameter<bool>("usePixel")),
0110       useStrip_(iConfig.getParameter<bool>("useStrip")),
0111       theMeVperADCPixel_(iConfig.getParameter<double>("MeVperADCPixel")),
0112       theMeVperADCStrip_(iConfig.getParameter<double>("MeVperADCStrip")),
0113       minTrackHits_(iConfig.getParameter<unsigned>("minTrackHits")),
0114       minTrackPt_(iConfig.getParameter<double>("minTrackPt")),
0115       minTrackPtPrescale_(iConfig.getParameter<double>("minTrackPtPrescale")),
0116       maxTrackEta_(iConfig.getParameter<double>("maxTrackEta")),
0117       calibrationPath_(iConfig.getParameter<string>("calibrationPath")),
0118       useCalibration_(iConfig.getParameter<bool>("useCalibration")),
0119       doShapeTest_(iConfig.getParameter<bool>("shapeTest")),
0120       usePixelShape_(not iConfig.getParameter<edm::InputTag>("clusterShapeCache").label().empty()),
0121       storeMomentumAtHit_(iConfig.getParameter<bool>("storeMomentumAtHit")),
0122       lowPtTracksPrescalePass_(iConfig.getParameter<uint32_t>("lowPtTracksPrescalePass")),
0123       lowPtTracksPrescaleFail_(iConfig.getParameter<uint32_t>("lowPtTracksPrescaleFail")),
0124       lowPtTracksEstimator_(iConfig.getParameter<edm::ParameterSet>("lowPtTracksEstimatorParameters")),
0125       lowPtTracksDeDxThreshold_(iConfig.getParameter<double>("lowPtTracksDeDxThreshold")),
0126       usePixelForPrescales_(iConfig.getParameter<bool>("usePixelForPrescales")),
0127       tracksToken_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"))),
0128       pixShapeCacheToken_(consumes<SiPixelClusterShapeCache>(iConfig.getParameter<edm::InputTag>("clusterShapeCache"))),
0129       tkGeomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()),
0130       clShapeToken_(
0131           esConsumes<ClusterShapeHitFilter, CkfComponentsRecord>(edm::ESInputTag("", "ClusterShapeHitFilter"))) {
0132   produces<reco::DeDxHitInfoCollection>();
0133   produces<reco::DeDxHitInfoAss>();
0134   produces<edm::ValueMap<int>>("prescale");
0135   if (storeMomentumAtHit_)
0136     produces<edm::ValueMap<std::vector<float>>>("momentumAtHit");
0137
0138   if (!usePixel_ && !useStrip_)
0139     edm::LogError("DeDxHitsProducer") << "No Pixel Hits NOR Strip Hits will be saved.  Running this module is useless";
0140 }
0141
0142 DeDxHitInfoProducer::~DeDxHitInfoProducer() = default;
0143
0144 // ------------ method called once each job just before starting event loop  ------------
0145 void DeDxHitInfoProducer::beginRun(edm::Run const& iRun, const edm::EventSetup& iSetup) {
0146   tkGeom_ = iSetup.getHandle(tkGeomToken_);
0147   if (useCalibration_ && calibGains_.empty()) {
0148     offsetDU_ = tkGeom_->offsetDU(GeomDetEnumerators::PixelBarrel);  //index start at the first pixel
0149
0150     deDxTools::makeCalibrationMap(calibrationPath_, *tkGeom_, calibGains_, offsetDU_);
0151   }
0152 }
0153
0154 void DeDxHitInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0155   edm::Handle<reco::TrackCollection> trackCollectionHandle;
0156   iEvent.getByToken(tracksToken_, trackCollectionHandle);
0157   const TrackCollection& trackCollection(*trackCollectionHandle.product());
0158
0159   clShape_ = iSetup.getHandle(clShapeToken_);
0160   if (usePixelShape_)
0161     pixShapeCache_ = iEvent.getHandle(pixShapeCacheToken_);
0162
0163   // creates the output collection
0164   auto resultdedxHitColl = std::make_unique<reco::DeDxHitInfoCollection>();
0165
0166   std::vector<int> indices;
0167   std::vector<int> prescales;
0168   std::vector<std::vector<float>> hitMomenta;
0169   uint64_t state[2] = {iEvent.id().event(), iEvent.id().luminosityBlock()};
0170   for (unsigned int j = 0; j < trackCollection.size(); j++) {
0171     const reco::Track& track = trackCollection[j];
0172
0173     //track selection
0174     bool passPt = (track.pt() >= minTrackPt_), passLowDeDx = false, passHighDeDx = false, pass = passPt;
0175     if (!pass && (track.pt() >= minTrackPtPrescale_)) {
0176       if (lowPtTracksPrescalePass_ > 0) {
0177         passHighDeDx = ((xorshift128p(state) % lowPtTracksPrescalePass_) == 0);
0178       }
0179       if (lowPtTracksPrescaleFail_ > 0) {
0180         passLowDeDx = ((xorshift128p(state) % lowPtTracksPrescaleFail_) == 0);
0181       }
0182       pass = passHighDeDx || passLowDeDx;
0183     }
0184     if (!pass || std::abs(track.eta()) > maxTrackEta_ || track.numberOfValidHits() < minTrackHits_) {
0185       indices.push_back(-1);
0186       continue;
0187     }
0188
0189     reco::DeDxHitInfo hitDeDxInfo;
0190     std::vector<float> hitMomentum;
0191     auto const& trajParams = track.extra()->trajParams();
0192     auto hb = track.recHitsBegin();
0193     for (unsigned int h = 0; h < track.recHitsSize(); h++) {
0194       auto recHit = *(hb + h);
0195       if (!trackerHitRTTI::isFromDet(*recHit))
0196         continue;
0197
0198       const auto& traj = trajParams[h];
0199       processHit(recHit, traj.momentum().mag(), traj.direction(), hitDeDxInfo, hitMomentum, traj.position());
0200     }
0201     assert(!storeMomentumAtHit_ || hitMomentum.size() == hitDeDxInfo.size());
0202
0203     if (!passPt) {
0204       std::vector<DeDxHit> hits;
0205       hits.reserve(hitDeDxInfo.size());
0206       for (unsigned int i = 0, n = hitDeDxInfo.size(); i < n; ++i) {
0207         if (hitDeDxInfo.detId(i).subdetId() <= 2 && usePixelForPrescales_) {
0208           hits.push_back(DeDxHit(hitDeDxInfo.charge(i) / hitDeDxInfo.pathlength(i) * theMeVperADCPixel_, 0, 0, 0));
0209         } else if (hitDeDxInfo.detId(i).subdetId() > 2) {
0210           if (doShapeTest_ && !deDxTools::shapeSelection(*hitDeDxInfo.stripCluster(i)))
0211             continue;
0212           hits.push_back(DeDxHit(hitDeDxInfo.charge(i) / hitDeDxInfo.pathlength(i) * theMeVperADCStrip_, 0, 0, 0));
0213         }
0214       }
0215
0216       // In case we have a pixel only track, but usePixelForPrescales_ is false
0217       if (hits.empty()) {
0218         indices.push_back(-1);
0219         continue;
0220       }
0221       std::sort(hits.begin(), hits.end(), std::less<DeDxHit>());
0222       if (lowPtTracksEstimator_.dedx(hits).first < lowPtTracksDeDxThreshold_) {
0223         if (passLowDeDx) {
0224           prescales.push_back(lowPtTracksPrescaleFail_);
0225         } else {
0226           indices.push_back(-1);
0227           continue;
0228         }
0229       } else {
0230         if (passHighDeDx) {
0231           prescales.push_back(lowPtTracksPrescalePass_);
0232         } else {
0233           indices.push_back(-1);
0234           continue;
0235         }
0236       }
0237     } else {
0238       prescales.push_back(1);
0239     }
0240     indices.push_back(resultdedxHitColl->size());
0241     resultdedxHitColl->push_back(hitDeDxInfo);
0242     if (storeMomentumAtHit_)
0243       hitMomenta.push_back(hitMomentum);
0244   }
0245   ///////////////////////////////////////
0246
0247   edm::OrphanHandle<reco::DeDxHitInfoCollection> dedxHitCollHandle = iEvent.put(std::move(resultdedxHitColl));
0248
0249   //create map passing the handle to the matched collection
0250   auto dedxMatch = std::make_unique<reco::DeDxHitInfoAss>(dedxHitCollHandle);
0251   reco::DeDxHitInfoAss::Filler filler(*dedxMatch);
0252   filler.insert(trackCollectionHandle, indices.begin(), indices.end());
0253   filler.fill();
0254   iEvent.put(std::move(dedxMatch));
0255
0256   auto dedxPrescale = std::make_unique<edm::ValueMap<int>>();
0257   edm::ValueMap<int>::Filler pfiller(*dedxPrescale);
0258   pfiller.insert(dedxHitCollHandle, prescales.begin(), prescales.end());
0259   pfiller.fill();
0260   iEvent.put(std::move(dedxPrescale), "prescale");
0261
0262   if (storeMomentumAtHit_) {
0263     auto dedxMomenta = std::make_unique<edm::ValueMap<std::vector<float>>>();
0264     edm::ValueMap<std::vector<float>>::Filler mfiller(*dedxMomenta);
0265     mfiller.insert(dedxHitCollHandle, hitMomenta.begin(), hitMomenta.end());
0266     mfiller.fill();
0267     iEvent.put(std::move(dedxMomenta), "momentumAtHit");
0268   }
0269 }
0270
0271 void DeDxHitInfoProducer::processRec(reco::DeDxHitInfo& hitDeDxInfo,
0272                                      const SiStripRecHit2D& recHit,
0273                                      const LocalPoint& lpos,
0274                                      const LocalVector& ldir,
0275                                      const float& cos) {
0276   uint8_t type(0);
0277   int meas;
0278   float pred;
0279   const auto& usable = clShape_->getSizes(recHit, {}, ldir, meas, pred);
0280   if (usable && meas <= int(std::abs(pred)) + 4)
0281     type |= (1 << reco::DeDxHitInfo::Complete);
0282   if (clShape_->isCompatible(recHit, ldir))
0283     type |= (1 << reco::DeDxHitInfo::Compatible);
0284   if (usable)
0285     type |= (1 << reco::DeDxHitInfo::Calibration);
0286
0287   int NSaturating(0);
0288   const auto& detId = recHit.geographicalId();
0289   const auto* detUnit = recHit.detUnit();
0290   if (detUnit == nullptr)
0291     detUnit = tkGeom_->idToDet(detId);
0292   const auto pathLen = detUnit->surface().bounds().thickness() / cos;
0293   float chargeAbs = deDxTools::getCharge(&(recHit.stripCluster()), NSaturating, *detUnit, calibGains_, offsetDU_);
0294   hitDeDxInfo.addHit(chargeAbs, pathLen, detId, lpos, type, recHit.stripCluster());
0295 }
0296
0297 void DeDxHitInfoProducer::processHit(const TrackingRecHit* recHit,
0298                                      const float trackMomentum,
0299                                      const LocalVector& trackDirection,
0300                                      reco::DeDxHitInfo& hitDeDxInfo,
0301                                      std::vector<float>& hitMomentum,
0302                                      const LocalPoint& hitLocalPos) {
0303   auto const& thit = static_cast<BaseTrackerRecHit const&>(*recHit);
0304   if (!thit.isValid())
0305     return;
0306
0307   //make sure cosine is not 0
0308   float cosine = trackDirection.z() / trackDirection.mag();
0309   float cosineAbs = std::max(0.00000001f, std::abs(cosine));
0310
0311   auto const& clus = thit.firstClusterRef();
0312   if (!clus.isValid())
0313     return;
0314
0315   const auto* detUnit = recHit->detUnit();
0316   if (detUnit == nullptr) {
0317     detUnit = tkGeom_->idToDet(thit.geographicalId());
0318   }
0319   float pathLen = detUnit->surface().bounds().thickness() / cosineAbs;
0320
0321   if (clus.isPixel()) {
0322     if (!usePixel_)
0323       return;
0324
0325     uint8_t type(0);
0326     const auto& pixelDet = *dynamic_cast<const PixelGeomDetUnit*>(detUnit);
0327     const auto& pixelRecHit = *dynamic_cast<const SiPixelRecHit*>(recHit);
0328     ClusterData data;
0329     ClusterShape().determineShape(pixelDet, clus.pixelCluster(), data);
0330     if (data.isComplete)
0331       type |= (1 << reco::DeDxHitInfo::Complete);
0332     if (usePixelShape_ && clShape_->isCompatible(pixelRecHit, trackDirection, *pixShapeCache_))
0333       type |= (1 << reco::DeDxHitInfo::Compatible);
0334     if (data.isComplete && data.isStraight && data.hasBigPixelsOnlyInside)
0335       type |= (1 << reco::DeDxHitInfo::Calibration);
0336
0337     float chargeAbs = clus.pixelCluster().charge();
0338     hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, type, clus.pixelCluster());
0339   } else if (clus.isStrip() && !thit.isMatched()) {
0340     if (!useStrip_)
0341       return;
0342
0343     processRec(hitDeDxInfo, {thit.geographicalId(), clus}, hitLocalPos, trackDirection, cosineAbs);
0344   } else if (clus.isStrip() && thit.isMatched()) {
0345     if (!useStrip_)
0346       return;
0347     const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0348     if (!matchedHit)
0349       return;
0350     if (storeMomentumAtHit_)
0351       hitMomentum.push_back(trackMomentum);
0352
0353     processRec(hitDeDxInfo, matchedHit->monoHit(), hitLocalPos, trackDirection, cosineAbs);
0354     processRec(hitDeDxInfo, matchedHit->stereoHit(), hitLocalPos, trackDirection, cosineAbs);
0355   }
0356   if (storeMomentumAtHit_ && (clus.isPixel() || clus.isStrip()))
0357     hitMomentum.push_back(trackMomentum);
0358 }
0359
0360 //define this as a plug-in
0361 DEFINE_FWK_MODULE(DeDxHitInfoProducer);