Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:31:00

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 }