File indexing completed on 2024-04-06 12:31:00
0001
0002
0003
0004
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/global/EDProducer.h"
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 #include "DataFormats/Common/interface/ValueMap.h"
0014 #include "DataFormats/Common/interface/View.h"
0015 #include "DataFormats/Math/interface/deltaPhi.h"
0016
0017 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0018 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0019 #include "DataFormats/TrackReco/interface/Track.h"
0020 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0021 #include "MagneticField/Engine/interface/MagneticField.h"
0022 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0023 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0024 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0025 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0026 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0027 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0028 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0029 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0030 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0031 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0032
0033 #include <memory>
0034 #include <random>
0035
0036 #include "CLHEP/Units/SystemOfUnits.h"
0037 #include "FWCore/Utilities/interface/isFinite.h"
0038 #include "FWCore/Utilities/interface/transform.h"
0039 #include "SimTracker/TrackAssociation/interface/ResolutionModel.h"
0040
0041 class TrackTimeValueMapProducer : public edm::global::EDProducer<> {
0042 public:
0043 TrackTimeValueMapProducer(const edm::ParameterSet &);
0044 ~TrackTimeValueMapProducer() override {}
0045
0046 void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0047
0048 private:
0049
0050 const edm::EDGetTokenT<edm::View<reco::Track>> tracks_;
0051 const std::string tracksName_;
0052 const edm::EDGetTokenT<TrackingParticleCollection> trackingParticles_;
0053 const edm::EDGetTokenT<TrackingVertexCollection> trackingVertices_;
0054 const edm::EDGetTokenT<std::vector<PileupSummaryInfo>> pileupSummaryInfo_;
0055 const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> builderToken_;
0056
0057 const std::vector<edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator>> associators_;
0058
0059 const float etaMin_, etaMax_, ptMin_, pMin_, etaMaxForPtThreshold_;
0060
0061 std::vector<std::unique_ptr<const ResolutionModel>> resolutions_;
0062
0063 float extractTrackVertexTime(const TrackingParticle &, const reco::TransientTrack &) const;
0064 };
0065
0066 DEFINE_FWK_MODULE(TrackTimeValueMapProducer);
0067
0068 namespace {
0069 constexpr float m_pion = 139.57061e-3;
0070 const std::string resolution("Resolution");
0071
0072 template <typename ParticleType, typename T>
0073 void writeValueMap(edm::Event &iEvent,
0074 const edm::Handle<edm::View<ParticleType>> &handle,
0075 const std::vector<T> &values,
0076 const std::string &label) {
0077 std::unique_ptr<edm::ValueMap<T>> valMap(new edm::ValueMap<T>());
0078 typename edm::ValueMap<T>::Filler filler(*valMap);
0079 filler.insert(handle, values.begin(), values.end());
0080 filler.fill();
0081 iEvent.put(std::move(valMap), label);
0082 }
0083 }
0084
0085 TrackTimeValueMapProducer::TrackTimeValueMapProducer(const edm::ParameterSet &conf)
0086 : tracks_(consumes<edm::View<reco::Track>>(conf.getParameter<edm::InputTag>("trackSrc"))),
0087 tracksName_(conf.getParameter<edm::InputTag>("trackSrc").label()),
0088 trackingParticles_(consumes<TrackingParticleCollection>(conf.getParameter<edm::InputTag>("trackingParticleSrc"))),
0089 trackingVertices_(consumes<TrackingVertexCollection>(conf.getParameter<edm::InputTag>("trackingVertexSrc"))),
0090 pileupSummaryInfo_(
0091 consumes<std::vector<PileupSummaryInfo>>(conf.getParameter<edm::InputTag>("pileupSummaryInfo"))),
0092 builderToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0093 associators_(edm::vector_transform(
0094 conf.getParameter<std::vector<edm::InputTag>>("associators"),
0095 [this](const edm::InputTag &tag) { return this->consumes<reco::TrackToTrackingParticleAssociator>(tag); })),
0096 etaMin_(conf.getParameter<double>("etaMin")),
0097 etaMax_(conf.getParameter<double>("etaMax")),
0098 ptMin_(conf.getParameter<double>("ptMin")),
0099 pMin_(conf.getParameter<double>("pMin")),
0100 etaMaxForPtThreshold_(conf.getParameter<double>("etaMaxForPtThreshold")) {
0101
0102 const std::vector<edm::ParameterSet> &resos = conf.getParameterSetVector("resolutionModels");
0103 for (const auto &reso : resos) {
0104 const std::string &name = reso.getParameter<std::string>("modelName");
0105 resolutions_.emplace_back(ResolutionModelFactory::get()->create(name, reso));
0106
0107
0108 produces<edm::ValueMap<float>>(tracksName_ + name);
0109 produces<edm::ValueMap<float>>(tracksName_ + name + resolution);
0110 }
0111 }
0112
0113 void TrackTimeValueMapProducer::produce(edm::StreamID sid, edm::Event &evt, const edm::EventSetup &es) const {
0114
0115 std::vector<edm::Handle<reco::TrackToTrackingParticleAssociator>> associators;
0116 for (const auto &token : associators_) {
0117 associators.emplace_back();
0118 auto &back = associators.back();
0119 evt.getByToken(token, back);
0120 }
0121
0122 std::vector<float> generalTrackTimes;
0123
0124
0125 edm::Handle<edm::View<reco::Track>> TrackCollectionH;
0126 evt.getByToken(tracks_, TrackCollectionH);
0127 const edm::View<reco::Track> &TrackCollection = *TrackCollectionH;
0128
0129
0130 edm::Handle<TrackingParticleCollection> TPCollectionH;
0131 evt.getByToken(trackingParticles_, TPCollectionH);
0132
0133 edm::Handle<std::vector<PileupSummaryInfo>> pileupSummaryH;
0134 evt.getByToken(pileupSummaryInfo_, pileupSummaryH);
0135
0136
0137 auto const &builder = es.getData(builderToken_);
0138
0139
0140 std::vector<reco::RecoToSimCollection> associatedTracks;
0141 associatedTracks.reserve(associators.size());
0142 for (const auto &associator : associators) {
0143 associatedTracks.emplace_back(associator->associateRecoToSim(TrackCollectionH, TPCollectionH));
0144 }
0145
0146 double sumSimTime = 0.;
0147 double sumSimTimeSq = 0.;
0148 int nsim = 0;
0149 for (const PileupSummaryInfo &puinfo : *pileupSummaryH) {
0150 if (puinfo.getBunchCrossing() == 0) {
0151 for (const float &time : puinfo.getPU_times()) {
0152 double simtime = time;
0153 sumSimTime += simtime;
0154 sumSimTimeSq += simtime * simtime;
0155 ++nsim;
0156 }
0157 break;
0158 }
0159 }
0160
0161 double meanSimTime = sumSimTime / double(nsim);
0162 double varSimTime = sumSimTimeSq / double(nsim) - meanSimTime * meanSimTime;
0163 double rmsSimTime = std::sqrt(std::max(0.1 * 0.1, varSimTime));
0164 std::normal_distribution<float> gausSimTime(meanSimTime, rmsSimTime);
0165
0166
0167 unsigned int runNum_uint = static_cast<unsigned int>(evt.id().run());
0168 unsigned int lumiNum_uint = static_cast<unsigned int>(evt.id().luminosityBlock());
0169 unsigned int evNum_uint = static_cast<unsigned int>(evt.id().event());
0170 unsigned int tkChi2_uint = uint32_t(TrackCollection.empty() ? 0 : TrackCollection.refAt(0)->chi2() / 0.01);
0171 std::uint32_t seed = tkChi2_uint + (lumiNum_uint << 10) + (runNum_uint << 20) + evNum_uint;
0172 std::mt19937 rng(seed);
0173
0174 for (unsigned itk = 0; itk < TrackCollection.size(); ++itk) {
0175 const auto tkref = TrackCollection.refAt(itk);
0176 reco::RecoToSimCollection::const_iterator track_tps = associatedTracks.back().end();
0177
0178 for (const auto &association : associatedTracks) {
0179 track_tps = association.find(tkref);
0180 if (track_tps != association.end())
0181 break;
0182 }
0183
0184 if (track_tps != associatedTracks.back().end() && track_tps->val.size() == 1) {
0185 reco::TransientTrack tt = builder.build(*tkref);
0186 float time = extractTrackVertexTime(*track_tps->val[0].first, tt);
0187 generalTrackTimes.push_back(time);
0188 } else {
0189 float rndtime = gausSimTime(rng);
0190 generalTrackTimes.push_back(rndtime);
0191 if (track_tps != associatedTracks.back().end() && track_tps->val.size() > 1) {
0192 LogDebug("TooManyTracks") << "track matched to " << track_tps->val.size() << " tracking particles!"
0193 << std::endl;
0194 }
0195 }
0196 }
0197
0198 for (const auto &reso : resolutions_) {
0199 const std::string &name = reso->name();
0200 std::vector<float> times, resos;
0201
0202 times.reserve(TrackCollection.size());
0203 resos.reserve(TrackCollection.size());
0204
0205 for (unsigned i = 0; i < TrackCollection.size(); ++i) {
0206 const reco::Track &tk = TrackCollection[i];
0207 const float absEta = std::abs(tk.eta());
0208 bool inAcceptance = absEta < etaMax_ && absEta >= etaMin_ && tk.p() > pMin_ &&
0209 (absEta > etaMaxForPtThreshold_ || tk.pt() > ptMin_);
0210 if (inAcceptance) {
0211 const float resolution = reso->getTimeResolution(tk);
0212 std::normal_distribution<float> gausGeneralTime(generalTrackTimes[i], resolution);
0213 times.push_back(gausGeneralTime(rng));
0214 resos.push_back(resolution);
0215 } else {
0216 times.push_back(0.0f);
0217 resos.push_back(-1.);
0218 }
0219 }
0220
0221 writeValueMap(evt, TrackCollectionH, times, tracksName_ + name);
0222 writeValueMap(evt, TrackCollectionH, resos, tracksName_ + name + resolution);
0223 }
0224 }
0225
0226 float TrackTimeValueMapProducer::extractTrackVertexTime(const TrackingParticle &tp,
0227 const reco::TransientTrack &tt) const {
0228 int pdgid = tp.pdgId();
0229 const auto &tvertex = tp.parentVertex();
0230 math::XYZTLorentzVectorD result = tvertex->position();
0231
0232
0233 if (tvertex->nSourceTracks() && tvertex->sourceTracks()[0]->pdgId() == pdgid) {
0234 auto pvertex = tvertex->sourceTracks()[0]->parentVertex();
0235 result = pvertex->position();
0236 while (pvertex->nSourceTracks() && pvertex->sourceTracks()[0]->pdgId() == pdgid) {
0237 pvertex = pvertex->sourceTracks()[0]->parentVertex();
0238 result = pvertex->position();
0239 }
0240 }
0241
0242 float time = result.T() * CLHEP::second;
0243
0244 GlobalPoint result_pos(result.x(), result.y(), result.z());
0245 const auto &tkstate = tt.trajectoryStateClosestToPoint(result_pos);
0246 float tkphi = tkstate.momentum().phi();
0247 float tkz = tkstate.position().z();
0248 float dphi = reco::deltaPhi(tkphi, tt.track().phi());
0249 float dz = tkz - tt.track().vz();
0250
0251 float radius = 100. * tt.track().pt() / (0.3 * tt.field()->inTesla(GlobalPoint(0, 0, 0)).z());
0252 float pathlengthrphi = tt.track().charge() * dphi * radius;
0253
0254 float pathlength = std::sqrt(pathlengthrphi * pathlengthrphi + dz * dz);
0255 float p = tt.track().p();
0256
0257 float speed = std::sqrt(1. / (1. + m_pion / p)) * CLHEP::c_light / CLHEP::cm;
0258 float dt = pathlength / speed;
0259
0260 return time - dt;
0261 }