File indexing completed on 2023-03-17 11:22:14
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
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
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;
0093 state[1] = x ^ y ^ (x >> 17) ^ (y >> 26);
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
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);
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
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
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
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
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
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
0302 DEFINE_FWK_MODULE(DeDxHitInfoProducer);