TrackAssociation/plugins/TrackTimeValueMapProducer.cc

0001 // This producer assigns vertex times (with a specified resolution) to tracks.
0002 // The times are produced as valuemaps associated to tracks, so the track
0003 // dataformat doesn't need to be modified.
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   // inputs
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   // tracking particle associators by order of preference
0057   const std::vector<edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator>> associators_;
0058   // eta bounds
0059   const float etaMin_, etaMax_, ptMin_, pMin_, etaMaxForPtThreshold_;
0060   // options
0061   std::vector<std::unique_ptr<const ResolutionModel>> resolutions_;
0062   // functions
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 }  // namespace
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   // setup resolution models
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     // times and time resolutions for general tracks
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   // get sim track associators
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   // get track collections
0125   edm::Handle<edm::View<reco::Track>> TrackCollectionH;
0126   evt.getByToken(tracks_, TrackCollectionH);
0127   const edm::View<reco::Track> &TrackCollection = *TrackCollectionH;
0128
0129   // get tracking particle collections
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   // transient track builder
0137   auto const &builder = es.getData(builderToken_);
0138
0139   // associate the reco tracks / gsf Tracks
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   // get event-based seed for RNG
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   // account for secondary vertices...
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   // correct for time of flight from track reference position
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;  // speed in cm/ns
0258   float dt = pathlength / speed;
0259
0260   return time - dt;
0261 }