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 "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0040
0041 using namespace reco;
0042 using namespace std;
0043 using namespace edm;
0044
0045 class DeDxHitInfoProducer : public edm::stream::EDProducer<> {
0046 public:
0047   explicit DeDxHitInfoProducer(const edm::ParameterSet&);
0048   ~DeDxHitInfoProducer() override;
0049
0050 private:
0051   void beginRun(edm::Run const& run, const edm::EventSetup&) override;
0052   void produce(edm::Event&, const edm::EventSetup&) override;
0053
0054   void makeCalibrationMap(const TrackerGeometry& tkGeom_);
0055   void processHit(const TrackingRecHit* recHit,
0056                   const float trackMomentum,
0057                   const float cosine,
0058                   reco::DeDxHitInfo& hitDeDxInfo,
0059                   const LocalPoint& hitLocalPos);
0060
0061   // ----------member data ---------------------------
0062   const bool usePixel_;
0063   const bool useStrip_;
0064   const float theMeVperADCPixel_;
0065   const float theMeVperADCStrip_;
0066
0067   const unsigned int minTrackHits_;
0068   const float minTrackPt_;
0069   const float minTrackPtPrescale_;
0070   const float maxTrackEta_;
0071
0072   const std::string calibrationPath_;
0073   const bool useCalibration_;
0074   const bool doShapeTest_;
0075
0076   const unsigned int lowPtTracksPrescalePass_, lowPtTracksPrescaleFail_;
0077   GenericTruncatedAverageDeDxEstimator lowPtTracksEstimator_;
0078   const float lowPtTracksDeDxThreshold_;
0079   const bool usePixelForPrescales_;
0080
0081   const edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0082   edm::ESHandle<TrackerGeometry> tkGeom_;
0083   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0084
0085   std::vector<std::vector<float>> calibGains_;
0086   unsigned int offsetDU_;
0087
0088   uint64_t xorshift128p(uint64_t state[2]) {
0089     uint64_t x = state[0];
0090     uint64_t const y = state[1];
0091     state[0] = y;
0092     x ^= x << 23;                              // a
0093     state[1] = x ^ y ^ (x >> 17) ^ (y >> 26);  // b, c
0094     return state[1] + y;
0095   }
0096 };
0097
0098 DeDxHitInfoProducer::DeDxHitInfoProducer(const edm::ParameterSet& iConfig)
0099     : usePixel_(iConfig.getParameter<bool>("usePixel")),
0100       useStrip_(iConfig.getParameter<bool>("useStrip")),
0101       theMeVperADCPixel_(iConfig.getParameter<double>("MeVperADCPixel")),
0102       theMeVperADCStrip_(iConfig.getParameter<double>("MeVperADCStrip")),
0103       minTrackHits_(iConfig.getParameter<unsigned>("minTrackHits")),
0104       minTrackPt_(iConfig.getParameter<double>("minTrackPt")),
0105       minTrackPtPrescale_(iConfig.getParameter<double>("minTrackPtPrescale")),
0106       maxTrackEta_(iConfig.getParameter<double>("maxTrackEta")),
0107       calibrationPath_(iConfig.getParameter<string>("calibrationPath")),
0108       useCalibration_(iConfig.getParameter<bool>("useCalibration")),
0109       doShapeTest_(iConfig.getParameter<bool>("shapeTest")),
0110       lowPtTracksPrescalePass_(iConfig.getParameter<uint32_t>("lowPtTracksPrescalePass")),
0111       lowPtTracksPrescaleFail_(iConfig.getParameter<uint32_t>("lowPtTracksPrescaleFail")),
0112       lowPtTracksEstimator_(iConfig.getParameter<edm::ParameterSet>("lowPtTracksEstimatorParameters")),
0113       lowPtTracksDeDxThreshold_(iConfig.getParameter<double>("lowPtTracksDeDxThreshold")),
0114       usePixelForPrescales_(iConfig.getParameter<bool>("usePixelForPrescales")),
0115       tracksToken_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"))),
0116       tkGeomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()) {
0117   produces<reco::DeDxHitInfoCollection>();
0118   produces<reco::DeDxHitInfoAss>();
0119   produces<edm::ValueMap<int>>("prescale");
0120
0121   if (!usePixel_ && !useStrip_)
0122     edm::LogError("DeDxHitsProducer") << "No Pixel Hits NOR Strip Hits will be saved.  Running this module is useless";
0123 }
0124
0125 DeDxHitInfoProducer::~DeDxHitInfoProducer() = default;
0126
0127 // ------------ method called once each job just before starting event loop  ------------
0128 void DeDxHitInfoProducer::beginRun(edm::Run const& iRun, const edm::EventSetup& iSetup) {
0129   tkGeom_ = iSetup.getHandle(tkGeomToken_);
0130   if (useCalibration_ && calibGains_.empty()) {
0131     offsetDU_ = tkGeom_->offsetDU(GeomDetEnumerators::PixelBarrel);  //index start at the first pixel
0132
0133     deDxTools::makeCalibrationMap(calibrationPath_, *tkGeom_, calibGains_, offsetDU_);
0134   }
0135 }
0136
0137 void DeDxHitInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0138   edm::Handle<reco::TrackCollection> trackCollectionHandle;
0139   iEvent.getByToken(tracksToken_, trackCollectionHandle);
0140   const TrackCollection& trackCollection(*trackCollectionHandle.product());
0141
0142   // creates the output collection
0143   auto resultdedxHitColl = std::make_unique<reco::DeDxHitInfoCollection>();
0144
0145   std::vector<int> indices;
0146   std::vector<int> prescales;
0147   uint64_t state[2] = {iEvent.id().event(), iEvent.id().luminosityBlock()};
0148   for (unsigned int j = 0; j < trackCollection.size(); j++) {
0149     const reco::Track& track = trackCollection[j];
0150
0151     //track selection
0152     bool passPt = (track.pt() >= minTrackPt_), passLowDeDx = false, passHighDeDx = false, pass = passPt;
0153     if (!pass && (track.pt() >= minTrackPtPrescale_)) {
0154       if (lowPtTracksPrescalePass_ > 0) {
0155         passHighDeDx = ((xorshift128p(state) % lowPtTracksPrescalePass_) == 0);
0156       }
0157       if (lowPtTracksPrescaleFail_ > 0) {
0158         passLowDeDx = ((xorshift128p(state) % lowPtTracksPrescaleFail_) == 0);
0159       }
0160       pass = passHighDeDx || passLowDeDx;
0161     }
0162     if (!pass || std::abs(track.eta()) > maxTrackEta_ || track.numberOfValidHits() < minTrackHits_) {
0163       indices.push_back(-1);
0164       continue;
0165     }
0166
0167     reco::DeDxHitInfo hitDeDxInfo;
0168     auto const& trajParams = track.extra()->trajParams();
0169     auto hb = track.recHitsBegin();
0170     for (unsigned int h = 0; h < track.recHitsSize(); h++) {
0171       auto recHit = *(hb + h);
0172       if (!trackerHitRTTI::isFromDet(*recHit))
0173         continue;
0174
0175       auto trackDirection = trajParams[h].direction();
0176       float cosine = trackDirection.z() / trackDirection.mag();
0177
0178       processHit(recHit, track.p(), cosine, hitDeDxInfo, trajParams[h].position());
0179     }
0180
0181     if (!passPt) {
0182       std::vector<DeDxHit> hits;
0183       hits.reserve(hitDeDxInfo.size());
0184       for (unsigned int i = 0, n = hitDeDxInfo.size(); i < n; ++i) {
0185         if (hitDeDxInfo.detId(i).subdetId() <= 2 && usePixelForPrescales_) {
0186           hits.push_back(DeDxHit(hitDeDxInfo.charge(i) / hitDeDxInfo.pathlength(i) * theMeVperADCPixel_, 0, 0, 0));
0187         } else if (hitDeDxInfo.detId(i).subdetId() > 2) {
0188           if (doShapeTest_ && !deDxTools::shapeSelection(*hitDeDxInfo.stripCluster(i)))
0189             continue;
0190           hits.push_back(DeDxHit(hitDeDxInfo.charge(i) / hitDeDxInfo.pathlength(i) * theMeVperADCStrip_, 0, 0, 0));
0191         }
0192       }
0193
0194       // In case we have a pixel only track, but usePixelForPrescales_ is false
0195       if (hits.empty()) {
0196         indices.push_back(-1);
0197         continue;
0198       }
0199       std::sort(hits.begin(), hits.end(), std::less<DeDxHit>());
0200       if (lowPtTracksEstimator_.dedx(hits).first < lowPtTracksDeDxThreshold_) {
0201         if (passLowDeDx) {
0202           prescales.push_back(lowPtTracksPrescaleFail_);
0203         } else {
0204           indices.push_back(-1);
0205           continue;
0206         }
0207       } else {
0208         if (passHighDeDx) {
0209           prescales.push_back(lowPtTracksPrescalePass_);
0210         } else {
0211           indices.push_back(-1);
0212           continue;
0213         }
0214       }
0215     } else {
0216       prescales.push_back(1);
0217     }
0218     indices.push_back(resultdedxHitColl->size());
0219     resultdedxHitColl->push_back(hitDeDxInfo);
0220   }
0221   ///////////////////////////////////////
0222
0223   edm::OrphanHandle<reco::DeDxHitInfoCollection> dedxHitCollHandle = iEvent.put(std::move(resultdedxHitColl));
0224
0225   //create map passing the handle to the matched collection
0226   auto dedxMatch = std::make_unique<reco::DeDxHitInfoAss>(dedxHitCollHandle);
0227   reco::DeDxHitInfoAss::Filler filler(*dedxMatch);
0228   filler.insert(trackCollectionHandle, indices.begin(), indices.end());
0229   filler.fill();
0230   iEvent.put(std::move(dedxMatch));
0231
0232   auto dedxPrescale = std::make_unique<edm::ValueMap<int>>();
0233   edm::ValueMap<int>::Filler pfiller(*dedxPrescale);
0234   pfiller.insert(dedxHitCollHandle, prescales.begin(), prescales.end());
0235   pfiller.fill();
0236   iEvent.put(std::move(dedxPrescale), "prescale");
0237 }
0238
0239 void DeDxHitInfoProducer::processHit(const TrackingRecHit* recHit,
0240                                      const float trackMomentum,
0241                                      const float cosine,
0242                                      reco::DeDxHitInfo& hitDeDxInfo,
0243                                      const LocalPoint& hitLocalPos) {
0244   auto const& thit = static_cast<BaseTrackerRecHit const&>(*recHit);
0245   if (!thit.isValid())
0246     return;
0247
0248   //make sure cosine is not 0
0249   float cosineAbs = std::max(0.00000001f, std::abs(cosine));
0250
0251   auto const& clus = thit.firstClusterRef();
0252   if (!clus.isValid())
0253     return;
0254
0255   const auto* detUnit = recHit->detUnit();
0256   if (detUnit == nullptr) {
0257     detUnit = tkGeom_->idToDet(thit.geographicalId());
0258   }
0259   float pathLen = detUnit->surface().bounds().thickness() / cosineAbs;
0260
0261   if (clus.isPixel()) {
0262     if (!usePixel_)
0263       return;
0264
0265     float chargeAbs = clus.pixelCluster().charge();
0266     hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, clus.pixelCluster());
0267   } else if (clus.isStrip() && !thit.isMatched()) {
0268     if (!useStrip_)
0269       return;
0270
0271     int NSaturating = 0;
0272     float chargeAbs = deDxTools::getCharge(&(clus.stripCluster()), NSaturating, *detUnit, calibGains_, offsetDU_);
0273     hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, clus.stripCluster());
0274   } else if (clus.isStrip() && thit.isMatched()) {
0275     if (!useStrip_)
0276       return;
0277     const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0278     if (!matchedHit)
0279       return;
0280
0281     const auto* detUnitM = matchedHit->monoHit().detUnit();
0282     if (detUnitM == nullptr)
0283       detUnitM = tkGeom_->idToDet(matchedHit->monoHit().geographicalId());
0284     int NSaturating = 0;
0285     auto pathLenM = detUnitM->surface().bounds().thickness() / cosineAbs;
0286     float chargeAbs =
0287         deDxTools::getCharge(&(matchedHit->monoHit().stripCluster()), NSaturating, *detUnitM, calibGains_, offsetDU_);
0288     hitDeDxInfo.addHit(chargeAbs, pathLenM, thit.geographicalId(), hitLocalPos, matchedHit->monoHit().stripCluster());
0289
0290     const auto* detUnitS = matchedHit->stereoHit().detUnit();
0291     if (detUnitS == nullptr)
0292       detUnitS = tkGeom_->idToDet(matchedHit->stereoHit().geographicalId());
0293     NSaturating = 0;
0294     auto pathLenS = detUnitS->surface().bounds().thickness() / cosineAbs;
0295     chargeAbs =
0296         deDxTools::getCharge(&(matchedHit->stereoHit().stripCluster()), NSaturating, *detUnitS, calibGains_, offsetDU_);
0297     hitDeDxInfo.addHit(chargeAbs, pathLenS, thit.geographicalId(), hitLocalPos, matchedHit->stereoHit().stripCluster());
0298   }
0299 }
0300
0301 //define this as a plug-in
0302 DEFINE_FWK_MODULE(DeDxHitInfoProducer);