File indexing completed on 2024-04-06 12:11:20
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 #include <memory>
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/stream/EDProducer.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/Framework/interface/ConsumesCollector.h"
0030
0031 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0032 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0033 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0034 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0035 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0036
0037 #include "DataFormats/Common/interface/ValueMap.h"
0038 #include "DataFormats/TrackReco/interface/DeDxData.h"
0039 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
0040 #include "DataFormats/TrackReco/interface/DeDxHit.h"
0041 #include "DataFormats/TrackReco/interface/Track.h"
0042 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0043
0044 #include "RecoTracker/DeDx/interface/BaseDeDxEstimator.h"
0045 #include "RecoTracker/DeDx/interface/GenericAverageDeDxEstimator.h"
0046 #include "RecoTracker/DeDx/interface/TruncatedAverageDeDxEstimator.h"
0047 #include "RecoTracker/DeDx/interface/MedianDeDxEstimator.h"
0048 #include "RecoTracker/DeDx/interface/UnbinnedFitDeDxEstimator.h"
0049 #include "RecoTracker/DeDx/interface/ProductDeDxDiscriminator.h"
0050 #include "RecoTracker/DeDx/interface/SmirnovDeDxDiscriminator.h"
0051 #include "RecoTracker/DeDx/interface/ASmirnovDeDxDiscriminator.h"
0052 #include "RecoTracker/DeDx/interface/BTagLikeDeDxDiscriminator.h"
0053
0054 #include "RecoTracker/DeDx/interface/DeDxTools.h"
0055
0056 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0057
0058 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0059 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0060 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0061 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHitCollection.h"
0062
0063
0064
0065
0066
0067 class FastTrackDeDxProducer : public edm::stream::EDProducer<> {
0068 public:
0069 explicit FastTrackDeDxProducer(const edm::ParameterSet&);
0070 ~FastTrackDeDxProducer() override = default;
0071 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0072
0073 private:
0074 void beginRun(edm::Run const&, const edm::EventSetup&) override;
0075 void produce(edm::Event&, const edm::EventSetup&) override;
0076
0077 void makeCalibrationMap(const TrackerGeometry& tkGeom);
0078 void processHit(const FastTrackerRecHit& recHit,
0079 float trackMomentum,
0080 float& cosine,
0081 reco::DeDxHitCollection& dedxHits,
0082 int& NClusterSaturating);
0083
0084
0085
0086 std::unique_ptr<BaseDeDxEstimator> m_estimator;
0087
0088 edm::EDGetTokenT<reco::TrackCollection> m_tracksTag;
0089
0090 float meVperADCPixel;
0091 float meVperADCStrip;
0092
0093 unsigned int MaxNrStrips;
0094
0095 std::string m_calibrationPath;
0096
0097 std::vector<std::vector<float>> calibGains;
0098 unsigned int offsetDU_;
0099
0100 edm::EDGetTokenT<edm::PSimHitContainer> simHitsToken;
0101 edm::EDGetTokenT<FastTrackerRecHitRefCollection> simHit2RecHitMapToken;
0102 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken;
0103
0104 bool usePixel;
0105 bool useStrip;
0106 bool useCalibration;
0107 bool shapetest;
0108 bool convertFromGeV2MeV;
0109 bool nothick;
0110 };
0111
0112 using namespace reco;
0113 using namespace std;
0114 using namespace edm;
0115
0116 void FastTrackDeDxProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0117 edm::ParameterSetDescription desc;
0118 desc.add<string>("estimator", "generic");
0119 desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0120 desc.add<bool>("UsePixel", false);
0121 desc.add<bool>("UseStrip", true);
0122 desc.add<double>("MeVperADCStrip", 3.61e-06 * 265);
0123 desc.add<double>("MeVperADCPixel", 3.61e-06);
0124 desc.add<bool>("ShapeTest", true);
0125 desc.add<bool>("UseCalibration", false);
0126 desc.add<string>("calibrationPath", "");
0127 desc.add<string>("Record", "SiStripDeDxMip_3D_Rcd");
0128 desc.add<string>("ProbabilityMode", "Accumulation");
0129 desc.add<double>("fraction", 0.4);
0130 desc.add<double>("exponent", -2.0);
0131 desc.add<bool>("convertFromGeV2MeV", true);
0132 desc.add<bool>("nothick", false);
0133 desc.add<edm::InputTag>("simHits");
0134 desc.add<edm::InputTag>("simHit2RecHitMap");
0135 descriptions.add("FastTrackDeDxProducer", desc);
0136 }
0137
0138 FastTrackDeDxProducer::FastTrackDeDxProducer(const edm::ParameterSet& iConfig)
0139 : simHitsToken(consumes<edm::PSimHitContainer>(iConfig.getParameter<edm::InputTag>("simHits"))),
0140 simHit2RecHitMapToken(
0141 consumes<FastTrackerRecHitRefCollection>(iConfig.getParameter<edm::InputTag>("simHit2RecHitMap"))) {
0142 produces<ValueMap<DeDxData>>();
0143
0144 auto cCollector = consumesCollector();
0145 string estimatorName = iConfig.getParameter<string>("estimator");
0146 if (estimatorName == "median")
0147 m_estimator = std::unique_ptr<BaseDeDxEstimator>(new MedianDeDxEstimator(iConfig));
0148 else if (estimatorName == "generic")
0149 m_estimator = std::unique_ptr<BaseDeDxEstimator>(new GenericAverageDeDxEstimator(iConfig));
0150 else if (estimatorName == "truncated")
0151 m_estimator = std::unique_ptr<BaseDeDxEstimator>(new TruncatedAverageDeDxEstimator(iConfig));
0152
0153 else if (estimatorName == "productDiscrim")
0154 m_estimator = std::unique_ptr<BaseDeDxEstimator>(new ProductDeDxDiscriminator(iConfig, cCollector));
0155 else if (estimatorName == "btagDiscrim")
0156 m_estimator = std::unique_ptr<BaseDeDxEstimator>(new BTagLikeDeDxDiscriminator(iConfig, cCollector));
0157 else if (estimatorName == "smirnovDiscrim")
0158 m_estimator = std::unique_ptr<BaseDeDxEstimator>(new SmirnovDeDxDiscriminator(iConfig, cCollector));
0159 else if (estimatorName == "asmirnovDiscrim")
0160 m_estimator = std::unique_ptr<BaseDeDxEstimator>(new ASmirnovDeDxDiscriminator(iConfig, cCollector));
0161 else
0162 throw cms::Exception("fastsim::SimplifiedGeometry::FastTrackDeDxProducer.cc") << " estimator name does not exist";
0163
0164
0165
0166
0167 m_tracksTag = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0168
0169 usePixel = iConfig.getParameter<bool>("UsePixel");
0170 useStrip = iConfig.getParameter<bool>("UseStrip");
0171 meVperADCPixel = iConfig.getParameter<double>("MeVperADCPixel");
0172 meVperADCStrip = iConfig.getParameter<double>("MeVperADCStrip");
0173
0174 shapetest = iConfig.getParameter<bool>("ShapeTest");
0175 useCalibration = iConfig.getParameter<bool>("UseCalibration");
0176 m_calibrationPath = iConfig.getParameter<string>("calibrationPath");
0177
0178 convertFromGeV2MeV = iConfig.getParameter<bool>("convertFromGeV2MeV");
0179 nothick = iConfig.getParameter<bool>("nothick");
0180
0181 if (!usePixel && !useStrip)
0182 throw cms::Exception("fastsim::SimplifiedGeometry::FastTrackDeDxProducer.cc")
0183 << " neither pixel hits nor strips hits will be used to compute de/dx";
0184
0185 if (useCalibration) {
0186 tkGeomToken = esConsumes<edm::Transition::BeginRun>();
0187 }
0188 }
0189
0190
0191 void FastTrackDeDxProducer::beginRun(edm::Run const& run, const edm::EventSetup& iSetup) {
0192 if (useCalibration && calibGains.empty()) {
0193 auto const& tkGeom = iSetup.getData(tkGeomToken);
0194 offsetDU_ = tkGeom.offsetDU(GeomDetEnumerators::PixelBarrel);
0195
0196 deDxTools::makeCalibrationMap(m_calibrationPath, tkGeom, calibGains, offsetDU_);
0197 }
0198
0199 m_estimator->beginRun(run, iSetup);
0200 }
0201
0202 void FastTrackDeDxProducer::produce(edm::Event& iEvent, const edm::EventSetup&) {
0203 auto trackDeDxEstimateAssociation = std::make_unique<ValueMap<DeDxData>>();
0204 ValueMap<DeDxData>::Filler filler(*trackDeDxEstimateAssociation);
0205
0206 edm::Handle<reco::TrackCollection> trackCollectionHandle;
0207 iEvent.getByToken(m_tracksTag, trackCollectionHandle);
0208 const auto& trackCollection = *trackCollectionHandle;
0209 std::vector<DeDxData> dedxEstimate(trackCollection.size());
0210
0211 for (unsigned int j = 0; j < trackCollection.size(); j++) {
0212 const reco::TrackRef track = reco::TrackRef(trackCollectionHandle.product(), j);
0213
0214 int NClusterSaturating = 0;
0215 DeDxHitCollection dedxHits;
0216
0217 auto const& trajParams = track->extra()->trajParams();
0218 assert(trajParams.size() == track->recHitsSize());
0219
0220 auto hb = track->recHitsBegin();
0221 dedxHits.reserve(track->recHitsSize() / 2);
0222
0223 for (unsigned int h = 0; h < track->recHitsSize(); h++) {
0224 const FastTrackerRecHit& recHit = static_cast<const FastTrackerRecHit&>(*(*(hb + h)));
0225 if (!recHit.isValid())
0226 continue;
0227 auto trackDirection = trajParams[h].direction();
0228 float cosine = trackDirection.z() / trackDirection.mag();
0229 processHit(recHit, track->p(), cosine, dedxHits, NClusterSaturating);
0230 }
0231
0232 sort(dedxHits.begin(), dedxHits.end(), less<DeDxHit>());
0233 std::pair<float, float> val_and_error = m_estimator->dedx(dedxHits);
0234
0235
0236 val_and_error.second = NClusterSaturating;
0237 dedxEstimate[j] = DeDxData(val_and_error.first, val_and_error.second, dedxHits.size());
0238 }
0239
0240 filler.insert(trackCollectionHandle, dedxEstimate.begin(), dedxEstimate.end());
0241
0242 filler.fill();
0243 iEvent.put(std::move(trackDeDxEstimateAssociation));
0244 }
0245
0246 void FastTrackDeDxProducer::processHit(const FastTrackerRecHit& recHit,
0247 float trackMomentum,
0248 float& cosine,
0249 reco::DeDxHitCollection& dedxHits,
0250 int& NClusterSaturating) {
0251 if (!recHit.isValid())
0252 return;
0253
0254 auto const& thit = static_cast<BaseTrackerRecHit const&>(recHit);
0255 if (!thit.isValid())
0256 return;
0257 if (!thit.hasPositionAndError())
0258 return;
0259
0260 if (recHit.isPixel()) {
0261 if (!usePixel)
0262 return;
0263
0264 auto& detUnit = *(recHit.detUnit());
0265 float pathLen = detUnit.surface().bounds().thickness() / fabs(cosine);
0266 if (nothick)
0267 pathLen = 1.0;
0268 float charge = recHit.energyLoss() / pathLen;
0269 if (convertFromGeV2MeV)
0270 charge *= 1000;
0271 dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, recHit.geographicalId()));
0272 } else if (!recHit.isPixel()) {
0273 if (!useStrip)
0274 return;
0275 auto& detUnit = *(recHit.detUnit());
0276 float pathLen = detUnit.surface().bounds().thickness() / fabs(cosine);
0277 if (nothick)
0278 pathLen = 1.0;
0279 float dedxOfRecHit = recHit.energyLoss() / pathLen;
0280 if (convertFromGeV2MeV)
0281 dedxOfRecHit *= 1000;
0282 if (!shapetest) {
0283 dedxHits.push_back(DeDxHit(dedxOfRecHit, trackMomentum, pathLen, recHit.geographicalId()));
0284 }
0285 }
0286 }
0287
0288
0289 DEFINE_FWK_MODULE(FastTrackDeDxProducer);