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