File indexing completed on 2024-07-16 22:52:37
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 #include "RecoTracker/DeDx/plugins/DeDxEstimatorProducer.h"
0023 #include "FWCore/Framework/interface/ConsumesCollector.h"
0024
0025 using namespace reco;
0026 using namespace std;
0027 using namespace edm;
0028
0029 void DeDxEstimatorProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0030 edm::ParameterSetDescription desc;
0031 desc.add<string>("estimator", "generic");
0032 desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0033 desc.add<bool>("UseDeDxHits", false);
0034 desc.add<edm::InputTag>("pixelDeDxHits", {});
0035 desc.add<edm::InputTag>("stripDeDxHits", {});
0036 desc.add<bool>("UsePixel", false);
0037 desc.add<bool>("UseStrip", true);
0038 desc.add<double>("MeVperADCPixel", 3.61e-06);
0039 desc.add<double>("MeVperADCStrip", 3.61e-06 * 265);
0040 desc.add<bool>("ShapeTest", true);
0041 desc.add<bool>("UseCalibration", false);
0042 desc.add<string>("calibrationPath", "");
0043 desc.add<string>("Record", "SiStripDeDxMip_3D_Rcd");
0044 desc.add<string>("ProbabilityMode", "Accumulation");
0045 desc.add<double>("fraction", 0.4);
0046 desc.add<double>("exponent", -2.0);
0047 desc.add<bool>("truncate", true);
0048
0049 descriptions.add("DeDxEstimatorProducer", desc);
0050 }
0051
0052 DeDxEstimatorProducer::DeDxEstimatorProducer(const edm::ParameterSet& iConfig)
0053 : tkGeomToken(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()) {
0054 produces<ValueMap<DeDxData>>();
0055
0056 auto cCollector = consumesCollector();
0057 string estimatorName = iConfig.getParameter<string>("estimator");
0058 if (estimatorName == "median")
0059 m_estimator = std::make_unique<MedianDeDxEstimator>(iConfig);
0060 else if (estimatorName == "generic")
0061 m_estimator = std::make_unique<GenericAverageDeDxEstimator>(iConfig);
0062 else if (estimatorName == "truncated")
0063 m_estimator = std::make_unique<TruncatedAverageDeDxEstimator>(iConfig);
0064 else if (estimatorName == "genericTruncated")
0065 m_estimator = std::make_unique<GenericTruncatedAverageDeDxEstimator>(iConfig);
0066 else if (estimatorName == "unbinnedFit")
0067 m_estimator = std::make_unique<UnbinnedFitDeDxEstimator>(iConfig);
0068 else if (estimatorName == "likelihoodFit")
0069 m_estimator = std::make_unique<LikelihoodFitDeDxEstimator>(iConfig);
0070 else if (estimatorName == "productDiscrim")
0071 m_estimator = std::make_unique<ProductDeDxDiscriminator>(iConfig, cCollector);
0072 else if (estimatorName == "btagDiscrim")
0073 m_estimator = std::make_unique<BTagLikeDeDxDiscriminator>(iConfig, cCollector);
0074 else if (estimatorName == "smirnovDiscrim")
0075 m_estimator = std::make_unique<SmirnovDeDxDiscriminator>(iConfig, cCollector);
0076 else if (estimatorName == "asmirnovDiscrim")
0077 m_estimator = std::make_unique<ASmirnovDeDxDiscriminator>(iConfig, cCollector);
0078
0079
0080
0081
0082 m_tracksTag = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0083
0084 useDeDxHits = iConfig.getParameter<bool>("UseDeDxHits");
0085 m_pixelDeDxHitsTag = consumes<reco::TrackDeDxHitsCollection>(iConfig.getParameter<edm::InputTag>("pixelDeDxHits"));
0086 m_stripDeDxHitsTag = consumes<reco::TrackDeDxHitsCollection>(iConfig.getParameter<edm::InputTag>("stripDeDxHits"));
0087
0088 usePixel = iConfig.getParameter<bool>("UsePixel");
0089 useStrip = iConfig.getParameter<bool>("UseStrip");
0090 meVperADCPixel = iConfig.getParameter<double>("MeVperADCPixel");
0091 meVperADCStrip = iConfig.getParameter<double>("MeVperADCStrip");
0092
0093 shapetest = iConfig.getParameter<bool>("ShapeTest");
0094 useCalibration = iConfig.getParameter<bool>("UseCalibration");
0095 m_calibrationPath = iConfig.getParameter<string>("calibrationPath");
0096
0097 if (!usePixel && !useStrip)
0098 edm::LogWarning("DeDxHitsProducer")
0099 << "Pixel Hits AND Strip Hits will not be used to estimate dEdx --> BUG, Please Update the config file";
0100 }
0101
0102 DeDxEstimatorProducer::~DeDxEstimatorProducer() {}
0103
0104
0105 void DeDxEstimatorProducer::beginRun(edm::Run const& run, const edm::EventSetup& iSetup) {
0106 tkGeom = &iSetup.getData(tkGeomToken);
0107 if (useCalibration && calibGains.empty()) {
0108 m_off = tkGeom->offsetDU(GeomDetEnumerators::PixelBarrel);
0109
0110 deDxTools::makeCalibrationMap(m_calibrationPath, *tkGeom, calibGains, m_off);
0111 }
0112
0113 m_estimator->beginRun(run, iSetup);
0114 }
0115
0116 void DeDxEstimatorProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0117 auto trackDeDxEstimateAssociation = std::make_unique<ValueMap<DeDxData>>();
0118 ValueMap<DeDxData>::Filler filler(*trackDeDxEstimateAssociation);
0119
0120 edm::Handle<reco::TrackCollection> trackCollectionHandle;
0121 iEvent.getByToken(m_tracksTag, trackCollectionHandle);
0122 const auto& pixelDeDxHits = iEvent.getHandle(m_pixelDeDxHitsTag);
0123 const auto& stripDeDxHits = iEvent.getHandle(m_stripDeDxHitsTag);
0124
0125 std::vector<DeDxData> dedxEstimate(trackCollectionHandle->size());
0126
0127 for (unsigned int j = 0; j < trackCollectionHandle->size(); j++) {
0128 const reco::TrackRef track = reco::TrackRef(trackCollectionHandle, j);
0129
0130 int NClusterSaturating = 0;
0131 DeDxHitCollection dedxHits;
0132
0133 if (useDeDxHits) {
0134 if (usePixel)
0135 dedxHits = (*pixelDeDxHits)[track];
0136 if (useStrip) {
0137 const auto& hits = (*stripDeDxHits)[track];
0138 dedxHits.insert(dedxHits.end(), hits.begin(), hits.end());
0139 }
0140 } else {
0141 auto const& trajParams = track->extra()->trajParams();
0142 assert(trajParams.size() == track->recHitsSize());
0143 auto hb = track->recHitsBegin();
0144 dedxHits.reserve(track->recHitsSize() / 2);
0145 for (unsigned int h = 0; h < track->recHitsSize(); h++) {
0146 auto recHit = *(hb + h);
0147 if (!trackerHitRTTI::isFromDet(*recHit))
0148 continue;
0149
0150 auto trackDirection = trajParams[h].direction();
0151 float cosine = trackDirection.z() / trackDirection.mag();
0152 processHit(recHit, track->p(), cosine, dedxHits, NClusterSaturating);
0153 }
0154 }
0155
0156 sort(dedxHits.begin(), dedxHits.end(), less<DeDxHit>());
0157 std::pair<float, float> val_and_error = m_estimator->dedx(dedxHits);
0158
0159
0160
0161 val_and_error.second = NClusterSaturating;
0162 dedxEstimate[j] = DeDxData(val_and_error.first, val_and_error.second, dedxHits.size());
0163 }
0164
0165
0166 filler.insert(trackCollectionHandle, dedxEstimate.begin(), dedxEstimate.end());
0167
0168
0169 filler.fill();
0170 iEvent.put(std::move(trackDeDxEstimateAssociation));
0171 }
0172
0173 void DeDxEstimatorProducer::processHit(const TrackingRecHit* recHit,
0174 float trackMomentum,
0175 float& cosine,
0176 reco::DeDxHitCollection& dedxHits,
0177 int& NClusterSaturating) {
0178 auto const& thit = static_cast<BaseTrackerRecHit const&>(*recHit);
0179 if (!thit.isValid())
0180 return;
0181
0182 auto const& clus = thit.firstClusterRef();
0183 if (!clus.isValid())
0184 return;
0185
0186 if (clus.isPixel()) {
0187 if (!usePixel)
0188 return;
0189
0190 const auto* detUnit = recHit->detUnit();
0191 if (detUnit == nullptr)
0192 detUnit = tkGeom->idToDet(thit.geographicalId());
0193 float pathLen = detUnit->surface().bounds().thickness() / fabs(cosine);
0194 float chargeAbs = clus.pixelCluster().charge();
0195 float charge = meVperADCPixel * chargeAbs / pathLen;
0196 dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, thit.geographicalId()));
0197 } else if (clus.isStrip() && !thit.isMatched()) {
0198 if (!useStrip)
0199 return;
0200
0201 const auto* detUnit = recHit->detUnit();
0202 if (detUnit == nullptr)
0203 detUnit = tkGeom->idToDet(thit.geographicalId());
0204 int NSaturating = 0;
0205 float pathLen = detUnit->surface().bounds().thickness() / fabs(cosine);
0206 float chargeAbs = deDxTools::getCharge(&(clus.stripCluster()), NSaturating, *detUnit, calibGains, m_off);
0207 float charge = meVperADCStrip * chargeAbs / pathLen;
0208 if (!shapetest || (shapetest && deDxTools::shapeSelection(clus.stripCluster()))) {
0209 dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, thit.geographicalId()));
0210 if (NSaturating > 0)
0211 NClusterSaturating++;
0212 }
0213 } else if (clus.isStrip() && thit.isMatched()) {
0214 if (!useStrip)
0215 return;
0216 const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0217 if (!matchedHit)
0218 return;
0219 const GluedGeomDet* gdet = static_cast<const GluedGeomDet*>(matchedHit->det());
0220 if (gdet == nullptr)
0221 gdet = static_cast<const GluedGeomDet*>(tkGeom->idToDet(thit.geographicalId()));
0222
0223 auto& detUnitM = *(gdet->monoDet());
0224 int NSaturating = 0;
0225 float pathLen = detUnitM.surface().bounds().thickness() / fabs(cosine);
0226 float chargeAbs = deDxTools::getCharge(&(matchedHit->monoCluster()), NSaturating, detUnitM, calibGains, m_off);
0227 float charge = meVperADCStrip * chargeAbs / pathLen;
0228 if (!shapetest || (shapetest && deDxTools::shapeSelection(matchedHit->monoCluster()))) {
0229 dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, matchedHit->monoId()));
0230 if (NSaturating > 0)
0231 NClusterSaturating++;
0232 }
0233
0234 auto& detUnitS = *(gdet->stereoDet());
0235 NSaturating = 0;
0236 pathLen = detUnitS.surface().bounds().thickness() / fabs(cosine);
0237 chargeAbs = deDxTools::getCharge(&(matchedHit->stereoCluster()), NSaturating, detUnitS, calibGains, m_off);
0238 charge = meVperADCStrip * chargeAbs / pathLen;
0239 if (!shapetest || (shapetest && deDxTools::shapeSelection(matchedHit->stereoCluster()))) {
0240 dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, matchedHit->stereoId()));
0241 if (NSaturating > 0)
0242 NClusterSaturating++;
0243 }
0244 }
0245 }
0246
0247
0248 DEFINE_FWK_MODULE(DeDxEstimatorProducer);